Implement frexp

19 files changed, 99 insertions, 1 deletions

diff --git a/floattypes.h b/floattypes.h
index 8e70b27..426b23a 100644
--- a/floattypes.h
+++ b/floattypes.h
@@ -46,6 +46,7 @@ namespace std {
   float fma(float x, float y, float z) { return ::fmaf(x, y, z); }
   float fmax(float x, float y) { return ::fmaxf(x, y); }
   float fmin(float x, float y) { return ::fminf(x, y); }
+  float frexp(float x, int* r) { return ::frexp(x, r); }
   float hypot(float x, float y) { return ::hypotf(x, y); }
   int ilogb(float x) { return ::ilogbf(x); }
   float log1p(float x) { return ::log1pf(x); }
@@ -70,6 +71,7 @@ namespace std {
   double fma(double x, double y, double z) { return ::fma(x, y, z); }
   double fmax(double x, double y) { return ::fmax(x, y); }
   double fmin(double x, double y) { return ::fmin(x, y); }
+  double frexp(double x, int* r) { return ::frexp(x, r); }
   double hypot(double x, double y) { return ::hypot(x, y); }
   int ilogb(double x) { return ::ilogb(x); }
   double log1p(double x) { return ::log1p(x); }
diff --git a/mathfuncs_base.h b/mathfuncs_base.h
index 20dba73..06cb590 100644
--- a/mathfuncs_base.h
+++ b/mathfuncs_base.h
@@ -67,6 +67,7 @@ namespace vecmathlib {
     static realvec_t vml_fma(realvec_t x, realvec_t y, realvec_t z);
     static realvec_t vml_fmax(realvec_t x, realvec_t y);
     static realvec_t vml_fmin(realvec_t x, realvec_t y);
+    static realvec_t vml_frexp(realvec_t x, intvec_t& r);
     static intvec_t vml_ilogb(realvec_t x);
     static boolvec_t vml_ieee_isfinite(realvec_t x);
     static boolvec_t vml_ieee_isinf(realvec_t x);
diff --git a/mathfuncs_fabs.h b/mathfuncs_fabs.h
index d6e78f0..07b5664 100644
--- a/mathfuncs_fabs.h
+++ b/mathfuncs_fabs.h
@@ -51,6 +51,28 @@ namespace vecmathlib {
   }
   
   template<typename realvec_t>
+  realvec_t mathfuncs<realvec_t>::vml_frexp(realvec_t x,
+                                            typename realvec_t::intvec_t& ir)
+  {
+    intvec_t e = lsr(as_int(x) & IV(FP::exponent_mask), FP::mantissa_bits);
+    ir = e - IV(FP::exponent_offset - 1);
+    ir = ifthen(convert_bool(e), ir, IV(std::numeric_limits<int_t>::min()));
+#if defined VML_HAVE_INF
+    ir = ifthen(isinf(x), IV(std::numeric_limits<int_t>::max()), ir);
+#endif
+#if defined VML_HAVE_NAN
+    ir = ifthen(isnan(x), IV(std::numeric_limits<int_t>::min()), ir);
+#endif
+    realvec_t r =
+      as_float((as_int(x) & IV(FP::signbit_mask | FP::mantissa_mask)) |
+               IV(FP::as_int(R(0.5)) & FP::exponent_mask));
+    boolvec_t iszero = x == RV(0.0);
+    ir = ifthen(iszero, IV(I(0)), ir);
+    r = ifthen(iszero, copysign(RV(R(0.0)), r), r);
+    return r;
+  }
+  
+  template<typename realvec_t>
   typename realvec_t::intvec_t mathfuncs<realvec_t>::vml_ilogb(realvec_t x)
   {
     intvec_t e = lsr(as_int(x) & IV(FP::exponent_mask), FP::mantissa_bits);
diff --git a/test.cc b/test.cc
index 960db9a..2ee0fb4 100644
--- a/test.cc
+++ b/test.cc
@@ -1024,6 +1024,30 @@ struct vecmathlib_test {
   }
   
   // Change signature: "int" -> "int_t"
+  static real_t local_frexp0(real_t x)
+  {
+    int r;
+    return std::frexp(x, &r);
+  }
+  static int_t local_frexp1(real_t x)
+  {
+    int r;
+    std::frexp(x, &r);
+    if (r==FP_ILOGB0) r=std::numeric_limits<int_t>::min();
+    if (r==FP_ILOGBNAN) r=std::numeric_limits<int_t>::max();
+    return r;
+  }
+  static realvec_t local_vfrexp0(realvec_t x)
+  {
+    intvec_t r;
+    return vecmathlib::frexp(x, r);
+  }
+  static intvec_t local_vfrexp1(realvec_t x)
+  {
+    intvec_t r;
+    vecmathlib::frexp(x, r);
+    return r;
+  }
   static int_t local_ilogb(real_t x)
   {
     int r = std::ilogb(x);
@@ -1034,7 +1058,7 @@ struct vecmathlib_test {
   static real_t local_ldexp(real_t x, int_t n) { return ldexp(x, n); }
   static void test_fabs()
   {
-    cout << "   testing + - + - * == != < <= > >= copysign fabs fdim fma fmax fmin ilogb isfinite isinf isnan isnormal ldexp nextafter signbit...\n" << flush;
+    cout << "   testing + - + - * == != < <= > >= copysign fabs fdim fma fmax fmin frexp ilogb isfinite isinf isnan isnormal ldexp nextafter signbit...\n" << flush;
     
     const real_t eps = FP::epsilon();
     const real_t int_min = R(std::numeric_limits<int_t>::min());
@@ -1126,6 +1150,8 @@ struct vecmathlib_test {
                            x, y, z, R(2.0)*accuracy());
       check_real<RV,RV>("fmax", std::fmax, vecmathlib::fmax, x, y, 0.0);
       check_real<RV,RV>("fmin", std::fmin, vecmathlib::fmin, x, y, 0.0);
+      check_real<RV>("frexp0", local_frexp0, local_vfrexp0, x, 0.0);
+      check_int<RV>("frexp1", local_frexp1, local_vfrexp1, x);
       check_int<RV>("ilogb",
                     local_ilogb, (intvec_t(*)(realvec_t))vecmathlib::ilogb, x);
 #if defined VML_HAVE_INF && defined VML_HAVE_NAN
diff --git a/vec_altivec_float4.h b/vec_altivec_float4.h
index 813141e..c5eb28b 100644
--- a/vec_altivec_float4.h
+++ b/vec_altivec_float4.h
@@ -445,6 +445,7 @@ namespace vecmathlib {
     realvec fmax(realvec y) const { return vec_max(v, y.v); }
     realvec fmin(realvec y) const { return vec_min(v, y.v); }
     realvec fmod(realvec y) const { return MF::vml_fmod(*this, y); }
+    realvec frexp(intvec_t& r) const { return MF::vml_frexp(*this, r); }
     realvec hypot(realvec y) const { return MF::vml_hypot(*this, y); }
     intvec_t ilogb() const { return MF::vml_ilogb(*this); }
     boolvec_t isfinite() const { return MF::vml_isfinite(*this); }
diff --git a/vec_avx_double4.h b/vec_avx_double4.h
index 15f7edf..9102e67 100644
--- a/vec_avx_double4.h
+++ b/vec_avx_double4.h
@@ -564,6 +564,7 @@ namespace vecmathlib {
     realvec fmax(realvec y) const { return _mm256_max_pd(v, y.v); }
     realvec fmin(realvec y) const { return _mm256_min_pd(v, y.v); }
     realvec fmod(realvec y) const { return MF::vml_fmod(*this, y); }
+    realvec frexp(intvec_t& r) const { return MF::vml_frexp(*this, r); }
     realvec hypot(realvec y) const { return MF::vml_hypot(*this, y); }
     intvec_t ilogb() const { return MF::vml_ilogb(*this); }
     boolvec_t isfinite() const { return MF::vml_isfinite(*this); }
diff --git a/vec_avx_float8.h b/vec_avx_float8.h
index e7e1187..6431227 100644
--- a/vec_avx_float8.h
+++ b/vec_avx_float8.h
@@ -556,6 +556,7 @@ namespace vecmathlib {
     realvec fmax(realvec y) const { return _mm256_max_ps(v, y.v); }
     realvec fmin(realvec y) const { return _mm256_min_ps(v, y.v); }
     realvec fmod(realvec y) const { return MF::vml_fmod(*this, y); }
+    realvec frexp(intvec_t& r) const { return MF::vml_frexp(*this, r); }
     realvec hypot(realvec y) const { return MF::vml_hypot(*this, y); }
     intvec_t ilogb() const { return MF::vml_ilogb(*this); }
     boolvec_t isfinite() const { return MF::vml_isfinite(*this); }
diff --git a/vec_base.h b/vec_base.h
index ee7864a..104ac97 100644
--- a/vec_base.h
+++ b/vec_base.h
@@ -358,6 +358,13 @@ namespace vecmathlib {
   }
   
   template<typename real_t, int size>
+  inline realvec<real_t, size> frexp(realvec<real_t, size> x,
+                                     intvec<real_t, size>& r)
+  {
+    return x.frexp(r);
+  }
+  
+  template<typename real_t, int size>
   inline realvec<real_t, size> hypot(realvec<real_t, size> x,
                                      realvec<real_t, size> y)
   {
diff --git a/vec_builtin.h b/vec_builtin.h
index a240a6d..f7db8ca 100644
--- a/vec_builtin.h
+++ b/vec_builtin.h
@@ -532,6 +532,7 @@ namespace vecmathlib {
     realvec_t fmax(realvec_t y) const { return MF::vml_fmax(*this, y); }
     realvec_t fmin(realvec_t y) const { return MF::vml_fmin(*this, y); }
     realvec_t fmod(realvec_t y) const { return MF::vml_fmod(*this, y); }
+    realvec frexp(intvec_t& r) const { return MF::vml_frexp(*this, r); }
     realvec_t hypot(realvec_t y) const { return MF::vml_hypot(*this, y); }
     intvec_t ilogb() const { return MF::vml_ilogb(*this); }
     boolvec_t isfinite() const { return MF::vml_isfinite(*this); }
diff --git a/vec_neon_float2.h b/vec_neon_float2.h
index 0030c88..ba14d86 100644
--- a/vec_neon_float2.h
+++ b/vec_neon_float2.h
@@ -459,6 +459,7 @@ namespace vecmathlib {
     realvec fmax(realvec y) const { return vmax_f32(v, y.v); }
     realvec fmin(realvec y) const { return vmin_f32(v, y.v); }
     realvec fmod(realvec y) const { return MF::vml_fmod(*this, y); }
+    realvec frexp(intvec_t& r) const { return MF::vml_frexp(*this, r); }
     realvec hypot(realvec y) const { return MF::vml_hypot(*this, y); }
     intvec_t ilogb() const { return MF::vml_ilogb(*this); }
     boolvec_t isfinite() const { return MF::vml_isfinite(*this); }
diff --git a/vec_neon_float4.h b/vec_neon_float4.h
index cf72a6d..2951420 100644
--- a/vec_neon_float4.h
+++ b/vec_neon_float4.h
@@ -473,6 +473,7 @@ namespace vecmathlib {
     realvec fmax(realvec y) const { return vmaxq_f32(v, y.v); }
     realvec fmin(realvec y) const { return vminq_f32(v, y.v); }
     realvec fmod(realvec y) const { return MF::vml_fmod(*this, y); }
+    realvec frexp(intvec_t& r) const { return MF::vml_frexp(*this, r); }
     realvec hypot(realvec y) const { return MF::vml_hypot(*this, y); }
     intvec_t ilogb() const { return MF::vml_ilogb(*this); }
     boolvec_t isfinite() const { return MF::vml_isfinite(*this); }
diff --git a/vec_pseudo.h b/vec_pseudo.h
index 0764d94..c3658c0 100644
--- a/vec_pseudo.h
+++ b/vec_pseudo.h
@@ -746,6 +746,19 @@ namespace vecmathlib {
     realpseudovec fmax(realpseudovec y) const { return map(std::fmax, y); }
     realpseudovec fmin(realpseudovec y) const { return map(std::fmin, y); }
     realpseudovec fmod(realpseudovec y) const { return map(std::fmod, y); }
+    realpseudovec frexp(intvec_t& ires) const
+    {
+      realvec_t res;
+      for (int d=0; d<size; ++d) {
+        int ir;
+        real_t r = std::frexp(v[d], &ir);
+        if (ir == FP_ILOGB0) ir = std::numeric_limits<int_t>::min();
+        else if (ir == FP_ILOGBNAN) ir = std::numeric_limits<int_t>::max();
+        res.v[d] = r;
+        ires.v[d] = ir;
+      }
+      return res;
+    }
     realpseudovec hypot(realpseudovec y) const { return map(std::hypot, y); }
     intvec_t ilogb() const
     {
@@ -1203,6 +1216,13 @@ namespace vecmathlib {
   }
   
   template<typename real_t, int size>
+  inline realpseudovec<real_t, size> frexp(realpseudovec<real_t, size> x,
+                                           intpseudovec<real_t, size>& r)
+  {
+    return x.frexp(r);
+  }
+  
+  template<typename real_t, int size>
   inline realpseudovec<real_t, size> hypot(realpseudovec<real_t, size> x,
                                            realpseudovec<real_t, size> y)
   {
diff --git a/vec_qpx_double4.h b/vec_qpx_double4.h
index e7a1c05..afcbcbb 100644
--- a/vec_qpx_double4.h
+++ b/vec_qpx_double4.h
@@ -560,6 +560,7 @@ namespace vecmathlib {
     realvec fmax(realvec y) const { return MF::vml_fmax(v, y.v); }
     realvec fmin(realvec y) const { return MF::vml_fmin(v, y.v); }
     realvec fmod(realvec y) const { return MF::vml_fmod(*this, y); }
+    realvec frexp(intvec_t& r) const { return MF::vml_frexp(*this, r); }
     realvec hypot(realvec y) const { return hypotd4(v, y.v); }
     intvec_t ilogb() const
     {
diff --git a/vec_sse_double1.h b/vec_sse_double1.h
index d4626dd..a30cd2d 100644
--- a/vec_sse_double1.h
+++ b/vec_sse_double1.h
@@ -421,6 +421,7 @@ namespace vecmathlib {
       return to_double(_mm_min_sd(from_double(v), from_double(y.v)));
     }
     realvec fmod(realvec y) const { return std::fmod(v, y.v); }
+    realvec frexp(intvec_t& r) const { return MF::vml_frexp(*this, r); }
     realvec hypot(realvec y) const { return MF::vml_hypot(*this, y); }
     intvec_t ilogb() const
     {
diff --git a/vec_sse_double2.h b/vec_sse_double2.h
index 2cbbd38..ee689cd 100644
--- a/vec_sse_double2.h
+++ b/vec_sse_double2.h
@@ -551,6 +551,7 @@ namespace vecmathlib {
     realvec fmax(realvec y) const { return _mm_max_pd(v, y.v); }
     realvec fmin(realvec y) const { return _mm_min_pd(v, y.v); }
     realvec fmod(realvec y) const { return MF::vml_fmod(*this, y); }
+    realvec frexp(intvec_t& r) const { return MF::vml_frexp(*this, r); }
     realvec hypot(realvec y) const { return MF::vml_hypot(*this, y); }
     intvec_t ilogb() const { return MF::vml_ilogb(*this); }
     boolvec_t isfinite() const { return MF::vml_isfinite(*this); }
diff --git a/vec_sse_float1.h b/vec_sse_float1.h
index 56ece24..12ab1ae 100644
--- a/vec_sse_float1.h
+++ b/vec_sse_float1.h
@@ -421,6 +421,7 @@ namespace vecmathlib {
       return to_float(_mm_min_ss(from_float(v), from_float(y.v)));
     }
     realvec fmod(realvec y) const { return std::fmod(v, y.v); }
+    realvec frexp(intvec_t& r) const { return MF::vml_frexp(*this, r); }
     realvec hypot(realvec y) const { return MF::vml_hypot(*this, y); }
     intvec_t ilogb() const
     {
diff --git a/vec_sse_float4.h b/vec_sse_float4.h
index 5259cb2..2163e09 100644
--- a/vec_sse_float4.h
+++ b/vec_sse_float4.h
@@ -545,6 +545,7 @@ namespace vecmathlib {
     realvec fmax(realvec y) const { return _mm_max_ps(v, y.v); }
     realvec fmin(realvec y) const { return _mm_min_ps(v, y.v); }
     realvec fmod(realvec y) const { return MF::vml_fmod(*this, y); }
+    realvec frexp(intvec_t& r) const { return MF::vml_frexp(*this, r); }
     realvec hypot(realvec y) const { return MF::vml_hypot(*this, y); }
     intvec_t ilogb() const { return MF::vml_ilogb(*this); }
     boolvec_t isfinite() const { return MF::vml_isfinite(*this); }
diff --git a/vec_test.h b/vec_test.h
index f0bbbea..2d05b67 100644
--- a/vec_test.h
+++ b/vec_test.h
@@ -688,6 +688,7 @@ namespace vecmathlib {
     realtestvec fmax(realtestvec y) const { return MF::vml_fmax(*this, y); }
     realtestvec fmin(realtestvec y) const { return MF::vml_fmin(*this, y); }
     realtestvec fmod(realtestvec y) const { return MF::vml_fmod(*this, y); }
+    realtestvec frexp(intvec_t& r) const { return MF::vml_frexp(*this, r); }
     realtestvec hypot(realtestvec y) const { return MF::vml_hypot(*this, y); }
     intvec_t ilogb() const { return MF::vml_ilogb(*this); }
     boolvec_t isfinite() const { return MF::vml_isfinite(*this); }
@@ -1118,6 +1119,13 @@ namespace vecmathlib {
   }
   
   template<typename real_t, int size>
+  inline realtestvec<real_t, size> frexp(realtestvec<real_t, size> x,
+                                         inttestvec<real_t, size>& r)
+  {
+    return x.frexp(r);
+  }
+  
+  template<typename real_t, int size>
   inline realtestvec<real_t, size> hypot(realtestvec<real_t, size> x,
                                          realtestvec<real_t, size> y)
   {
diff --git a/vec_vsx_double2.h b/vec_vsx_double2.h
index 024056d..562969a 100644
--- a/vec_vsx_double2.h
+++ b/vec_vsx_double2.h
@@ -560,6 +560,7 @@ namespace vecmathlib {
     realvec fmax(realvec y) const { return vec_max(v, y.v); }
     realvec fmin(realvec y) const { return vec_min(v, y.v); }
     realvec fmod(realvec y) const { return MF::vml_fmod(*this, y); }
+    realvec frexp(intvec_t& r) const { return MF::vml_frexp(*this, r); }
     realvec hypot(realvec y) const { return MF::vml_hypot(*this, y); }
     intvec_t ilogb() const { return MF::vml_ilogb(*this); }
     boolvec_t isfinite() const { return MF::vml_isfinite(*this); }