Correct sqrt and convert_*

3 files changed, 13 insertions, 48 deletions

diff --git a/mathfuncs_convert.h b/mathfuncs_convert.h
index 4dfffef..19f7025 100644
--- a/mathfuncs_convert.h
+++ b/mathfuncs_convert.h
@@ -24,21 +24,24 @@ namespace vecmathlib {
     // Convert lower bits
     intvec_t xlo = x & IV((U(1) << lobits) - 1);
     // exponent for the equivalent floating point number
-    int_t exponent_lo = FP::exponent_offset + lobits;
+    int_t exponent_lo = (FP::exponent_offset + lobits) << FP::mantissa_bits;
     xlo |= exponent_lo;
     // subtract hidden mantissa bit
     realvec_t flo = as_float(xlo) - RV(FP::as_float(exponent_lo));
     
     // Convert upper bits
     // make unsigned by subtracting largest negative number
+    // (only do this for the high bits, since they have sufficient
+    // precision to handle the overflow)
     x ^= FP::sign_mask;
     intvec_t xhi = lsr(x, lobits);
     // exponent for the equivalent floating point number
-    int_t exponent_hi = FP::exponent_offset + FP::bits;
+    int_t exponent_hi = (FP::exponent_offset + 2*lobits) << FP::mantissa_bits;
     xhi |= exponent_hi;
-    // add largest negative number, and subtract hidden mantissa bit
-    realvec_t
-      fhi = as_float(xhi) - RV(FP::as_float(exponent_hi) + R(FP::sign_mask));
+    // subtract hidden mantissa bit
+    realvec_t fhi = as_float(xhi) - RV(FP::as_float(exponent_hi));
+    // add largest negative number again
+    fhi -= RV(R(FP::sign_mask));
     
     // Combine results
     return flo + fhi;
@@ -76,7 +79,7 @@ namespace vecmathlib {
     // Handle negative numbers
     ix = ifthen(is_negative, -ix, ix);
     // Handle zero
-    ix = ifthen(is_zero, IV(0), ix);
+    ix = ifthen(is_zero, IV(I(0)), ix);
     
     return ix;
   }
diff --git a/mathfuncs_sqrt.h b/mathfuncs_sqrt.h
index 2c43478..e33d1cb 100644
--- a/mathfuncs_sqrt.h
+++ b/mathfuncs_sqrt.h
@@ -15,7 +15,6 @@ namespace vecmathlib {
   template<typename realvec_t>
   realvec_t mathfuncs<realvec_t>::vml_sqrt(realvec_t x)
   {
-    std::cout << "mathfuncs::sqrt\n";
     // Handle special case: zero
     boolvec_t is_zero = x <= RV(0.0);
     x = ifthen(is_zero, RV(1.0), x);
@@ -28,18 +27,15 @@ namespace vecmathlib {
     // TODO: divide by M_SQRT2 if ilogb_x % 2 == 1 ?
 #else
     real_t correction =
-      std::scalbn(4.0/3.0, FP::exponent_offset >> 1) *
-      (FP::exponent_offset & 1 ? 1.0 / M_SQRT1_2 : 1.0);
+      std::scalbn(FP::exponent_offset & 1 ? M_SQRT2 : 1.0,
+                  FP::exponent_offset >> 1);
     realvec_t r = lsr(x.as_int(), 1).as_float() * RV(correction);
-    real_t one=1.0;
-    std::cout << "x=" << x << " r=" << r << " eo=" << FP::exponent_offset << " corr=" << correction << " ione=" << RV(one).as_int() << "\n";
 #endif
     
     // Iterate
     // nmax iterations give an accuracy of 2^nmax binary digits. 5
     // iterations suffice for double precision with its 53 digits.
     int const nmax = 5;
-    std::cout << "bits" << FP::bits << "nmax" << nmax << "\n";
     for (int n=1; n<nmax; ++n) {
       // Step
       assert(all(x > RV(0.0)));
diff --git a/vec_double_avx.h b/vec_double_avx.h
index 258ffd0..cbea5c2 100644
--- a/vec_double_avx.h
+++ b/vec_double_avx.h
@@ -395,27 +395,7 @@ namespace vecmathlib {
     
     
     intvec_t as_int() const { return _mm256_castpd_si256(v); }
-    intvec_t convert_int() const
-    {
-#if 0
-      __m128i iv0123 = _mm256_cvtpd_epi32(v);
-      __m128i iv2301 = _mm_shuffle_ps(iv0123, iv0123, 0b10110001);
-      __m256i iv01232301 =
-        _mm256_insertf128_si256(_mm256_castsi128_si256(iv0123), iv2301, 1);
-      __m256i zero = _mm256_setzero_ps();
-      return _mm256_unpacklo_ps(iv01232301, zero);
-#else
-      realvec x = _mm256_floor_pd(v);
-      intvec_t ix = x.as_int();
-      boolvec_t sign = ix.as_bool();
-      intvec_t exponent = (ix & exponent_mask) >> mantissa_bits;
-      ix &= mantissa_mask;
-      ix |= U(1) << mantissa_bits;  // add hidden bit
-      ix <<= exponent - exponent_offset + 52; // ???
-      ix = ifthen(sign, -ix, ix);
-      return ix;
-#endif
-    }
+    intvec_t convert_int() const { return MF::vml_convert_int(*this); }
     
     
     
@@ -566,21 +546,7 @@ namespace vecmathlib {
   
   inline auto intvec<double,4>::convert_float() const -> realvec_t
   {
-    intvec x = v;
-    uint_t signbitmask = U(1) << (bits-1);
-    // make unsigned
-    x += signbitmask;
-    // convert lower 52 bits
-    intvec xlo = x & IV((U(1) << mantissa_bits) - 1);
-    int_t exponent_0 = exponent_offset + mantissa_bits;
-    xlo = xlo | exponent_0;
-    realvec_t flo = xlo.as_float() - FP::as_float(exponent_0);
-    // convert upper 22 bits
-    intvec xhi = x.lsr(U(mantissa_bits));
-    int_t exponent_52 = exponent_0 + mantissa_bits;
-    xhi = xhi | exponent_52;
-    realvec_t fhi = xhi.as_float() - FP::as_float(exponent_52);
-    return flo + fhi - R(signbitmask);
+    return MF::vml_convert_float(*this);
   }
   
 } // namespace vecmathlib