try to fix #170

Signed-off-by: Benjamin Schultzer <[email protected]>

Author: Schultzer

src/math/jnf.rs

@@ -75,14 +75,14 @@ pub fn jnf(n: i32, mut x: f32) -> f32 {
             }
             temp = 0.5 * x;
             b = temp;
-            a = 1.0;
+            let mut a = 1;
             i = 2;
             while i <= nm1 + 1 {
-                a *= i as f32; /* a = n! */
+                a *= i; /* a = n! */
                 b *= temp; /* b = (x/2)^n */
                 i += 1;
             }
-            b = b / a;
+            b /= a as f32;
         } else {
             /* use backward recurrence */
             /*                      x      x^2      x^2
@@ -124,7 +124,7 @@ pub fn jnf(n: i32, mut x: f32) -> f32 {
             let mut k: i32;
 
             nf = (nm1 as f32) + 1.0;
-            w = 2.0 * (nf as f32) / x;
+            w = 2.0 * nf / x;
             h = 2.0 / x;
             z = w + h;
             q0 = w;
@@ -169,10 +169,10 @@ pub fn jnf(n: i32, mut x: f32) -> f32 {
                     b = 2.0 * (i as f32) * b / x - a;
                     a = temp;
                     /* scale b to avoid spurious overflow */
-                    let x1p60 = f32::from_bits(0x5d800000); // 0x1p60 == 2^60
-                    if b > x1p60 {
+                    // let x1p60 = f32::from_bits(0x5d800000); // 0x1p60 == 2^60
+                    if b > 1.152922e+18 {
                         a /= b;
-                        t /= b;
+                        t = t / b;
                         b = 1.0;
                     }
                     i -= 1;
@@ -182,8 +182,9 @@ pub fn jnf(n: i32, mut x: f32) -> f32 {
             w = j1f(x);
             if fabsf(z) >= fabsf(w) {
                 b = t * z / b;
+                // panic!("{:?} {:?} {:?} {:?} {:?} {:?}", t, z, w, b, f32::from_bits(0x5d800000), 1.152922e+18);
             } else {
-                b = t * w / a;
+                b = t * (w / a);
             }
         }
     }

src/math/logf.rs

@@ -1,66 +1,69 @@
-/* origin: FreeBSD /usr/src/lib/msun/src/e_logf.c */
-/*
- * Conversion to float by Ian Lance Taylor, Cygnus Support, [email protected].
- */
-/*
- * ====================================================
- * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
- *
- * Developed at SunPro, a Sun Microsystems, Inc. business.
- * Permission to use, copy, modify, and distribute this
- * software is freely granted, provided that this notice
- * is preserved.
- * ====================================================
- */
+const LOGF_TABLE_BITS: u32 = 4;
+const T: [(f64, f64); 16] = [(1.398907162146528e0, -3.3569133332882284e-1),
+                             (1.3403141896637998e0, -2.929040563774074e-1),
+                             (1.286432210124115e0, -2.518726580937369e-1),
+                             (1.2367150214269895e0, -2.1245868807117255e-1),
+                             (1.1906977166711752e0, -1.7453945183745634e-1),
+                             (1.1479821020556429e0, -1.380057072319758e-1),
+                             (1.1082251448272158e0, -1.0275976698545139e-1),
+                             (1.0711297413057381e0, -6.871392447020525e-2),
+                             (1.036437278977283e0, -3.57891387398228e-2),
+                             (1e0, 0e0),
+                             (9.492859795739057e-1, 5.204517742929496e-2),
+                             (8.951049428609004e-1, 1.1081431298787942e-1),
+                             (8.476821620351103e-1, 1.652495223695143e-1),
+                             (8.050314851692001e-1, 2.1687389031699977e-1),
+                             (7.664671008843108e-1, 2.659635028121397e-1),
+                             (7.31428603316328e-1, 3.127556664073557e-1)];
 
-const LN2_HI: f32 = 6.9313812256e-01; /* 0x3f317180 */
-const LN2_LO: f32 = 9.0580006145e-06; /* 0x3717f7d1 */
-/* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */
-const LG1: f32 = 0.66666662693; /*  0xaaaaaa.0p-24*/
-const LG2: f32 = 0.40000972152; /*  0xccce13.0p-25 */
-const LG3: f32 = 0.28498786688; /*  0x91e9ee.0p-25 */
-const LG4: f32 = 0.24279078841; /*  0xf89e26.0p-26 */
+const A: [f64; 3] = [-2.5089342214237154e-1, 3.33456765744066e-1, -4.999997485802103e-1];
+const LN2: f64 = 6.931471805599453e-1;
+const N: u32 =  (1 << LOGF_TABLE_BITS);
+const OFF: u32 = 0x3f330000;
 
 #[inline]
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
-pub fn logf(mut x: f32) -> f32 {
-    let x1p25 = f32::from_bits(0x4c000000); // 0x1p25f === 2 ^ 25
-
+pub fn logf(x: f32) -> f32 {
     let mut ix = x.to_bits();
-    let mut k = 0i32;
 
-    if (ix < 0x00800000) || ((ix >> 31) != 0) {
-        /* x < 2**-126  */
-        if ix << 1 == 0 {
-            return -1. / (x * x); /* log(+-0)=-inf */
-        }
-        if (ix >> 31) != 0 {
-            return (x - x) / 0.; /* log(-#) = NaN */
-        }
-        /* subnormal number, scale up x */
-        k -= 25;
-        x *= x1p25;
-        ix = x.to_bits();
-    } else if ix >= 0x7f800000 {
-        return x;
-    } else if ix == 0x3f800000 {
+    /* Fix sign of zero with downward rounding when x==1.  */
+    if ix == 0x3f800000 {
         return 0.;
     }
+    if ix - 0x00800000 >= 0x7f800000 - 0x00800000 {
+    /* x < 0x1p-126 or inf or nan.  */
+    if ix * 2 == 0 {
+        return  -1. / 0.;
+    }
+    if ix == 0x7f800000 { /* log(inf) == inf.  */
+        return x;
+    }
+    if (ix & 0x80000000) != 0 || ix * 2 >= 0xff000000 {
+        return (x - x) / (x - x);
+    }
+    /* x is subnormal, normalize it.  */
+    ix = f32::to_bits(x * 8.388608e6);
+        ix -= 23 << 23;
+    }
+
+    /* x = 2^k z; where z is in range [OFF,2*OFF] and exact.
+    The range is split into N subintervals.
+    The ith subinterval contains z and c is near its center.  */
+	let tmp = ix - OFF;
+	let i = (tmp >> (23 - LOGF_TABLE_BITS)) % N;
+	let k = tmp as i32 >> 23; /* arithmetic shift */
+	let iz = ix - (tmp & 0x1ff << 23);
+	let (invc, logc) = T[i as usize];
+	let z = f32::from_bits(iz) as f64;
 
-    /* reduce x into [sqrt(2)/2, sqrt(2)] */
-    ix += 0x3f800000 - 0x3f3504f3;
-    k += ((ix >> 23) as i32) - 0x7f;
-    ix = (ix & 0x007fffff) + 0x3f3504f3;
-    x = f32::from_bits(ix);
+	/* log(x) = log1p(z/c-1) + log(c) + k*Ln2 */
+	let r = z * invc - 1.;
+	let y0 = logc + k as f64 * LN2;
 
-    let f = x - 1.;
-    let s = f / (2. + f);
-    let z = s * s;
-    let w = z * z;
-    let t1 = w * (LG2 + w * LG4);
-    let t2 = z * (LG1 + w * LG3);
-    let r = t2 + t1;
-    let hfsq = 0.5 * f * f;
-    let dk = k as f32;
-    s * (hfsq + r) + dk * LN2_LO - hfsq + f + dk * LN2_HI
+	/* Pipelined polynomial evaluation to approximate log1p(r).  */
+	let r2 = r * r;
+	let mut y = A[1] * r + A[2];
+	y = A[0] * r2 + y;
+	y = y * r2 + (y0 + r);
+	return y as f32;
 }