Implement ceil floor fmod pow remainder round

10 files changed, 146 insertions, 29 deletions

diff --git a/mathfuncs.h b/mathfuncs.h
index f693150..697a868 100644
--- a/mathfuncs.h
+++ b/mathfuncs.h
@@ -11,6 +11,7 @@
 #include "mathfuncs_exp.h"
 #include "mathfuncs_fabs.h"
 #include "mathfuncs_log.h"
+#include "mathfuncs_pow.h"
 #include "mathfuncs_rcp.h"
 #include "mathfuncs_sqrt.h"
 
diff --git a/mathfuncs_base.h b/mathfuncs_base.h
index 59be31a..5fdcafd 100644
--- a/mathfuncs_base.h
+++ b/mathfuncs_base.h
@@ -43,8 +43,11 @@ namespace vecmathlib {
     static realvec_t vml_atanh(realvec_t x);
     
     // convert
+    static realvec_t vml_ceil(realvec_t x);
     static realvec_t vml_convert_float(intvec_t x);
     static intvec_t vml_convert_int(realvec_t x);
+    static realvec_t vml_floor(realvec_t x);
+    static realvec_t vml_round(realvec_t x);
     
     // fabs
     static realvec_t vml_copysign(realvec_t x, realvec_t y);
@@ -65,8 +68,13 @@ namespace vecmathlib {
     static realvec_t vml_log1p(realvec_t x);
     static realvec_t vml_log2(realvec_t x);
     
+    // pow
+    static realvec_t vml_pow(realvec_t x, realvec_t y);
+    
     // rcp
+    static realvec_t vml_fmod(realvec_t x, realvec_t y);
     static realvec_t vml_rcp(realvec_t x);
+    static realvec_t vml_remainder(realvec_t x, realvec_t y);
     
     // sqrt
     static realvec_t vml_rsqrt(realvec_t x);
diff --git a/mathfuncs_convert.h b/mathfuncs_convert.h
index 19f7025..852c94e 100644
--- a/mathfuncs_convert.h
+++ b/mathfuncs_convert.h
@@ -84,6 +84,31 @@ namespace vecmathlib {
     return ix;
   }
   
+  
+  
+  template<typename realvec_t>
+  realvec_t mathfuncs<realvec_t>::vml_round(realvec_t x)
+  {
+    realvec_t r = fabs(x);
+    real_t offset = RV(std::scalbn(R(1.0), FP::mantissa_bits));
+    r += offset;
+#warning "TODO: don't optimise this away!"
+    r -= offset;
+    return copysign(r, x);
+  }
+  
+  template<typename realvec_t>
+  realvec_t mathfuncs<realvec_t>::vml_ceil(realvec_t x)
+  {
+    return round(x + RV(0.5));
+  }
+  
+  template<typename realvec_t>
+  realvec_t mathfuncs<realvec_t>::vml_floor(realvec_t x)
+  {
+    return round(x - RV(0.5));
+  }
+  
 }; // namespace vecmathlib
 
 #endif  // #ifndef MATHFUNCS_CONVERT_H
diff --git a/mathfuncs_pow.h b/mathfuncs_pow.h
new file mode 100644
index 0000000..985117e
--- /dev/null
+++ b/mathfuncs_pow.h
@@ -0,0 +1,26 @@
+// -*-C++-*-
+
+#ifndef MATHFUNCS_POW_H
+#define MATHFUNCS_POW_H
+
+#include "mathfuncs_base.h"
+
+#include <cassert>
+#include <cmath>
+
+
+
+namespace vecmathlib {
+  
+  template<typename realvec_t>
+  realvec_t mathfuncs<realvec_t>::vml_pow(realvec_t x, realvec_t y)
+  {
+    realvec_t r = exp(log(fabs(x)) * y);
+    // The result is negative if x<0 and if y is integer and odd
+    realvec_t sign = fmod(x, RV(2.0)) + RV(1.0);
+    return ifthen(x == RV(0.0), RV(0.0), copysign(r, sign));
+  }
+  
+}; // namespace vecmathlib
+
+#endif  // #ifndef MATHFUNCS_POW_H
diff --git a/mathfuncs_rcp.h b/mathfuncs_rcp.h
index 1704d8f..2b931c9 100644
--- a/mathfuncs_rcp.h
+++ b/mathfuncs_rcp.h
@@ -43,6 +43,22 @@ namespace vecmathlib {
     return r;
   }
   
+  
+  
+  template<typename realvec_t>
+  realvec_t mathfuncs<realvec_t>::vml_remainder(realvec_t x, realvec_t y)
+  {
+    return x - round(x / y) * y;
+  }
+  
+  template<typename realvec_t>
+  realvec_t mathfuncs<realvec_t>::vml_fmod(realvec_t x, realvec_t y)
+  {
+    realvec_t r = fabs(x);
+    y = fabs(y);
+    return copysign(r - floor(r / y) * y, x);
+  }
+  
 }; // namespace vecmathlib
 
 #endif  // #ifndef MATHFUNCS_RCP_H
diff --git a/mathfuncs_sqrt.h b/mathfuncs_sqrt.h
index e33d1cb..c1ea3f7 100644
--- a/mathfuncs_sqrt.h
+++ b/mathfuncs_sqrt.h
@@ -27,7 +27,7 @@ namespace vecmathlib {
     // TODO: divide by M_SQRT2 if ilogb_x % 2 == 1 ?
 #else
     real_t correction =
-      std::scalbn(FP::exponent_offset & 1 ? M_SQRT2 : 1.0,
+      std::scalbn(R(FP::exponent_offset & 1 ? M_SQRT2 : 1.0),
                   FP::exponent_offset >> 1);
     realvec_t r = lsr(x.as_int(), 1).as_float() * RV(correction);
 #endif
diff --git a/test.cc b/test.cc
index dda725e..f75c981 100644
--- a/test.cc
+++ b/test.cc
@@ -185,13 +185,16 @@ struct vecmathlib_test {
   
   static void test_convert()
   {
-    cout << "   testing convert_float convert_int...\n";
+    cout << "   testing ceil convert_float convert_int floor round...\n";
     for (int i=0; i<imax; ++i) {
       realvec_t const x = random(real_t(-10.0), real_t(+10.0));
       intvec_t const n = random(int_t(-10), int_t(+10));
+      check("ceil", ceil, vecmathlib::ceil, x, accuracy);
       check("convert_float",
             FP::convert_float, vecmathlib::convert_float, n, accuracy);
       check("convert_int", FP::convert_int, vecmathlib::convert_int, x);
+      check("floor", floor, vecmathlib::floor, x, accuracy);
+      check("round", round, vecmathlib::round, x, accuracy);
     }
   }
   
@@ -254,31 +257,26 @@ struct vecmathlib_test {
     }
   }
   
-  // static real_t pown(real_t const x, long long const n)
-  // {
-  //   return pow(x, real_t(n));
-  // }
-  // void test_pow()
-  // {
-  //   for (int i=0; i<imax; ++i) {
-  //     real_t const x = random(0.001, 1000.0);
-  //     real_t const y = random(-10.0, +10.0);
-  //     check("pow", pow, vecmathlib::pow, x, y, accuracy);
-  //   }
-  //   for (int i=0; i<imax; ++i) {
-  //     real_t const x = random(-1000.0, +1000.0);
-  //     long long const n = llrandom(-10, +10);
-  //     check("pown", pown, vecmathlib::pown, x, n, accuracy);
-  //   }
-  // }
+  static void test_pow()
+  {
+    cout << "   testing pow...\n";
+    for (int i=0; i<imax; ++i) {
+      realvec_t const x = random(real_t(0.001), real_t(1000.0));
+      realvec_t const y = random(real_t(-10.0), real_t(+10.0));
+      check("pow", pow, vecmathlib::pow, x, y, accuracy);
+    }
+  }
   
   static real_t rcp(real_t x) { return real_t(1.0)/x; }
   static void test_rcp()
   {
-    cout << "   testing rcp...\n";
+    cout << "   testing fmod rcp remainder...\n";
     for (int i=0; i<imax; ++i) {
       realvec_t const x = random(real_t(-10.0), real_t(+10.0));
+      realvec_t const y = random(real_t(-10.0), real_t(+10.0));
+      check("fmod", fmod, vecmathlib::fmod, x, y, accuracy);
       check("rcp", rcp, vecmathlib::rcp, x, accuracy);
+      check("remainder", remainder, vecmathlib::remainder, x, y, accuracy);
     }
   }
   
@@ -329,8 +327,8 @@ struct vecmathlib_test {
     test_asin();
     test_asinh();
     test_exp();
-    // test_log();
-    // test_pow();
+    test_log();
+    test_pow();
     test_rcp();
     // test_sin();
     // test_sinh();
diff --git a/vec_base.h b/vec_base.h
index 91316a9..6fa5f4f 100644
--- a/vec_base.h
+++ b/vec_base.h
@@ -168,7 +168,13 @@ namespace vecmathlib {
   {
     return x.atanh();
   }
-  
+    
+  template<typename real_t, int size>
+  inline realvec<real_t, size> ceil(realvec<real_t, size> x)
+  {
+    return x.ceil();
+  }
+
   template<typename real_t, int size>
   inline realvec<real_t, size> copysign(realvec<real_t, size> x,
                                           realvec<real_t, size> y)
@@ -213,6 +219,13 @@ namespace vecmathlib {
   }
   
   template<typename real_t, int size>
+  inline realvec<real_t, size> fmod(realvec<real_t, size> x,
+                                    realvec<real_t, size> y)
+  {
+    return x.fmod(y);
+  }
+  
+  template<typename real_t, int size>
   inline intvec<real_t, size> ilogb(realvec<real_t, size> x)
   {
     return x.ilogb();
@@ -243,12 +256,32 @@ namespace vecmathlib {
   }
   
   template<typename real_t, int size>
+  inline realvec<real_t, size> pow(realvec<real_t, size> x,
+                                   realvec<real_t, size> y)
+  {
+    return x.pow(y);
+  }
+  
+  template<typename real_t, int size>
   inline realvec<real_t, size> rcp(realvec<real_t, size> x)
   {
     return x.rcp();
   }
   
   template<typename real_t, int size>
+  inline realvec<real_t, size> remainder(realvec<real_t, size> x,
+                                         realvec<real_t, size> y)
+  {
+    return x.remainder(y);
+  }
+  
+  template<typename real_t, int size>
+  inline realvec<real_t, size> round(realvec<real_t, size> x)
+  {
+    return x.round();
+  }
+  
+  template<typename real_t, int size>
   inline realvec<real_t, size> rsqrt(realvec<real_t, size> x)
   {
     return x.rsqrt();
diff --git a/vec_double_avx.h b/vec_double_avx.h
index cca990d..21c6d46 100644
--- a/vec_double_avx.h
+++ b/vec_double_avx.h
@@ -458,6 +458,7 @@ namespace vecmathlib {
     realvec asinh() const { return MF::vml_asinh(*this); }
     realvec atan() const { return MF::vml_atan(*this); }
     realvec atanh() const { return MF::vml_atanh(*this); }
+    realvec ceil() const { return _mm256_ceil_pd(v); }
     realvec copysign(realvec y) const { return MF::vml_copysign(*this, y); }
     realvec exp() const { return MF::vml_exp(*this); }
     realvec exp10() const { return MF::vml_exp10(*this); }
@@ -465,13 +466,17 @@ namespace vecmathlib {
     realvec expm1() const { return MF::vml_expm1(*this); }
     realvec fabs() const { return MF::vml_fabs(*this); }
     realvec floor() const { return _mm256_floor_pd(v); }
+    realvec fmod(realvec y) const { return MF::vml_fmod(*this, y); }
     intvec_t ilogb() const { return MF::vml_ilogb(*this); }
     realvec log() const { return MF::vml_log(*this); }
     realvec log10() const { return MF::vml_log10(*this); }
     realvec log1p() const { return MF::vml_log1p(*this); }
     realvec log2() const { return MF::vml_log2(*this); }
+    realvec pow(realvec y) const { return MF::vml_pow(*this, y); }
     realvec rcp() const { return _mm256_div_pd(_mm256_set1_pd(1.0), v); }
     // realvec rcp() const { return MF::vml_rcp(*this); }
+    realvec remainder(realvec y) const { return MF::vml_remainder(*this, y); }
+    realvec round() const { return _mm256_round_pd(v, _MM_FROUND_NINT); }
     realvec rsqrt() const { return MF::vml_rsqrt(*this); }
     boolvec_t signbit() const { return v; }
     realvec scalbn(intvec_t n) const { return MF::vml_scalbn(*this, n); }
diff --git a/vec_float.h b/vec_float.h
index d89d5e3..8404974 100644
--- a/vec_float.h
+++ b/vec_float.h
@@ -247,22 +247,27 @@ namespace vecmathlib {
     realvec asinh() const { return MF::vml_asinh(*this); }
     realvec atan() const { return MF::vml_atan(*this); }
     realvec atanh() const { return MF::vml_atanh(*this); }
-    realvec copysign(realvec x) const { return MF::vml_copysign(v, x.v); }
+    realvec ceil() const { return MF::vml_ceil(*this); }
+    realvec copysign(realvec y) const { return MF::vml_copysign(*this, y); }
     realvec exp() const { return MF::vml_exp(*this); }
     realvec exp10() const { return MF::vml_exp10(*this); }
     realvec exp2() const { return MF::vml_exp2(*this); }
     realvec expm1() const { return MF::vml_expm1(*this); }
-    realvec fabs() const { return MF::vml_fabs(v); }
-    realvec floor() const { return std::floor(v); }
-    intvec_t ilogb() const { return MF::vml_ilogb(v); }
+    realvec fabs() const { return MF::vml_fabs(*this); }
+    realvec floor() const { return MF::vml_floor(*this); }
+    realvec fmod(realvec y) const { return MF::vml_fmod(*this, y); }
+    intvec_t ilogb() const { return MF::vml_ilogb(*this); }
     realvec log() const { return MF::vml_log(*this); }
     realvec log10() const { return MF::vml_log10(*this); }
     realvec log1p() const { return MF::vml_log1p(*this); }
     realvec log2() const { return MF::vml_log2(*this); }
+    realvec pow(realvec y) const { return MF::vml_pow(*this, y); }
     realvec rcp() const { return MF::vml_rcp(*this); }
+    realvec remainder(realvec y) const { return MF::vml_remainder(*this, y); }
+    realvec round() const { return MF::vml_round(*this); }
     realvec rsqrt() const { return MF::vml_rsqrt(*this); }
-    realvec scalbn(intvec_t n) const { return MF::vml_scalbn(v, n.v); }
-    boolvec_t signbit() const { return MF::vml_signbit(v); }
+    realvec scalbn(intvec_t n) const { return MF::vml_scalbn(*this, n); }
+    boolvec_t signbit() const { return MF::vml_signbit(*this); }
     realvec sqrt() const { return MF::vml_sqrt(*this); }
   };