libm - cargo Package Compare versions

Comparing version 0.1.0 to 0.1.1

src/math/acos.rs

+/* origin: FreeBSD /usr/src/lib/msun/src/e_acos.c */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+/* acos(x)
+ * Method :
+ *      acos(x)  = pi/2 - asin(x)
+ *      acos(-x) = pi/2 + asin(x)
+ * For |x|<=0.5
+ *      acos(x) = pi/2 - (x + x*x^2*R(x^2))     (see asin.c)
+ * For x>0.5
+ *      acos(x) = pi/2 - (pi/2 - 2asin(sqrt((1-x)/2)))
+ *              = 2asin(sqrt((1-x)/2))
+ *              = 2s + 2s*z*R(z)        ...z=(1-x)/2, s=sqrt(z)
+ *              = 2f + (2c + 2s*z*R(z))
+ *     where f=hi part of s, and c = (z-f*f)/(s+f) is the correction term
+ *     for f so that f+c ~ sqrt(z).
+ * For x<-0.5
+ *      acos(x) = pi - 2asin(sqrt((1-|x|)/2))
+ *              = pi - 0.5*(s+s*z*R(z)), where z=(1-|x|)/2,s=sqrt(z)
+ *
+ * Special cases:
+ *      if x is NaN, return x itself;
+ *      if |x|>1, return NaN with invalid signal.
+ *
+ * Function needed: sqrt
+ */
+
+use super::sqrt;
+
+const PIO2_HI: f64 = 1.57079632679489655800e+00; /* 0x3FF921FB, 0x54442D18 */
+const PIO2_LO: f64 = 6.12323399573676603587e-17; /* 0x3C91A626, 0x33145C07 */
+const PS0: f64 = 1.66666666666666657415e-01; /* 0x3FC55555, 0x55555555 */
+const PS1: f64 = -3.25565818622400915405e-01; /* 0xBFD4D612, 0x03EB6F7D */
+const PS2: f64 = 2.01212532134862925881e-01; /* 0x3FC9C155, 0x0E884455 */
+const PS3: f64 = -4.00555345006794114027e-02; /* 0xBFA48228, 0xB5688F3B */
+const PS4: f64 = 7.91534994289814532176e-04; /* 0x3F49EFE0, 0x7501B288 */
+const PS5: f64 = 3.47933107596021167570e-05; /* 0x3F023DE1, 0x0DFDF709 */
+const QS1: f64 = -2.40339491173441421878e+00; /* 0xC0033A27, 0x1C8A2D4B */
+const QS2: f64 = 2.02094576023350569471e+00; /* 0x40002AE5, 0x9C598AC8 */
+const QS3: f64 = -6.88283971605453293030e-01; /* 0xBFE6066C, 0x1B8D0159 */
+const QS4: f64 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
+
+#[inline]
+fn r(z: f64) -> f64 {
+    let p: f64 = z * (PS0 + z * (PS1 + z * (PS2 + z * (PS3 + z * (PS4 + z * PS5)))));
+    let q: f64 = 1.0 + z * (QS1 + z * (QS2 + z * (QS3 + z * QS4)));
+    return p / q;
+}
+
+#[inline]
+pub fn acos(x: f64) -> f64 {
+    let x1p_120f = f64::from_bits(0x3870000000000000); // 0x1p-120 === 2 ^ -120
+    let z: f64;
+    let w: f64;
+    let s: f64;
+    let c: f64;
+    let df: f64;
+    let hx: u32;
+    let ix: u32;
+
+    hx = (x.to_bits() >> 32) as u32;
+    ix = hx & 0x7fffffff;
+    /* |x| >= 1 or nan */
+    if ix >= 0x3ff00000 {
+        let lx: u32 = x.to_bits() as u32;
+
+        if (ix - 0x3ff00000 | lx) == 0 {
+            /* acos(1)=0, acos(-1)=pi */
+            if (hx >> 31) != 0 {
+                return 2. * PIO2_HI + x1p_120f;
+            }
+            return 0.;
+        }
+        return 0. / (x - x);
+    }
+    /* |x| < 0.5 */
+    if ix < 0x3fe00000 {
+        if ix <= 0x3c600000 {
+            /* |x| < 2**-57 */
+            return PIO2_HI + x1p_120f;
+        }
+        return PIO2_HI - (x - (PIO2_LO - x * r(x * x)));
+    }
+    /* x < -0.5 */
+    if (hx >> 31) != 0 {
+        z = (1.0 + x) * 0.5;
+        s = sqrt(z);
+        w = r(z) * s - PIO2_LO;
+        return 2. * (PIO2_HI - (s + w));
+    }
+    /* x > 0.5 */
+    z = (1.0 - x) * 0.5;
+    s = sqrt(z);
+    // Set the low 4 bytes to zero
+    df = f64::from_bits(s.to_bits() & 0xff_ff_ff_ff_00_00_00_00);
+
+    c = (z - df * df) / (s + df);
+    w = r(z) * s + c;
+    return 2. * (df + w);
+}

src/math/acosf.rs

+use super::sqrtf::sqrtf;
+
+const PIO2_HI: f32 = 1.5707962513e+00; /* 0x3fc90fda */
+const PIO2_LO: f32 = 7.5497894159e-08; /* 0x33a22168 */
+const P_S0: f32 = 1.6666586697e-01;
+const P_S1: f32 = -4.2743422091e-02;
+const P_S2: f32 = -8.6563630030e-03;
+const Q_S1: f32 = -7.0662963390e-01;
+
+fn r(z: f32) -> f32 {
+    let p = z * (P_S0 + z * (P_S1 + z * P_S2));
+    let q = 1. + z * Q_S1;
+    p / q
+}
+
+#[inline]
+pub fn acosf(x: f32) -> f32 {
+    let x1p_120 = f32::from_bits(0x03800000); // 0x1p-120 === 2 ^ (-120)
+
+    let z: f32;
+    let w: f32;
+    let s: f32;
+
+    let mut hx = x.to_bits();
+    let ix = hx & 0x7fffffff;
+    /* |x| >= 1 or nan */
+    if ix >= 0x3f800000 {
+        if ix == 0x3f800000 {
+            if (hx >> 31) != 0 {
+                return 2. * PIO2_HI + x1p_120;
+            }
+            return 0.;
+        }
+        return 0. / (x - x);
+    }
+    /* |x| < 0.5 */
+    if ix < 0x3f000000 {
+        if ix <= 0x32800000 {
+            /* |x| < 2**-26 */
+            return PIO2_HI + x1p_120;
+        }
+        return PIO2_HI - (x - (PIO2_LO - x * r(x * x)));
+    }
+    /* x < -0.5 */
+    if (hx >> 31) != 0 {
+        z = (1. + x) * 0.5;
+        s = sqrtf(z);
+        w = r(z) * s - PIO2_LO;
+        return 2. * (PIO2_HI - (s + w));
+    }
+    /* x > 0.5 */
+    z = (1. - x) * 0.5;
+    s = sqrtf(z);
+    hx = s.to_bits();
+    let df = f32::from_bits(hx & 0xfffff000);
+    let c = (z - df * df) / (s + df);
+    w = r(z) * s + c;
+    2. * (df + w)
+}

src/math/asin.rs

+/* origin: FreeBSD /usr/src/lib/msun/src/e_asin.c */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+/* asin(x)
+ * Method :
+ *      Since  asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
+ *      we approximate asin(x) on [0,0.5] by
+ *              asin(x) = x + x*x^2*R(x^2)
+ *      where
+ *              R(x^2) is a rational approximation of (asin(x)-x)/x^3
+ *      and its remez error is bounded by
+ *              |(asin(x)-x)/x^3 - R(x^2)| < 2^(-58.75)
+ *
+ *      For x in [0.5,1]
+ *              asin(x) = pi/2-2*asin(sqrt((1-x)/2))
+ *      Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2;
+ *      then for x>0.98
+ *              asin(x) = pi/2 - 2*(s+s*z*R(z))
+ *                      = pio2_hi - (2*(s+s*z*R(z)) - pio2_lo)
+ *      For x<=0.98, let pio4_hi = pio2_hi/2, then
+ *              f = hi part of s;
+ *              c = sqrt(z) - f = (z-f*f)/(s+f)         ...f+c=sqrt(z)
+ *      and
+ *              asin(x) = pi/2 - 2*(s+s*z*R(z))
+ *                      = pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo)
+ *                      = pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c))
+ *
+ * Special cases:
+ *      if x is NaN, return x itself;
+ *      if |x|>1, return NaN with invalid signal.
+ *
+ */
+
+use super::{fabs, get_high_word, get_low_word, sqrt, with_set_low_word};
+
+const PIO2_HI: f64 = 1.57079632679489655800e+00; /* 0x3FF921FB, 0x54442D18 */
+const PIO2_LO: f64 = 6.12323399573676603587e-17; /* 0x3C91A626, 0x33145C07 */
+/* coefficients for R(x^2) */
+const P_S0: f64 = 1.66666666666666657415e-01; /* 0x3FC55555, 0x55555555 */
+const P_S1: f64 = -3.25565818622400915405e-01; /* 0xBFD4D612, 0x03EB6F7D */
+const P_S2: f64 = 2.01212532134862925881e-01; /* 0x3FC9C155, 0x0E884455 */
+const P_S3: f64 = -4.00555345006794114027e-02; /* 0xBFA48228, 0xB5688F3B */
+const P_S4: f64 = 7.91534994289814532176e-04; /* 0x3F49EFE0, 0x7501B288 */
+const P_S5: f64 = 3.47933107596021167570e-05; /* 0x3F023DE1, 0x0DFDF709 */
+const Q_S1: f64 = -2.40339491173441421878e+00; /* 0xC0033A27, 0x1C8A2D4B */
+const Q_S2: f64 = 2.02094576023350569471e+00; /* 0x40002AE5, 0x9C598AC8 */
+const Q_S3: f64 = -6.88283971605453293030e-01; /* 0xBFE6066C, 0x1B8D0159 */
+const Q_S4: f64 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
+
+fn comp_r(z: f64) -> f64 {
+    let p = z * (P_S0 + z * (P_S1 + z * (P_S2 + z * (P_S3 + z * (P_S4 + z * P_S5)))));
+    let q = 1.0 + z * (Q_S1 + z * (Q_S2 + z * (Q_S3 + z * Q_S4)));
+    return p / q;
+}
+
+pub fn asin(mut x: f64) -> f64 {
+    let z: f64;
+    let r: f64;
+    let s: f64;
+    let hx: u32;
+    let ix: u32;
+
+    hx = get_high_word(x);
+    ix = hx & 0x7fffffff;
+    /* |x| >= 1 or nan */
+    if ix >= 0x3ff00000 {
+        let lx: u32;
+        lx = get_low_word(x);
+        if (ix - 0x3ff00000 | lx) == 0 {
+            /* asin(1) = +-pi/2 with inexact */
+            return x * PIO2_HI + f64::from_bits(0x3870000000000000);
+        } else {
+            return 0.0 / (x - x);
+        }
+    }
+    /* |x| < 0.5 */
+    if ix < 0x3fe00000 {
+        /* if 0x1p-1022 <= |x| < 0x1p-26, avoid raising underflow */
+        if ix < 0x3e500000 && ix >= 0x00100000 {
+            return x;
+        } else {
+            return x + x * comp_r(x * x);
+        }
+    }
+    /* 1 > |x| >= 0.5 */
+    z = (1.0 - fabs(x)) * 0.5;
+    s = sqrt(z);
+    r = comp_r(z);
+    if ix >= 0x3fef3333 {
+        /* if |x| > 0.975 */
+        x = PIO2_HI - (2. * (s + s * r) - PIO2_LO);
+    } else {
+        let f: f64;
+        let c: f64;
+        /* f+c = sqrt(z) */
+        f = with_set_low_word(s, 0);
+        c = (z - f * f) / (s + f);
+        x = 0.5 * PIO2_HI - (2.0 * s * r - (PIO2_LO - 2.0 * c) - (0.5 * PIO2_HI - 2.0 * f));
+    }
+    if hx >> 31 != 0 {
+        return -x;
+    } else {
+        return x;
+    }
+}

src/math/asinf.rs

+use super::fabsf::fabsf;
+use super::sqrt::sqrt;
+
+const PIO2: f64 = 1.570796326794896558e+00;
+
+/* coefficients for R(x^2) */
+const P_S0: f32 = 1.6666586697e-01;
+const P_S1: f32 = -4.2743422091e-02;
+const P_S2: f32 = -8.6563630030e-03;
+const Q_S1: f32 = -7.0662963390e-01;
+
+fn r(z: f32) -> f32 {
+    let p = z * (P_S0 + z * (P_S1 + z * P_S2));
+    let q = 1. + z * Q_S1;
+    p / q
+}
+
+#[inline]
+pub fn asinf(mut x: f32) -> f32 {
+    let x1p_120 = f64::from_bits(0x3870000000000000); // 0x1p-120 === 2 ^ (-120)
+
+    let hx = x.to_bits();
+    let ix = hx & 0x7fffffff;
+
+    if ix >= 0x3f800000 {
+        /* |x| >= 1 */
+        if ix == 0x3f800000 {
+            /* |x| == 1 */
+            return ((x as f64) * PIO2 + x1p_120) as f32; /* asin(+-1) = +-pi/2 with inexact */
+        }
+        return 0. / (x - x); /* asin(|x|>1) is NaN */
+    }
+
+    if ix < 0x3f000000 {
+        /* |x| < 0.5 */
+        /* if 0x1p-126 <= |x| < 0x1p-12, avoid raising underflow */
+        if (ix < 0x39800000) && (ix >= 0x00800000) {
+            return x;
+        }
+        return x + x * r(x * x);
+    }
+
+    /* 1 > |x| >= 0.5 */
+    let z = (1. - fabsf(x)) * 0.5;
+    let s = sqrt(z as f64);
+    x = (PIO2 - 2. * (s + s * (r(z) as f64))) as f32;
+    if (hx >> 31) != 0 {
+        -x
+    } else {
+        x
+    }
+}

src/math/atanf.rs

+use super::fabsf;
+
+const ATAN_HI: [f32; 4] = [
+    4.6364760399e-01, /* atan(0.5)hi 0x3eed6338 */
+    7.8539812565e-01, /* atan(1.0)hi 0x3f490fda */
+    9.8279368877e-01, /* atan(1.5)hi 0x3f7b985e */
+    1.5707962513e+00, /* atan(inf)hi 0x3fc90fda */
+];
+
+const ATAN_LO: [f32; 4] = [
+    5.0121582440e-09, /* atan(0.5)lo 0x31ac3769 */
+    3.7748947079e-08, /* atan(1.0)lo 0x33222168 */
+    3.4473217170e-08, /* atan(1.5)lo 0x33140fb4 */
+    7.5497894159e-08, /* atan(inf)lo 0x33a22168 */
+];
+
+const A_T: [f32; 5] = [
+    3.3333328366e-01,
+    -1.9999158382e-01,
+    1.4253635705e-01,
+    -1.0648017377e-01,
+    6.1687607318e-02,
+];
+
+#[inline]
+pub fn atanf(mut x: f32) -> f32 {
+    let x1p_120 = f32::from_bits(0x03800000); // 0x1p-120 === 2 ^ (-120)
+
+    let z: f32;
+
+    let mut ix = x.to_bits();
+    let sign = (ix >> 31) != 0;
+    ix &= 0x7fffffff;
+
+    if ix >= 0x4c800000 {
+        /* if |x| >= 2**26 */
+        if x.is_nan() {
+            return x;
+        }
+        z = ATAN_HI[3] + x1p_120;
+        return if sign { -z } else { z };
+    }
+    let id = if ix < 0x3ee00000 {
+        /* |x| < 0.4375 */
+        if ix < 0x39800000 {
+            /* |x| < 2**-12 */
+            if ix < 0x00800000 {
+                /* raise underflow for subnormal x */
+                force_eval!(x * x);
+            }
+            return x;
+        }
+        -1
+    } else {
+        x = fabsf(x);
+        if ix < 0x3f980000 {
+            /* |x| < 1.1875 */
+            if ix < 0x3f300000 {
+                /*  7/16 <= |x| < 11/16 */
+                x = (2. * x - 1.) / (2. + x);
+                0
+            } else {
+                /* 11/16 <= |x| < 19/16 */
+                x = (x - 1.) / (x + 1.);
+                1
+            }
+        } else {
+            if ix < 0x401c0000 {
+                /* |x| < 2.4375 */
+                x = (x - 1.5) / (1. + 1.5 * x);
+                2
+            } else {
+                /* 2.4375 <= |x| < 2**26 */
+                x = -1. / x;
+                3
+            }
+        }
+    };
+    /* end of argument reduction */
+    z = x * x;
+    let w = z * z;
+    /* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
+    let s1 = z * (A_T[0] + w * (A_T[2] + w * A_T[4]));
+    let s2 = w * (A_T[1] + w * A_T[3]);
+    if id < 0 {
+        return x - x * (s1 + s2);
+    }
+    let id = id as usize;
+    let z = ATAN_HI[id] - ((x * (s1 + s2) - ATAN_LO[id]) - x);
+    if sign {
+        -z
+    } else {
+        z
+    }
+}

src/math/cbrt.rs

+/* origin: FreeBSD /usr/src/lib/msun/src/s_cbrt.c */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ *
+ * Optimized by Bruce D. Evans.
+ */
+/* cbrt(x)
+ * Return cube root of x
+ */
+
+use core::f64;
+
+const B1: u32 = 715094163; /* B1 = (1023-1023/3-0.03306235651)*2**20 */
+const B2: u32 = 696219795; /* B2 = (1023-1023/3-54/3-0.03306235651)*2**20 */
+
+/* |1/cbrt(x) - p(x)| < 2**-23.5 (~[-7.93e-8, 7.929e-8]). */
+const P0: f64 = 1.87595182427177009643; /* 0x3ffe03e6, 0x0f61e692 */
+const P1: f64 = -1.88497979543377169875; /* 0xbffe28e0, 0x92f02420 */
+const P2: f64 = 1.621429720105354466140; /* 0x3ff9f160, 0x4a49d6c2 */
+const P3: f64 = -0.758397934778766047437; /* 0xbfe844cb, 0xbee751d9 */
+const P4: f64 = 0.145996192886612446982; /* 0x3fc2b000, 0xd4e4edd7 */
+
+#[inline]
+pub fn cbrt(x: f64) -> f64 {
+    let x1p54 = f64::from_bits(0x4350000000000000); // 0x1p54 === 2 ^ 54
+
+    let mut ui: u64 = x.to_bits();
+    let mut r: f64;
+    let s: f64;
+    let mut t: f64;
+    let w: f64;
+    let mut hx: u32 = (ui >> 32) as u32 & 0x7fffffff;
+
+    if hx >= 0x7ff00000 {
+        /* cbrt(NaN,INF) is itself */
+        return x + x;
+    }
+
+    /*
+     * Rough cbrt to 5 bits:
+     *    cbrt(2**e*(1+m) ~= 2**(e/3)*(1+(e%3+m)/3)
+     * where e is integral and >= 0, m is real and in [0, 1), and "/" and
+     * "%" are integer division and modulus with rounding towards minus
+     * infinity.  The RHS is always >= the LHS and has a maximum relative
+     * error of about 1 in 16.  Adding a bias of -0.03306235651 to the
+     * (e%3+m)/3 term reduces the error to about 1 in 32. With the IEEE
+     * floating point representation, for finite positive normal values,
+     * ordinary integer divison of the value in bits magically gives
+     * almost exactly the RHS of the above provided we first subtract the
+     * exponent bias (1023 for doubles) and later add it back.  We do the
+     * subtraction virtually to keep e >= 0 so that ordinary integer
+     * division rounds towards minus infinity; this is also efficient.
+     */
+    if hx < 0x00100000 {
+        /* zero or subnormal? */
+        ui = (x * x1p54).to_bits();
+        hx = (ui >> 32) as u32 & 0x7fffffff;
+        if hx == 0 {
+            return x; /* cbrt(0) is itself */
+        }
+        hx = hx / 3 + B2;
+    } else {
+        hx = hx / 3 + B1;
+    }
+    ui &= 1 << 63;
+    ui |= (hx as u64) << 32;
+    t = f64::from_bits(ui);
+
+    /*
+     * New cbrt to 23 bits:
+     *    cbrt(x) = t*cbrt(x/t**3) ~= t*P(t**3/x)
+     * where P(r) is a polynomial of degree 4 that approximates 1/cbrt(r)
+     * to within 2**-23.5 when |r - 1| < 1/10.  The rough approximation
+     * has produced t such than |t/cbrt(x) - 1| ~< 1/32, and cubing this
+     * gives us bounds for r = t**3/x.
+     *
+     * Try to optimize for parallel evaluation as in __tanf.c.
+     */
+    r = (t * t) * (t / x);
+    t = t * ((P0 + r * (P1 + r * P2)) + ((r * r) * r) * (P3 + r * P4));
+
+    /*
+     * Round t away from zero to 23 bits (sloppily except for ensuring that
+     * the result is larger in magnitude than cbrt(x) but not much more than
+     * 2 23-bit ulps larger).  With rounding towards zero, the error bound
+     * would be ~5/6 instead of ~4/6.  With a maximum error of 2 23-bit ulps
+     * in the rounded t, the infinite-precision error in the Newton
+     * approximation barely affects third digit in the final error
+     * 0.667; the error in the rounded t can be up to about 3 23-bit ulps
+     * before the final error is larger than 0.667 ulps.
+     */
+    ui = t.to_bits();
+    ui = (ui + 0x80000000) & 0xffffffffc0000000;
+    t = f64::from_bits(ui);
+
+    /* one step Newton iteration to 53 bits with error < 0.667 ulps */
+    s = t * t; /* t*t is exact */
+    r = x / s; /* error <= 0.5 ulps; |r| < |t| */
+    w = t + t; /* t+t is exact */
+    r = (r - t) / (w + r); /* r-t is exact; w+r ~= 3*t */
+    t = t + t * r; /* error <= 0.5 + 0.5/3 + epsilon */
+    t
+}

src/math/cbrtf.rs

+/* origin: FreeBSD /usr/src/lib/msun/src/s_cbrtf.c */
+/*
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, [e-mail].
+ * Debugged and optimized by Bruce D. Evans.
+ */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+/* cbrtf(x)
+ * Return cube root of x
+ */
+
+use core::f32;
+
+const B1: u32 = 709958130; /* B1 = (127-127.0/3-0.03306235651)*2**23 */
+const B2: u32 = 642849266; /* B2 = (127-127.0/3-24/3-0.03306235651)*2**23 */
+
+#[inline]
+pub fn cbrtf(x: f32) -> f32 {
+    let x1p24 = f32::from_bits(0x4b800000); // 0x1p24f === 2 ^ 24
+
+    let mut r: f64;
+    let mut t: f64;
+    let mut ui: u32 = x.to_bits();
+    let mut hx: u32 = ui & 0x7fffffff;
+
+    if hx >= 0x7f800000 {
+        /* cbrt(NaN,INF) is itself */
+        return x + x;
+    }
+
+    /* rough cbrt to 5 bits */
+    if hx < 0x00800000 {
+        /* zero or subnormal? */
+        if hx == 0 {
+            return x; /* cbrt(+-0) is itself */
+        }
+        ui = (x * x1p24).to_bits();
+        hx = ui & 0x7fffffff;
+        hx = hx / 3 + B2;
+    } else {
+        hx = hx / 3 + B1;
+    }
+    ui &= 0x80000000;
+    ui |= hx;
+
+    /*
+     * First step Newton iteration (solving t*t-x/t == 0) to 16 bits.  In
+     * double precision so that its terms can be arranged for efficiency
+     * without causing overflow or underflow.
+     */
+    t = f32::from_bits(ui) as f64;
+    r = t * t * t;
+    t = t * (x as f64 + x as f64 + r) / (x as f64 + r + r);
+
+    /*
+     * Second step Newton iteration to 47 bits.  In double precision for
+     * efficiency and accuracy.
+     */
+    r = t * t * t;
+    t = t * (x as f64 + x as f64 + r) / (x as f64 + r + r);
+
+    /* rounding to 24 bits is perfect in round-to-nearest mode */
+    t as f32
+}

src/math/ceil.rs

+use core::f64;
+
+const TOINT: f64 = 1. / f64::EPSILON;
+
+#[inline]
+pub fn ceil(x: f64) -> f64 {
+    let u: u64 = x.to_bits();
+    let e: i64 = (u >> 52 & 0x7ff) as i64;
+    let y: f64;
+
+    if e >= 0x3ff + 52 || x == 0. {
+        return x;
+    }
+    // y = int(x) - x, where int(x) is an integer neighbor of x
+    y = if (u >> 63) != 0 {
+        x - TOINT + TOINT - x
+    } else {
+        x + TOINT - TOINT - x
+    };
+    // special case because of non-nearest rounding modes
+    if e <= 0x3ff - 1 {
+        force_eval!(y);
+        return if (u >> 63) != 0 { -0. } else { 1. };
+    }
+    if y < 0. {
+        x + y + 1.
+    } else {
+        x + y
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    #[test]
+    fn sanity_check() {
+        assert_eq!(super::ceil(1.1), 2.0);
+        assert_eq!(super::ceil(2.9), 3.0);
+    }
+}

src/math/ceilf.rs

+use core::f32;
+
+#[inline]
+pub fn ceilf(x: f32) -> f32 {
+    let mut ui = x.to_bits();
+    let e = (((ui >> 23) & 0xff) - 0x7f) as i32;
+
+    if e >= 23 {
+        return x;
+    }
+    if e >= 0 {
+        let m = 0x007fffff >> e;
+        if (ui & m) == 0 {
+            return x;
+        }
+        force_eval!(x + f32::from_bits(0x7b800000));
+        if ui >> 31 == 0 {
+            ui += m;
+        }
+        ui &= !m;
+    } else {
+        force_eval!(x + f32::from_bits(0x7b800000));
+        if ui >> 31 != 0 {
+            return -0.0;
+        } else if ui << 1 != 0 {
+            return 1.0;
+        }
+    }
+    return f32::from_bits(ui);
+}

src/math/cosf.rs

+use super::{k_cosf, k_sinf, rem_pio2f};
+
+use core::f64::consts::FRAC_PI_2;
+
+/* Small multiples of pi/2 rounded to double precision. */
+const C1_PIO2: f64 = 1. * FRAC_PI_2; /* 0x3FF921FB, 0x54442D18 */
+const C2_PIO2: f64 = 2. * FRAC_PI_2; /* 0x400921FB, 0x54442D18 */
+const C3_PIO2: f64 = 3. * FRAC_PI_2; /* 0x4012D97C, 0x7F3321D2 */
+const C4_PIO2: f64 = 4. * FRAC_PI_2; /* 0x401921FB, 0x54442D18 */
+
+#[inline]
+pub fn cosf(x: f32) -> f32 {
+    let x64 = x as f64;
+
+    let x1p120 = f32::from_bits(0x7b800000); // 0x1p120f === 2 ^ 120
+
+    let mut ix = x.to_bits();
+    let sign = (ix >> 31) != 0;
+    ix &= 0x7fffffff;
+
+    if ix <= 0x3f490fda {
+        /* |x| ~<= pi/4 */
+        if ix < 0x39800000 {
+            /* |x| < 2**-12 */
+            /* raise inexact if x != 0 */
+            force_eval!(x + x1p120);
+            return 1.;
+        }
+        return k_cosf(x64);
+    }
+    if ix <= 0x407b53d1 {
+        /* |x| ~<= 5*pi/4 */
+        if ix > 0x4016cbe3 {
+            /* |x|  ~> 3*pi/4 */
+            return -k_cosf(if sign { x64 + C2_PIO2 } else { x64 - C2_PIO2 });
+        } else {
+            if sign {
+                return k_sinf(x64 + C1_PIO2);
+            } else {
+                return k_sinf(C1_PIO2 - x64);
+            }
+        }
+    }
+    if ix <= 0x40e231d5 {
+        /* |x| ~<= 9*pi/4 */
+        if ix > 0x40afeddf {
+            /* |x| ~> 7*pi/4 */
+            return k_cosf(if sign { x64 + C4_PIO2 } else { x64 - C4_PIO2 });
+        } else {
+            if sign {
+                return k_sinf(-x64 - C3_PIO2);
+            } else {
+                return k_sinf(x64 - C3_PIO2);
+            }
+        }
+    }
+
+    /* cos(Inf or NaN) is NaN */
+    if ix >= 0x7f800000 {
+        return x - x;
+    }
+
+    /* general argument reduction needed */
+    let (n, y) = rem_pio2f(x);
+    match n & 3 {
+        0 => k_cosf(y),
+        1 => k_sinf(-y),
+        2 => -k_cosf(y),
+        _ => k_sinf(y),
+    }
+}

src/math/exp.rs

+/* origin: FreeBSD /usr/src/lib/msun/src/e_exp.c */
+/*
+ * ====================================================
+ * Copyright (C) 2004 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+/* exp(x)
+ * Returns the exponential of x.
+ *
+ * Method
+ *   1. Argument reduction:
+ *      Reduce x to an r so that |r| <= 0.5*ln2 ~ 0.34658.
+ *      Given x, find r and integer k such that
+ *
+ *               x = k*ln2 + r,  |r| <= 0.5*ln2.
+ *
+ *      Here r will be represented as r = hi-lo for better
+ *      accuracy.
+ *
+ *   2. Approximation of exp(r) by a special rational function on
+ *      the interval [0,0.34658]:
+ *      Write
+ *          R(r**2) = r*(exp(r)+1)/(exp(r)-1) = 2 + r*r/6 - r**4/360 + ...
+ *      We use a special Remez algorithm on [0,0.34658] to generate
+ *      a polynomial of degree 5 to approximate R. The maximum error
+ *      of this polynomial approximation is bounded by 2**-59. In
+ *      other words,
+ *          R(z) ~ 2.0 + P1*z + P2*z**2 + P3*z**3 + P4*z**4 + P5*z**5
+ *      (where z=r*r, and the values of P1 to P5 are listed below)
+ *      and
+ *          |                  5          |     -59
+ *          | 2.0+P1*z+...+P5*z   -  R(z) | <= 2
+ *          |                             |
+ *      The computation of exp(r) thus becomes
+ *                              2*r
+ *              exp(r) = 1 + ----------
+ *                            R(r) - r
+ *                                 r*c(r)
+ *                     = 1 + r + ----------- (for better accuracy)
+ *                                2 - c(r)
+ *      where
+ *                              2       4             10
+ *              c(r) = r - (P1*r  + P2*r  + ... + P5*r   ).
+ *
+ *   3. Scale back to obtain exp(x):
+ *      From step 1, we have
+ *         exp(x) = 2^k * exp(r)
+ *
+ * Special cases:
+ *      exp(INF) is INF, exp(NaN) is NaN;
+ *      exp(-INF) is 0, and
+ *      for finite argument, only exp(0)=1 is exact.
+ *
+ * Accuracy:
+ *      according to an error analysis, the error is always less than
+ *      1 ulp (unit in the last place).
+ *
+ * Misc. info.
+ *      For IEEE double
+ *          if x >  709.782712893383973096 then exp(x) overflows
+ *          if x < -745.133219101941108420 then exp(x) underflows
+ */
+
+use super::scalbn;
+
+const HALF: [f64; 2] = [0.5, -0.5];
+const LN2HI: f64 = 6.93147180369123816490e-01; /* 0x3fe62e42, 0xfee00000 */
+const LN2LO: f64 = 1.90821492927058770002e-10; /* 0x3dea39ef, 0x35793c76 */
+const INVLN2: f64 = 1.44269504088896338700e+00; /* 0x3ff71547, 0x652b82fe */
+const P1: f64 = 1.66666666666666019037e-01; /* 0x3FC55555, 0x5555553E */
+const P2: f64 = -2.77777777770155933842e-03; /* 0xBF66C16C, 0x16BEBD93 */
+const P3: f64 = 6.61375632143793436117e-05; /* 0x3F11566A, 0xAF25DE2C */
+const P4: f64 = -1.65339022054652515390e-06; /* 0xBEBBBD41, 0xC5D26BF1 */
+const P5: f64 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
+
+#[inline]
+pub fn exp(mut x: f64) -> f64 {
+    let x1p1023 = f64::from_bits(0x7fe0000000000000); // 0x1p1023 === 2 ^ 1023
+    let x1p_149 = f64::from_bits(0x36a0000000000000); // 0x1p-149 === 2 ^ -149
+
+    let hi: f64;
+    let lo: f64;
+    let c: f64;
+    let xx: f64;
+    let y: f64;
+    let k: i32;
+    let sign: i32;
+    let mut hx: u32;
+
+    hx = (x.to_bits() >> 32) as u32;
+    sign = (hx >> 31) as i32;
+    hx &= 0x7fffffff; /* high word of |x| */
+
+    /* special cases */
+    if hx >= 0x4086232b {
+        /* if |x| >= 708.39... */
+        if x.is_nan() {
+            return x;
+        }
+        if x > 709.782712893383973096 {
+            /* overflow if x!=inf */
+            x *= x1p1023;
+            return x;
+        }
+        if x < -708.39641853226410622 {
+            /* underflow if x!=-inf */
+            force_eval!((-x1p_149 / x) as f32);
+            if x < -745.13321910194110842 {
+                return 0.;
+            }
+        }
+    }
+
+    /* argument reduction */
+    if hx > 0x3fd62e42 {
+        /* if |x| > 0.5 ln2 */
+        if hx >= 0x3ff0a2b2 {
+            /* if |x| >= 1.5 ln2 */
+            k = (INVLN2 * x + HALF[sign as usize]) as i32;
+        } else {
+            k = 1 - sign - sign;
+        }
+        hi = x - k as f64 * LN2HI; /* k*ln2hi is exact here */
+        lo = k as f64 * LN2LO;
+        x = hi - lo;
+    } else if hx > 0x3e300000 {
+        /* if |x| > 2**-28 */
+        k = 0;
+        hi = x;
+        lo = 0.;
+    } else {
+        /* inexact if x!=0 */
+        force_eval!(x1p1023 + x);
+        return 1. + x;
+    }
+
+    /* x is now in primary range */
+    xx = x * x;
+    c = x - xx * (P1 + xx * (P2 + xx * (P3 + xx * (P4 + xx * P5))));
+    y = 1. + (x * c / (2. - c) - lo + hi);
+    if k == 0 {
+        y
+    } else {
+        scalbn(y, k)
+    }
+}

src/math/exp2.rs

+// origin: FreeBSD /usr/src/lib/msun/src/s_exp2.c */
+//-
+// Copyright (c) 2005 David Schultz <[e-mail]>
+// All rights reserved.
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions
+// are met:
+// 1. Redistributions of source code must retain the above copyright
+//    notice, this list of conditions and the following disclaimer.
+// 2. Redistributions in binary form must reproduce the above copyright
+//    notice, this list of conditions and the following disclaimer in the
+//    documentation and/or other materials provided with the distribution.
+//
+// THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+// FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+// DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+// OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+// HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+// LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+// OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+// SUCH DAMAGE.
+
+use super::scalbn::scalbn;
+
+const TBLSIZE: usize = 256;
+
+#[cfg_attr(rustfmt, rustfmt_skip)]
+static TBL: [u64; TBLSIZE * 2] = [
+    //  exp2(z + eps)          eps
+    0x3fe6a09e667f3d5d, 0x3d39880000000000,
+    0x3fe6b052fa751744, 0x3cd8000000000000,
+    0x3fe6c012750bd9fe, 0xbd28780000000000,
+    0x3fe6cfdcddd476bf, 0x3d1ec00000000000,
+    0x3fe6dfb23c651a29, 0xbcd8000000000000,
+    0x3fe6ef9298593ae3, 0xbcbc000000000000,
+    0x3fe6ff7df9519386, 0xbd2fd80000000000,
+    0x3fe70f7466f42da3, 0xbd2c880000000000,
+    0x3fe71f75e8ec5fc3, 0x3d13c00000000000,
+    0x3fe72f8286eacf05, 0xbd38300000000000,
+    0x3fe73f9a48a58152, 0xbd00c00000000000,
+    0x3fe74fbd35d7ccfc, 0x3d2f880000000000,
+    0x3fe75feb564267f1, 0x3d03e00000000000,
+    0x3fe77024b1ab6d48, 0xbd27d00000000000,
+    0x3fe780694fde5d38, 0xbcdd000000000000,
+    0x3fe790b938ac1d00, 0x3ce3000000000000,
+    0x3fe7a11473eb0178, 0xbced000000000000,
+    0x3fe7b17b0976d060, 0x3d20400000000000,
+    0x3fe7c1ed0130c133, 0x3ca0000000000000,
+    0x3fe7d26a62ff8636, 0xbd26900000000000,
+    0x3fe7e2f336cf4e3b, 0xbd02e00000000000,
+    0x3fe7f3878491c3e8, 0xbd24580000000000,
+    0x3fe80427543e1b4e, 0x3d33000000000000,
+    0x3fe814d2add1071a, 0x3d0f000000000000,
+    0x3fe82589994ccd7e, 0xbd21c00000000000,
+    0x3fe8364c1eb942d0, 0x3d29d00000000000,
+    0x3fe8471a4623cab5, 0x3d47100000000000,
+    0x3fe857f4179f5bbc, 0x3d22600000000000,
+    0x3fe868d99b4491af, 0xbd32c40000000000,
+    0x3fe879cad931a395, 0xbd23000000000000,
+    0x3fe88ac7d98a65b8, 0xbd2a800000000000,
+    0x3fe89bd0a4785800, 0xbced000000000000,
+    0x3fe8ace5422aa223, 0x3d33280000000000,
+    0x3fe8be05bad619fa, 0x3d42b40000000000,
+    0x3fe8cf3216b54383, 0xbd2ed00000000000,
+    0x3fe8e06a5e08664c, 0xbd20500000000000,
+    0x3fe8f1ae99157807, 0x3d28280000000000,
+    0x3fe902fed0282c0e, 0xbd1cb00000000000,
+    0x3fe9145b0b91ff96, 0xbd05e00000000000,
+    0x3fe925c353aa2ff9, 0x3cf5400000000000,
+    0x3fe93737b0cdc64a, 0x3d17200000000000,
+    0x3fe948b82b5f98ae, 0xbd09000000000000,
+    0x3fe95a44cbc852cb, 0x3d25680000000000,
+    0x3fe96bdd9a766f21, 0xbd36d00000000000,
+    0x3fe97d829fde4e2a, 0xbd01000000000000,
+    0x3fe98f33e47a23a3, 0x3d2d000000000000,
+    0x3fe9a0f170ca0604, 0xbd38a40000000000,
+    0x3fe9b2bb4d53ff89, 0x3d355c0000000000,
+    0x3fe9c49182a3f15b, 0x3d26b80000000000,
+    0x3fe9d674194bb8c5, 0xbcec000000000000,
+    0x3fe9e86319e3238e, 0x3d17d00000000000,
+    0x3fe9fa5e8d07f302, 0x3d16400000000000,
+    0x3fea0c667b5de54d, 0xbcf5000000000000,
+    0x3fea1e7aed8eb8f6, 0x3d09e00000000000,
+    0x3fea309bec4a2e27, 0x3d2ad80000000000,
+    0x3fea42c980460a5d, 0xbd1af00000000000,
+    0x3fea5503b23e259b, 0x3d0b600000000000,
+    0x3fea674a8af46213, 0x3d38880000000000,
+    0x3fea799e1330b3a7, 0x3d11200000000000,
+    0x3fea8bfe53c12e8d, 0x3d06c00000000000,
+    0x3fea9e6b5579fcd2, 0xbd29b80000000000,
+    0x3feab0e521356fb8, 0x3d2b700000000000,
+    0x3feac36bbfd3f381, 0x3cd9000000000000,
+    0x3fead5ff3a3c2780, 0x3ce4000000000000,
+    0x3feae89f995ad2a3, 0xbd2c900000000000,
+    0x3feafb4ce622f367, 0x3d16500000000000,
+    0x3feb0e07298db790, 0x3d2fd40000000000,
+    0x3feb20ce6c9a89a9, 0x3d12700000000000,
+    0x3feb33a2b84f1a4b, 0x3d4d470000000000,
+    0x3feb468415b747e7, 0xbd38380000000000,
+    0x3feb59728de5593a, 0x3c98000000000000,
+    0x3feb6c6e29f1c56a, 0x3d0ad00000000000,
+    0x3feb7f76f2fb5e50, 0x3cde800000000000,
+    0x3feb928cf22749b2, 0xbd04c00000000000,
+    0x3feba5b030a10603, 0xbd0d700000000000,
+    0x3febb8e0b79a6f66, 0x3d0d900000000000,
+    0x3febcc1e904bc1ff, 0x3d02a00000000000,
+    0x3febdf69c3f3a16f, 0xbd1f780000000000,
+    0x3febf2c25bd71db8, 0xbd10a00000000000,
+    0x3fec06286141b2e9, 0xbd11400000000000,
+    0x3fec199bdd8552e0, 0x3d0be00000000000,
+    0x3fec2d1cd9fa64ee, 0xbd09400000000000,
+    0x3fec40ab5fffd02f, 0xbd0ed00000000000,
+    0x3fec544778fafd15, 0x3d39660000000000,
+    0x3fec67f12e57d0cb, 0xbd1a100000000000,
+    0x3fec7ba88988c1b6, 0xbd58458000000000,
+    0x3fec8f6d9406e733, 0xbd1a480000000000,
+    0x3feca3405751c4df, 0x3ccb000000000000,
+    0x3fecb720dcef9094, 0x3d01400000000000,
+    0x3feccb0f2e6d1689, 0x3cf0200000000000,
+    0x3fecdf0b555dc412, 0x3cf3600000000000,
+    0x3fecf3155b5bab3b, 0xbd06900000000000,
+    0x3fed072d4a0789bc, 0x3d09a00000000000,
+    0x3fed1b532b08c8fa, 0xbd15e00000000000,
+    0x3fed2f87080d8a85, 0x3d1d280000000000,
+    0x3fed43c8eacaa203, 0x3d01a00000000000,
+    0x3fed5818dcfba491, 0x3cdf000000000000,
+    0x3fed6c76e862e6a1, 0xbd03a00000000000,
+    0x3fed80e316c9834e, 0xbd0cd80000000000,
+    0x3fed955d71ff6090, 0x3cf4c00000000000,
+    0x3feda9e603db32ae, 0x3cff900000000000,
+    0x3fedbe7cd63a8325, 0x3ce9800000000000,
+    0x3fedd321f301b445, 0xbcf5200000000000,
+    0x3fede7d5641c05bf, 0xbd1d700000000000,
+    0x3fedfc97337b9aec, 0xbd16140000000000,
+    0x3fee11676b197d5e, 0x3d0b480000000000,
+    0x3fee264614f5a3e7, 0x3d40ce0000000000,
+    0x3fee3b333b16ee5c, 0x3d0c680000000000,
+    0x3fee502ee78b3fb4, 0xbd09300000000000,
+    0x3fee653924676d68, 0xbce5000000000000,
+    0x3fee7a51fbc74c44, 0xbd07f80000000000,
+    0x3fee8f7977cdb726, 0xbcf3700000000000,
+    0x3feea4afa2a490e8, 0x3ce5d00000000000,
+    0x3feeb9f4867ccae4, 0x3d161a0000000000,
+    0x3feecf482d8e680d, 0x3cf5500000000000,
+    0x3feee4aaa2188514, 0x3cc6400000000000,
+    0x3feefa1bee615a13, 0xbcee800000000000,
+    0x3fef0f9c1cb64106, 0xbcfa880000000000,
+    0x3fef252b376bb963, 0xbd2c900000000000,
+    0x3fef3ac948dd7275, 0x3caa000000000000,
+    0x3fef50765b6e4524, 0xbcf4f00000000000,
+    0x3fef6632798844fd, 0x3cca800000000000,
+    0x3fef7bfdad9cbe38, 0x3cfabc0000000000,
+    0x3fef91d802243c82, 0xbcd4600000000000,
+    0x3fefa7c1819e908e, 0xbd0b0c0000000000,
+    0x3fefbdba3692d511, 0xbcc0e00000000000,
+    0x3fefd3c22b8f7194, 0xbd10de8000000000,
+    0x3fefe9d96b2a23ee, 0x3cee430000000000,
+    0x3ff0000000000000, 0x0,
+    0x3ff00b1afa5abcbe, 0xbcb3400000000000,
+    0x3ff0163da9fb3303, 0xbd12170000000000,
+    0x3ff02168143b0282, 0x3cba400000000000,
+    0x3ff02c9a3e77806c, 0x3cef980000000000,
+    0x3ff037d42e11bbca, 0xbcc7400000000000,
+    0x3ff04315e86e7f89, 0x3cd8300000000000,
+    0x3ff04e5f72f65467, 0xbd1a3f0000000000,
+    0x3ff059b0d315855a, 0xbd02840000000000,
+    0x3ff0650a0e3c1f95, 0x3cf1600000000000,
+    0x3ff0706b29ddf71a, 0x3d15240000000000,
+    0x3ff07bd42b72a82d, 0xbce9a00000000000,
+    0x3ff0874518759bd0, 0x3ce6400000000000,
+    0x3ff092bdf66607c8, 0xbd00780000000000,
+    0x3ff09e3ecac6f383, 0xbc98000000000000,
+    0x3ff0a9c79b1f3930, 0x3cffa00000000000,
+    0x3ff0b5586cf988fc, 0xbcfac80000000000,
+    0x3ff0c0f145e46c8a, 0x3cd9c00000000000,
+    0x3ff0cc922b724816, 0x3d05200000000000,
+    0x3ff0d83b23395dd8, 0xbcfad00000000000,
+    0x3ff0e3ec32d3d1f3, 0x3d1bac0000000000,
+    0x3ff0efa55fdfa9a6, 0xbd04e80000000000,
+    0x3ff0fb66affed2f0, 0xbd0d300000000000,
+    0x3ff1073028d7234b, 0x3cf1500000000000,
+    0x3ff11301d0125b5b, 0x3cec000000000000,
+    0x3ff11edbab5e2af9, 0x3d16bc0000000000,
+    0x3ff12abdc06c31d5, 0x3ce8400000000000,
+    0x3ff136a814f2047d, 0xbd0ed00000000000,
+    0x3ff1429aaea92de9, 0x3ce8e00000000000,
+    0x3ff14e95934f3138, 0x3ceb400000000000,
+    0x3ff15a98c8a58e71, 0x3d05300000000000,
+    0x3ff166a45471c3df, 0x3d03380000000000,
+    0x3ff172b83c7d5211, 0x3d28d40000000000,
+    0x3ff17ed48695bb9f, 0xbd05d00000000000,
+    0x3ff18af9388c8d93, 0xbd1c880000000000,
+    0x3ff1972658375d66, 0x3d11f00000000000,
+    0x3ff1a35beb6fcba7, 0x3d10480000000000,
+    0x3ff1af99f81387e3, 0xbd47390000000000,
+    0x3ff1bbe084045d54, 0x3d24e40000000000,
+    0x3ff1c82f95281c43, 0xbd0a200000000000,
+    0x3ff1d4873168b9b2, 0x3ce3800000000000,
+    0x3ff1e0e75eb44031, 0x3ceac00000000000,
+    0x3ff1ed5022fcd938, 0x3d01900000000000,
+    0x3ff1f9c18438cdf7, 0xbd1b780000000000,
+    0x3ff2063b88628d8f, 0x3d2d940000000000,
+    0x3ff212be3578a81e, 0x3cd8000000000000,
+    0x3ff21f49917ddd41, 0x3d2b340000000000,
+    0x3ff22bdda2791323, 0x3d19f80000000000,
+    0x3ff2387a6e7561e7, 0xbd19c80000000000,
+    0x3ff2451ffb821427, 0x3d02300000000000,
+    0x3ff251ce4fb2a602, 0xbd13480000000000,
+    0x3ff25e85711eceb0, 0x3d12700000000000,
+    0x3ff26b4565e27d16, 0x3d11d00000000000,
+    0x3ff2780e341de00f, 0x3d31ee0000000000,
+    0x3ff284dfe1f5633e, 0xbd14c00000000000,
+    0x3ff291ba7591bb30, 0xbd13d80000000000,
+    0x3ff29e9df51fdf09, 0x3d08b00000000000,
+    0x3ff2ab8a66d10e9b, 0xbd227c0000000000,
+    0x3ff2b87fd0dada3a, 0x3d2a340000000000,
+    0x3ff2c57e39771af9, 0xbd10800000000000,
+    0x3ff2d285a6e402d9, 0xbd0ed00000000000,
+    0x3ff2df961f641579, 0xbcf4200000000000,
+    0x3ff2ecafa93e2ecf, 0xbd24980000000000,
+    0x3ff2f9d24abd8822, 0xbd16300000000000,
+    0x3ff306fe0a31b625, 0xbd32360000000000,
+    0x3ff31432edeea50b, 0xbd70df8000000000,
+    0x3ff32170fc4cd7b8, 0xbd22480000000000,
+    0x3ff32eb83ba8e9a2, 0xbd25980000000000,
+    0x3ff33c08b2641766, 0x3d1ed00000000000,
+    0x3ff3496266e3fa27, 0xbcdc000000000000,
+    0x3ff356c55f929f0f, 0xbd30d80000000000,
+    0x3ff36431a2de88b9, 0x3d22c80000000000,
+    0x3ff371a7373aaa39, 0x3d20600000000000,
+    0x3ff37f26231e74fe, 0xbd16600000000000,
+    0x3ff38cae6d05d838, 0xbd0ae00000000000,
+    0x3ff39a401b713ec3, 0xbd44720000000000,
+    0x3ff3a7db34e5a020, 0x3d08200000000000,
+    0x3ff3b57fbfec6e95, 0x3d3e800000000000,
+    0x3ff3c32dc313a8f2, 0x3cef800000000000,
+    0x3ff3d0e544ede122, 0xbd17a00000000000,
+    0x3ff3dea64c1234bb, 0x3d26300000000000,
+    0x3ff3ec70df1c4ecc, 0xbd48a60000000000,
+    0x3ff3fa4504ac7e8c, 0xbd3cdc0000000000,
+    0x3ff40822c367a0bb, 0x3d25b80000000000,
+    0x3ff4160a21f72e95, 0x3d1ec00000000000,
+    0x3ff423fb27094646, 0xbd13600000000000,
+    0x3ff431f5d950a920, 0x3d23980000000000,
+    0x3ff43ffa3f84b9eb, 0x3cfa000000000000,
+    0x3ff44e0860618919, 0xbcf6c00000000000,
+    0x3ff45c2042a7d201, 0xbd0bc00000000000,
+    0x3ff46a41ed1d0016, 0xbd12800000000000,
+    0x3ff4786d668b3326, 0x3d30e00000000000,
+    0x3ff486a2b5c13c00, 0xbd2d400000000000,
+    0x3ff494e1e192af04, 0x3d0c200000000000,
+    0x3ff4a32af0d7d372, 0xbd1e500000000000,
+    0x3ff4b17dea6db801, 0x3d07800000000000,
+    0x3ff4bfdad53629e1, 0xbd13800000000000,
+    0x3ff4ce41b817c132, 0x3d00800000000000,
+    0x3ff4dcb299fddddb, 0x3d2c700000000000,
+    0x3ff4eb2d81d8ab96, 0xbd1ce00000000000,
+    0x3ff4f9b2769d2d02, 0x3d19200000000000,
+    0x3ff508417f4531c1, 0xbd08c00000000000,
+    0x3ff516daa2cf662a, 0xbcfa000000000000,
+    0x3ff5257de83f51ea, 0x3d4a080000000000,
+    0x3ff5342b569d4eda, 0xbd26d80000000000,
+    0x3ff542e2f4f6ac1a, 0xbd32440000000000,
+    0x3ff551a4ca5d94db, 0x3d483c0000000000,
+    0x3ff56070dde9116b, 0x3d24b00000000000,
+    0x3ff56f4736b529de, 0x3d415a0000000000,
+    0x3ff57e27dbe2c40e, 0xbd29e00000000000,
+    0x3ff58d12d497c76f, 0xbd23080000000000,
+    0x3ff59c0827ff0b4c, 0x3d4dec0000000000,
+    0x3ff5ab07dd485427, 0xbcc4000000000000,
+    0x3ff5ba11fba87af4, 0x3d30080000000000,
+    0x3ff5c9268a59460b, 0xbd26c80000000000,
+    0x3ff5d84590998e3f, 0x3d469a0000000000,
+    0x3ff5e76f15ad20e1, 0xbd1b400000000000,
+    0x3ff5f6a320dcebca, 0x3d17700000000000,
+    0x3ff605e1b976dcb8, 0x3d26f80000000000,
+    0x3ff6152ae6cdf715, 0x3d01000000000000,
+    0x3ff6247eb03a5531, 0xbd15d00000000000,
+    0x3ff633dd1d1929b5, 0xbd12d00000000000,
+    0x3ff6434634ccc313, 0xbcea800000000000,
+    0x3ff652b9febc8efa, 0xbd28600000000000,
+    0x3ff6623882553397, 0x3d71fe0000000000,
+    0x3ff671c1c708328e, 0xbd37200000000000,
+    0x3ff68155d44ca97e, 0x3ce6800000000000,
+    0x3ff690f4b19e9471, 0xbd29780000000000,
+];
+
+// exp2(x): compute the base 2 exponential of x
+//
+// Accuracy: Peak error < 0.503 ulp for normalized results.
+//
+// Method: (accurate tables)
+//
+//   Reduce x:
+//     x = k + y, for integer k and |y| <= 1/2.
+//     Thus we have exp2(x) = 2**k * exp2(y).
+//
+//   Reduce y:
+//     y = i/TBLSIZE + z - eps[i] for integer i near y * TBLSIZE.
+//     Thus we have exp2(y) = exp2(i/TBLSIZE) * exp2(z - eps[i]),
+//     with |z - eps[i]| <= 2**-9 + 2**-39 for the table used.
+//
+//   We compute exp2(i/TBLSIZE) via table lookup and exp2(z - eps[i]) via
+//   a degree-5 minimax polynomial with maximum error under 1.3 * 2**-61.
+//   The values in exp2t[] and eps[] are chosen such that
+//   exp2t[i] = exp2(i/TBLSIZE + eps[i]), and eps[i] is a small offset such
+//   that exp2t[i] is accurate to 2**-64.
+//
+//   Note that the range of i is +-TBLSIZE/2, so we actually index the tables
+//   by i0 = i + TBLSIZE/2.  For cache efficiency, exp2t[] and eps[] are
+//   virtual tables, interleaved in the real table tbl[].
+//
+//   This method is due to Gal, with many details due to Gal and Bachelis:
+//
+//      Gal, S. and Bachelis, B.  An Accurate Elementary Mathematical Library
+//      for the IEEE Floating Point Standard.  TOMS 17(1), 26-46 (1991).
+#[inline]
+pub fn exp2(mut x: f64) -> f64 {
+    let redux = f64::from_bits(0x4338000000000000) / TBLSIZE as f64;
+    let p1 = f64::from_bits(0x3fe62e42fefa39ef);
+    let p2 = f64::from_bits(0x3fcebfbdff82c575);
+    let p3 = f64::from_bits(0x3fac6b08d704a0a6);
+    let p4 = f64::from_bits(0x3f83b2ab88f70400);
+    let p5 = f64::from_bits(0x3f55d88003875c74);
+
+    // double_t r, t, z;
+    // uint32_t ix, i0;
+    // union {double f; uint64_t i;} u = {x};
+    // union {uint32_t u; int32_t i;} k;
+    let x1p1023 = f64::from_bits(0x7fe0000000000000);
+    let x1p52 = f64::from_bits(0x4330000000000000);
+    let _0x1p_149 = f64::from_bits(0xb6a0000000000000);
+
+    /* Filter out exceptional cases. */
+    let ui = f64::to_bits(x);
+    let ix = ui >> 32 & 0x7fffffff;
+    if ix >= 0x408ff000 {
+        /* |x| >= 1022 or nan */
+        if ix >= 0x40900000 && ui >> 63 == 0 {
+            /* x >= 1024 or nan */
+            /* overflow */
+            x *= x1p1023;
+            return x;
+        }
+        if ix >= 0x7ff00000 {
+            /* -inf or -nan */
+            return -1.0 / x;
+        }
+        if ui >> 63 != 0 {
+            /* x <= -1022 */
+            /* underflow */
+            if x <= -1075.0 || x - x1p52 + x1p52 != x {
+                force_eval!((_0x1p_149 / x) as f32);
+            }
+            if x <= -1075.0 {
+                return 0.0;
+            }
+        }
+    } else if ix < 0x3c900000 {
+        /* |x| < 0x1p-54 */
+        return 1.0 + x;
+    }
+
+    /* Reduce x, computing z, i0, and k. */
+    let ui = f64::to_bits(x + redux);
+    let mut i0 = ui as u32;
+    i0 += TBLSIZE as u32 / 2;
+    let ku = i0 / TBLSIZE as u32 * TBLSIZE as u32;
+    let ki = ku as i32 / TBLSIZE as i32;
+    i0 %= TBLSIZE as u32;
+    let uf = f64::from_bits(ui) - redux;
+    let mut z = x - uf;
+
+    /* Compute r = exp2(y) = exp2t[i0] * p(z - eps[i]). */
+    let t = f64::from_bits(TBL[2 * i0 as usize]); /* exp2t[i0] */
+    z -= f64::from_bits(TBL[2 * i0 as usize + 1]); /* eps[i0]   */
+    let r = t + t * z * (p1 + z * (p2 + z * (p3 + z * (p4 + z * p5))));
+
+    scalbn(r, ki)
+}

src/math/exp2f.rs

+// origin: FreeBSD /usr/src/lib/msun/src/s_exp2f.c
+//-
+// Copyright (c) 2005 David Schultz <[e-mail]>
+// All rights reserved.
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions
+// are met:
+// 1. Redistributions of source code must retain the above copyright
+//    notice, this list of conditions and the following disclaimer.
+// 2. Redistributions in binary form must reproduce the above copyright
+//    notice, this list of conditions and the following disclaimer in the
+//    documentation and/or other materials provided with the distribution.
+//
+// THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
+// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
+// FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+// DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
+// OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
+// HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
+// LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
+// OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
+// SUCH DAMAGE.
+
+const TBLSIZE: usize = 16;
+
+static EXP2FT: [u64; TBLSIZE] = [
+    0x3fe6a09e667f3bcd,
+    0x3fe7a11473eb0187,
+    0x3fe8ace5422aa0db,
+    0x3fe9c49182a3f090,
+    0x3feae89f995ad3ad,
+    0x3fec199bdd85529c,
+    0x3fed5818dcfba487,
+    0x3feea4afa2a490da,
+    0x3ff0000000000000,
+    0x3ff0b5586cf9890f,
+    0x3ff172b83c7d517b,
+    0x3ff2387a6e756238,
+    0x3ff306fe0a31b715,
+    0x3ff3dea64c123422,
+    0x3ff4bfdad5362a27,
+    0x3ff5ab07dd485429,
+];
+
+// exp2f(x): compute the base 2 exponential of x
+//
+// Accuracy: Peak error < 0.501 ulp; location of peak: -0.030110927.
+//
+// Method: (equally-spaced tables)
+//
+//   Reduce x:
+//     x = k + y, for integer k and |y| <= 1/2.
+//     Thus we have exp2f(x) = 2**k * exp2(y).
+//
+//   Reduce y:
+//     y = i/TBLSIZE + z for integer i near y * TBLSIZE.
+//     Thus we have exp2(y) = exp2(i/TBLSIZE) * exp2(z),
+//     with |z| <= 2**-(TBLSIZE+1).
+//
+//   We compute exp2(i/TBLSIZE) via table lookup and exp2(z) via a
+//   degree-4 minimax polynomial with maximum error under 1.4 * 2**-33.
+//   Using double precision for everything except the reduction makes
+//   roundoff error insignificant and simplifies the scaling step.
+//
+//   This method is due to Tang, but I do not use his suggested parameters:
+//
+//      Tang, P.  Table-driven Implementation of the Exponential Function
+//      in IEEE Floating-Point Arithmetic.  TOMS 15(2), 144-157 (1989).
+#[inline]
+pub fn exp2f(mut x: f32) -> f32 {
+    let redux = f32::from_bits(0x4b400000) / TBLSIZE as f32;
+    let p1 = f32::from_bits(0x3f317218);
+    let p2 = f32::from_bits(0x3e75fdf0);
+    let p3 = f32::from_bits(0x3d6359a4);
+    let p4 = f32::from_bits(0x3c1d964e);
+
+    // double_t t, r, z;
+    // uint32_t ix, i0, k;
+
+    let x1p127 = f32::from_bits(0x7f000000);
+
+    /* Filter out exceptional cases. */
+    let ui = f32::to_bits(x);
+    let ix = ui & 0x7fffffff;
+    if ix > 0x42fc0000 {
+        /* |x| > 126 */
+        if ix > 0x7f800000 {
+            /* NaN */
+            return x;
+        }
+        if ui >= 0x43000000 && ui < 0x80000000 {
+            /* x >= 128 */
+            x *= x1p127;
+            return x;
+        }
+        if ui >= 0x80000000 {
+            /* x < -126 */
+            if ui >= 0xc3160000 || (ui & 0x0000ffff != 0) {
+                force_eval!(f32::from_bits(0x80000001) / x);
+            }
+            if ui >= 0xc3160000 {
+                /* x <= -150 */
+                return 0.0;
+            }
+        }
+    } else if ix <= 0x33000000 {
+        /* |x| <= 0x1p-25 */
+        return 1.0 + x;
+    }
+
+    /* Reduce x, computing z, i0, and k. */
+    let ui = f32::to_bits(x + redux);
+    let mut i0 = ui;
+    i0 += TBLSIZE as u32 / 2;
+    let k = i0 / TBLSIZE as u32;
+    let ukf = f64::from_bits(((0x3ff + k) as u64) << 52);
+    i0 &= TBLSIZE as u32 - 1;
+    let mut uf = f32::from_bits(ui);
+    uf -= redux;
+    let z: f64 = (x - uf) as f64;
+    /* Compute r = exp2(y) = exp2ft[i0] * p(z). */
+    let r: f64 = f64::from_bits(EXP2FT[i0 as usize]);
+    let t: f64 = r as f64 * z;
+    let r: f64 = r + t * (p1 as f64 + z * p2 as f64) + t * (z * z) * (p3 as f64 + z * p4 as f64);
+
+    /* Scale by 2**k */
+    (r * ukf) as f32
+}

src/math/expm1.rs

+use core::f64;
+
+const O_THRESHOLD: f64 = 7.09782712893383973096e+02; /* 0x40862E42, 0xFEFA39EF */
+const LN2_HI: f64 = 6.93147180369123816490e-01; /* 0x3fe62e42, 0xfee00000 */
+const LN2_LO: f64 = 1.90821492927058770002e-10; /* 0x3dea39ef, 0x35793c76 */
+const INVLN2: f64 = 1.44269504088896338700e+00; /* 0x3ff71547, 0x652b82fe */
+/* Scaled Q's: Qn_here = 2**n * Qn_above, for R(2*z) where z = hxs = x*x/2: */
+const Q1: f64 = -3.33333333333331316428e-02; /* BFA11111 111110F4 */
+const Q2: f64 = 1.58730158725481460165e-03; /* 3F5A01A0 19FE5585 */
+const Q3: f64 = -7.93650757867487942473e-05; /* BF14CE19 9EAADBB7 */
+const Q4: f64 = 4.00821782732936239552e-06; /* 3ED0CFCA 86E65239 */
+const Q5: f64 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
+
+#[allow(warnings)]
+pub fn expm1(mut x: f64) -> f64 {
+    let hi: f64;
+    let lo: f64;
+    let k: i32;
+    let c: f64;
+    let mut t: f64;
+    let mut y: f64;
+
+    let mut ui = x.to_bits();
+    let hx = ((ui >> 32) & 0x7fffffff) as u32;
+    let sign = (ui >> 63) as i32;
+
+    /* filter out huge and non-finite argument */
+    if hx >= 0x4043687A {
+        /* if |x|>=56*ln2 */
+        if x.is_nan() {
+            return x;
+        }
+        if sign != 0 {
+            return -1.0;
+        }
+        if x > O_THRESHOLD {
+            x *= f64::from_bits(0x7fe0000000000000);
+            return x;
+        }
+    }
+
+    /* argument reduction */
+    if hx > 0x3fd62e42 {
+        /* if  |x| > 0.5 ln2 */
+        if hx < 0x3FF0A2B2 {
+            /* and |x| < 1.5 ln2 */
+            if sign == 0 {
+                hi = x - LN2_HI;
+                lo = LN2_LO;
+                k = 1;
+            } else {
+                hi = x + LN2_HI;
+                lo = -LN2_LO;
+                k = -1;
+            }
+        } else {
+            k = (INVLN2 * x + if sign != 0 { -0.5 } else { 0.5 }) as i32;
+            t = k as f64;
+            hi = x - t * LN2_HI; /* t*ln2_hi is exact here */
+            lo = t * LN2_LO;
+        }
+        x = hi - lo;
+        c = (hi - x) - lo;
+    } else if hx < 0x3c900000 {
+        /* |x| < 2**-54, return x */
+        if hx < 0x00100000 {
+            force_eval!(x);
+        }
+        return x;
+    } else {
+        c = 0.0;
+        k = 0;
+    }
+
+    /* x is now in primary range */
+    let hfx = 0.5 * x;
+    let hxs = x * hfx;
+    let r1 = 1.0 + hxs * (Q1 + hxs * (Q2 + hxs * (Q3 + hxs * (Q4 + hxs * Q5))));
+    t = 3.0 - r1 * hfx;
+    let mut e = hxs * ((r1 - t) / (6.0 - x * t));
+    if k == 0 {
+        /* c is 0 */
+        return x - (x * e - hxs);
+    }
+    e = x * (e - c) - c;
+    e -= hxs;
+    /* exp(x) ~ 2^k (x_reduced - e + 1) */
+    if k == -1 {
+        return 0.5 * (x - e) - 0.5;
+    }
+    if k == 1 {
+        if x < -0.25 {
+            return -2.0 * (e - (x + 0.5));
+        }
+        return 1.0 + 2.0 * (x - e);
+    }
+    ui = ((0x3ff + k) as u64) << 52; /* 2^k */
+    let twopk = f64::from_bits(ui);
+    if k < 0 || k > 56 {
+        /* suffice to return exp(x)-1 */
+        y = x - e + 1.0;
+        if k == 1024 {
+            y = y * 2.0 * f64::from_bits(0x7fe0000000000000);
+        } else {
+            y = y * twopk;
+        }
+        return y - 1.0;
+    }
+    ui = ((0x3ff - k) as u64) << 52; /* 2^-k */
+    let uf = f64::from_bits(ui);
+    if k < 20 {
+        y = (x - e + (1.0 - uf)) * twopk;
+    } else {
+        y = (x - (e + uf) + 1.0) * twopk;
+    }
+    y
+}
+
+#[cfg(test)]
+mod tests {
+    #[test]
+    fn sanity_check() {
+        assert_eq!(super::expm1(1.1), 2.0041660239464334);
+    }
+}

src/math/expm1f.rs

+const O_THRESHOLD: f32 = 8.8721679688e+01; /* 0x42b17180 */
+const LN2_HI: f32 = 6.9313812256e-01; /* 0x3f317180 */
+const LN2_LO: f32 = 9.0580006145e-06; /* 0x3717f7d1 */
+const INV_LN2: f32 = 1.4426950216e+00; /* 0x3fb8aa3b */
+/*
+ * Domain [-0.34568, 0.34568], range ~[-6.694e-10, 6.696e-10]:
+ * |6 / x * (1 + 2 * (1 / (exp(x) - 1) - 1 / x)) - q(x)| < 2**-30.04
+ * Scaled coefficients: Qn_here = 2**n * Qn_for_q (see s_expm1.c):
+ */
+const Q1: f32 = -3.3333212137e-2; /* -0x888868.0p-28 */
+const Q2: f32 = 1.5807170421e-3; /*  0xcf3010.0p-33 */
+
+#[inline]
+pub fn expm1f(mut x: f32) -> f32 {
+    let x1p127 = f32::from_bits(0x7f000000); // 0x1p127f === 2 ^ 127
+
+    let mut hx = x.to_bits();
+    let sign = (hx >> 31) != 0;
+    hx &= 0x7fffffff;
+
+    /* filter out huge and non-finite argument */
+    if hx >= 0x4195b844 {
+        /* if |x|>=27*ln2 */
+        if hx > 0x7f800000 {
+            /* NaN */
+            return x;
+        }
+        if sign {
+            return -1.;
+        }
+        if x > O_THRESHOLD {
+            x *= x1p127;
+            return x;
+        }
+    }
+
+    let k: i32;
+    let hi: f32;
+    let lo: f32;
+    let mut c = 0f32;
+    /* argument reduction */
+    if hx > 0x3eb17218 {
+        /* if  |x| > 0.5 ln2 */
+        if hx < 0x3F851592 {
+            /* and |x| < 1.5 ln2 */
+            if !sign {
+                hi = x - LN2_HI;
+                lo = LN2_LO;
+                k = 1;
+            } else {
+                hi = x + LN2_HI;
+                lo = -LN2_LO;
+                k = -1;
+            }
+        } else {
+            k = (INV_LN2 * x + (if sign { -0.5 } else { 0.5 })) as i32;
+            let t = k as f32;
+            hi = x - t * LN2_HI; /* t*ln2_hi is exact here */
+            lo = t * LN2_LO;
+        }
+        x = hi - lo;
+        c = (hi - x) - lo;
+    } else if hx < 0x33000000 {
+        /* when |x|<2**-25, return x */
+        if hx < 0x00800000 {
+            force_eval!(x * x);
+        }
+        return x;
+    } else {
+        k = 0;
+    }
+
+    /* x is now in primary range */
+    let hfx = 0.5 * x;
+    let hxs = x * hfx;
+    let r1 = 1. + hxs * (Q1 + hxs * Q2);
+    let t = 3. - r1 * hfx;
+    let mut e = hxs * ((r1 - t) / (6. - x * t));
+    if k == 0 {
+        /* c is 0 */
+        return x - (x * e - hxs);
+    }
+    e = x * (e - c) - c;
+    e -= hxs;
+    /* exp(x) ~ 2^k (x_reduced - e + 1) */
+    if k == -1 {
+        return 0.5 * (x - e) - 0.5;
+    }
+    if k == 1 {
+        if x < -0.25 {
+            return -2. * (e - (x + 0.5));
+        }
+        return 1. + 2. * (x - e);
+    }
+    let twopk = f32::from_bits(((0x7f + k) << 23) as u32); /* 2^k */
+    if (k < 0) || (k > 56) {
+        /* suffice to return exp(x)-1 */
+        let mut y = x - e + 1.;
+        if k == 128 {
+            y = y * 2. * x1p127;
+        } else {
+            y = y * twopk;
+        }
+        return y - 1.;
+    }
+    let uf = f32::from_bits(((0x7f - k) << 23) as u32); /* 2^-k */
+    if k < 23 {
+        (x - e + (1. - uf)) * twopk
+    } else {
+        (x - (e + uf) + 1.) * twopk
+    }
+}

src/math/fdim.rs

+use core::f64;
+
+#[inline]
+pub fn fdim(x: f64, y: f64) -> f64 {
+    if x.is_nan() {
+        x
+    } else if y.is_nan() {
+        y
+    } else {
+        if x > y {
+            x - y
+        } else {
+            0.0
+        }
+    }
+}

src/math/fdimf.rs

+use core::f32;
+
+#[inline]
+pub fn fdimf(x: f32, y: f32) -> f32 {
+    if x.is_nan() {
+        x
+    } else if y.is_nan() {
+        y
+    } else {
+        if x > y {
+            x - y
+        } else {
+            0.0
+        }
+    }
+}

src/math/floorf.rs

+use core::f32;
+
+#[inline]
+pub fn floorf(x: f32) -> f32 {
+    let mut ui = x.to_bits();
+    let e = (((ui >> 23) & 0xff) - 0x7f) as i32;
+
+    if e >= 23 {
+        return x;
+    }
+    if e >= 0 {
+        let m: u32 = 0x007fffff >> e;
+        if (ui & m) == 0 {
+            return x;
+        }
+        force_eval!(x + f32::from_bits(0x7b800000));
+        if ui >> 31 != 0 {
+            ui += m;
+        }
+        ui &= !m;
+    } else {
+        force_eval!(x + f32::from_bits(0x7b800000));
+        if ui >> 31 == 0 {
+            ui = 0;
+        } else if ui << 1 != 0 {
+            return -1.0;
+        }
+    }
+    return f32::from_bits(ui);
+}

src/math/fma.rs

+use core::{f32, f64};
+
+use super::scalbn;
+
+const ZEROINFNAN: i32 = 0x7ff - 0x3ff - 52 - 1;
+
+struct Num {
+    m: u64,
+    e: i32,
+    sign: i32,
+}
+
+#[inline]
+fn normalize(x: f64) -> Num {
+    let x1p63: f64 = f64::from_bits(0x43e0000000000000); // 0x1p63 === 2 ^ 63
+
+    let mut ix: u64 = x.to_bits();
+    let mut e: i32 = (ix >> 52) as i32;
+    let sign: i32 = e & 0x800;
+    e &= 0x7ff;
+    if e == 0 {
+        ix = (x * x1p63).to_bits();
+        e = (ix >> 52) as i32 & 0x7ff;
+        e = if e != 0 { e - 63 } else { 0x800 };
+    }
+    ix &= (1 << 52) - 1;
+    ix |= 1 << 52;
+    ix <<= 1;
+    e -= 0x3ff + 52 + 1;
+    Num { m: ix, e, sign }
+}
+
+#[inline]
+fn mul(x: u64, y: u64) -> (u64, u64) {
+    let t1: u64;
+    let t2: u64;
+    let t3: u64;
+    let xlo: u64 = x as u32 as u64;
+    let xhi: u64 = x >> 32;
+    let ylo: u64 = y as u32 as u64;
+    let yhi: u64 = y >> 32;
+
+    t1 = xlo * ylo;
+    t2 = xlo * yhi + xhi * ylo;
+    t3 = xhi * yhi;
+    let lo = t1 + (t2 << 32);
+    let hi = t3 + (t2 >> 32) + (t1 > lo) as u64;
+    (hi, lo)
+}
+
+#[inline]
+pub fn fma(x: f64, y: f64, z: f64) -> f64 {
+    let x1p63: f64 = f64::from_bits(0x43e0000000000000); // 0x1p63 === 2 ^ 63
+    let x0_ffffff8p_63 = f64::from_bits(0x3bfffffff0000000); // 0x0.ffffff8p-63
+
+    /* normalize so top 10bits and last bit are 0 */
+    let nx = normalize(x);
+    let ny = normalize(y);
+    let nz = normalize(z);
+
+    if nx.e >= ZEROINFNAN || ny.e >= ZEROINFNAN {
+        return x * y + z;
+    }
+    if nz.e >= ZEROINFNAN {
+        if nz.e > ZEROINFNAN {
+            /* z==0 */
+            return x * y + z;
+        }
+        return z;
+    }
+
+    /* mul: r = x*y */
+    let zhi: u64;
+    let zlo: u64;
+    let (mut rhi, mut rlo) = mul(nx.m, ny.m);
+    /* either top 20 or 21 bits of rhi and last 2 bits of rlo are 0 */
+
+    /* align exponents */
+    let mut e: i32 = nx.e + ny.e;
+    let mut d: i32 = nz.e - e;
+    /* shift bits z<<=kz, r>>=kr, so kz+kr == d, set e = e+kr (== ez-kz) */
+    if d > 0 {
+        if d < 64 {
+            zlo = nz.m << d;
+            zhi = nz.m >> 64 - d;
+        } else {
+            zlo = 0;
+            zhi = nz.m;
+            e = nz.e - 64;
+            d -= 64;
+            if d == 0 {
+            } else if d < 64 {
+                rlo = rhi << 64 - d | rlo >> d | ((rlo << 64 - d) != 0) as u64;
+                rhi = rhi >> d;
+            } else {
+                rlo = 1;
+                rhi = 0;
+            }
+        }
+    } else {
+        zhi = 0;
+        d = -d;
+        if d == 0 {
+            zlo = nz.m;
+        } else if d < 64 {
+            zlo = nz.m >> d | ((nz.m << 64 - d) != 0) as u64;
+        } else {
+            zlo = 1;
+        }
+    }
+
+    /* add */
+    let mut sign: i32 = nx.sign ^ ny.sign;
+    let samesign: bool = (sign ^ nz.sign) == 0;
+    let mut nonzero: i32 = 1;
+    if samesign {
+        /* r += z */
+        rlo += zlo;
+        rhi += zhi + (rlo < zlo) as u64;
+    } else {
+        /* r -= z */
+        let t = rlo;
+        rlo -= zlo;
+        rhi = rhi - zhi - (t < rlo) as u64;
+        if (rhi >> 63) != 0 {
+            rlo = (-(rlo as i64)) as u64;
+            rhi = (-(rhi as i64)) as u64 - (rlo != 0) as u64;
+            sign = (sign == 0) as i32;
+        }
+        nonzero = (rhi != 0) as i32;
+    }
+
+    /* set rhi to top 63bit of the result (last bit is sticky) */
+    if nonzero != 0 {
+        e += 64;
+        d = rhi.leading_zeros() as i32 - 1;
+        /* note: d > 0 */
+        rhi = rhi << d | rlo >> 64 - d | ((rlo << d) != 0) as u64;
+    } else if rlo != 0 {
+        d = rlo.leading_zeros() as i32 - 1;
+        if d < 0 {
+            rhi = rlo >> 1 | (rlo & 1);
+        } else {
+            rhi = rlo << d;
+        }
+    } else {
+        /* exact +-0 */
+        return x * y + z;
+    }
+    e -= d;
+
+    /* convert to double */
+    let mut i: i64 = rhi as i64; /* i is in [1<<62,(1<<63)-1] */
+    if sign != 0 {
+        i = -i;
+    }
+    let mut r: f64 = i as f64; /* |r| is in [0x1p62,0x1p63] */
+
+    if e < -1022 - 62 {
+        /* result is subnormal before rounding */
+        if e == -1022 - 63 {
+            let mut c: f64 = x1p63;
+            if sign != 0 {
+                c = -c;
+            }
+            if r == c {
+                /* min normal after rounding, underflow depends
+                   on arch behaviour which can be imitated by
+                   a double to float conversion */
+                let fltmin: f32 = (x0_ffffff8p_63 * f32::MIN_POSITIVE as f64 * r) as f32;
+                return f64::MIN_POSITIVE / f32::MIN_POSITIVE as f64 * fltmin as f64;
+            }
+            /* one bit is lost when scaled, add another top bit to
+               only round once at conversion if it is inexact */
+            if (rhi << 53) != 0 {
+                i = (rhi >> 1 | (rhi & 1) | 1 << 62) as i64;
+                if sign != 0 {
+                    i = -i;
+                }
+                r = i as f64;
+                r = 2. * r - c; /* remove top bit */
+
+                /* raise underflow portably, such that it
+                   cannot be optimized away */
+                {
+                    let tiny: f64 = f64::MIN_POSITIVE / f32::MIN_POSITIVE as f64 * r;
+                    r += (tiny * tiny) * (r - r);
+                }
+            }
+        } else {
+            /* only round once when scaled */
+            d = 10;
+            i = ((rhi >> d | ((rhi << 64 - d) != 0) as u64) << d) as i64;
+            if sign != 0 {
+                i = -i;
+            }
+            r = i as f64;
+        }
+    }
+    scalbn(r, e)
+}

src/math/fmod.rs

+use core::u64;
+
+#[inline]
+pub fn fmod(x: f64, y: f64) -> f64 {
+    let mut uxi = x.to_bits();
+    let mut uyi = y.to_bits();
+    let mut ex = (uxi >> 52 & 0x7ff) as i64;
+    let mut ey = (uyi >> 52 & 0x7ff) as i64;
+    let sx = uxi >> 63;
+    let mut i;
+
+    if uyi << 1 == 0 || y.is_nan() || ex == 0x7ff {
+        return (x * y) / (x * y);
+    }
+    if uxi << 1 <= uyi << 1 {
+        if uxi << 1 == uyi << 1 {
+            return 0.0 * x;
+        }
+        return x;
+    }
+
+    /* normalize x and y */
+    if ex == 0 {
+        i = uxi << 12;
+        while i >> 63 == 0 {
+            ex -= 1;
+            i <<= 1;
+        }
+        uxi <<= -ex + 1;
+    } else {
+        uxi &= u64::MAX >> 12;
+        uxi |= 1 << 52;
+    }
+    if ey == 0 {
+        i = uyi << 12;
+        while i >> 63 == 0 {
+            ey -= 1;
+            i <<= 1;
+        }
+        uyi <<= -ey + 1;
+    } else {
+        uyi &= u64::MAX >> 12;
+        uyi |= 1 << 52;
+    }
+
+    /* x mod y */
+    while ex > ey {
+        i = uxi - uyi;
+        if i >> 63 == 0 {
+            if i == 0 {
+                return 0.0 * x;
+            }
+            uxi = i;
+        }
+        uxi <<= 1;
+        ex -= 1;
+    }
+    i = uxi - uyi;
+    if i >> 63 == 0 {
+        if i == 0 {
+            return 0.0 * x;
+        }
+        uxi = i;
+    }
+    while uxi >> 52 == 0 {
+        uxi <<= 1;
+        ex -= 1;
+    }
+
+    /* scale result */
+    if ex > 0 {
+        uxi -= 1 << 52;
+        uxi |= (ex as u64) << 52;
+    } else {
+        uxi >>= -ex + 1;
+    }
+    uxi |= (sx as u64) << 63;
+
+    f64::from_bits(uxi)
+}

src/math/k_cosf.rs

+/* |cos(x) - c(x)| < 2**-34.1 (~[-5.37e-11, 5.295e-11]). */
+const C0: f64 = -0.499999997251031003120; /* -0x1ffffffd0c5e81.0p-54 */
+const C1: f64 = 0.0416666233237390631894; /*  0x155553e1053a42.0p-57 */
+const C2: f64 = -0.00138867637746099294692; /* -0x16c087e80f1e27.0p-62 */
+const C3: f64 = 0.0000243904487962774090654; /*  0x199342e0ee5069.0p-68 */
+
+#[inline]
+pub(crate) fn k_cosf(x: f64) -> f32 {
+    let z = x * x;
+    let w = z * z;
+    let r = C2 + z * C3;
+    (((1.0 + z * C0) + w * C1) + (w * z) * r) as f32
+}

src/math/k_sinf.rs

+/* |sin(x)/x - s(x)| < 2**-37.5 (~[-4.89e-12, 4.824e-12]). */
+const S1: f64 = -0.166666666416265235595; /* -0x15555554cbac77.0p-55 */
+const S2: f64 = 0.0083333293858894631756; /*  0x111110896efbb2.0p-59 */
+const S3: f64 = -0.000198393348360966317347; /* -0x1a00f9e2cae774.0p-65 */
+const S4: f64 = 0.0000027183114939898219064; /*  0x16cd878c3b46a7.0p-71 */
+
+#[inline]
+pub(crate) fn k_sinf(x: f64) -> f32 {
+    let z = x * x;
+    let w = z * z;
+    let r = S3 + z * S4;
+    let s = z * x;
+    ((x + s * (S1 + z * S2)) + s * w * r) as f32
+}

src/math/k_tanf.rs

+/* |tan(x)/x - t(x)| < 2**-25.5 (~[-2e-08, 2e-08]). */
+const T: [f64; 6] = [
+    0.333331395030791399758,   /* 0x15554d3418c99f.0p-54 */
+    0.133392002712976742718,   /* 0x1112fd38999f72.0p-55 */
+    0.0533812378445670393523,  /* 0x1b54c91d865afe.0p-57 */
+    0.0245283181166547278873,  /* 0x191df3908c33ce.0p-58 */
+    0.00297435743359967304927, /* 0x185dadfcecf44e.0p-61 */
+    0.00946564784943673166728, /* 0x1362b9bf971bcd.0p-59 */
+];
+
+#[inline]
+pub(crate) fn k_tanf(x: f64, odd: bool) -> f32 {
+    let z = x * x;
+    /*
+     * Split up the polynomial into small independent terms to give
+     * opportunities for parallel evaluation.  The chosen splitting is
+     * micro-optimized for Athlons (XP, X64).  It costs 2 multiplications
+     * relative to Horner's method on sequential machines.
+     *
+     * We add the small terms from lowest degree up for efficiency on
+     * non-sequential machines (the lowest degree terms tend to be ready
+     * earlier).  Apart from this, we don't care about order of
+     * operations, and don't need to to care since we have precision to
+     * spare.  However, the chosen splitting is good for accuracy too,
+     * and would give results as accurate as Horner's method if the
+     * small terms were added from highest degree down.
+     */
+    let mut r = T[4] + z * T[5];
+    let t = T[2] + z * T[3];
+    let w = z * z;
+    let s = z * x;
+    let u = T[0] + z * T[1];
+    r = (x + s * u) + (s * w) * (t + w * r);
+    (if odd { -1. / r } else { r }) as f32
+}

src/math/log.rs

+/* origin: FreeBSD /usr/src/lib/msun/src/e_log.c */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+/* log(x)
+ * Return the logarithm of x
+ *
+ * Method :
+ *   1. Argument Reduction: find k and f such that
+ *                      x = 2^k * (1+f),
+ *         where  sqrt(2)/2 < 1+f < sqrt(2) .
+ *
+ *   2. Approximation of log(1+f).
+ *      Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s)
+ *               = 2s + 2/3 s**3 + 2/5 s**5 + .....,
+ *               = 2s + s*R
+ *      We use a special Remez algorithm on [0,0.1716] to generate
+ *      a polynomial of degree 14 to approximate R The maximum error
+ *      of this polynomial approximation is bounded by 2**-58.45. In
+ *      other words,
+ *                      2      4      6      8      10      12      14
+ *          R(z) ~ Lg1*s +Lg2*s +Lg3*s +Lg4*s +Lg5*s  +Lg6*s  +Lg7*s
+ *      (the values of Lg1 to Lg7 are listed in the program)
+ *      and
+ *          |      2          14          |     -58.45
+ *          | Lg1*s +...+Lg7*s    -  R(z) | <= 2
+ *          |                             |
+ *      Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
+ *      In order to guarantee error in log below 1ulp, we compute log
+ *      by
+ *              log(1+f) = f - s*(f - R)        (if f is not too large)
+ *              log(1+f) = f - (hfsq - s*(hfsq+R)).     (better accuracy)
+ *
+ *      3. Finally,  log(x) = k*ln2 + log(1+f).
+ *                          = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
+ *         Here ln2 is split into two floating point number:
+ *                      ln2_hi + ln2_lo,
+ *         where n*ln2_hi is always exact for |n| < 2000.
+ *
+ * Special cases:
+ *      log(x) is NaN with signal if x < 0 (including -INF) ;
+ *      log(+INF) is +INF; log(0) is -INF with signal;
+ *      log(NaN) is that NaN with no signal.
+ *
+ * Accuracy:
+ *      according to an error analysis, the error is always less than
+ *      1 ulp (unit in the last place).
+ *
+ * Constants:
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
+ * compiler will convert from decimal to binary accurately enough
+ * to produce the hexadecimal values shown.
+ */
+
+const LN2_HI: f64 = 6.93147180369123816490e-01; /* 3fe62e42 fee00000 */
+const LN2_LO: f64 = 1.90821492927058770002e-10; /* 3dea39ef 35793c76 */
+const LG1: f64 = 6.666666666666735130e-01; /* 3FE55555 55555593 */
+const LG2: f64 = 3.999999999940941908e-01; /* 3FD99999 9997FA04 */
+const LG3: f64 = 2.857142874366239149e-01; /* 3FD24924 94229359 */
+const LG4: f64 = 2.222219843214978396e-01; /* 3FCC71C5 1D8E78AF */
+const LG5: f64 = 1.818357216161805012e-01; /* 3FC74664 96CB03DE */
+const LG6: f64 = 1.531383769920937332e-01; /* 3FC39A09 D078C69F */
+const LG7: f64 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
+
+#[inline]
+pub fn log(mut x: f64) -> f64 {
+    let x1p54 = f64::from_bits(0x4350000000000000); // 0x1p54 === 2 ^ 54
+
+    let mut ui = x.to_bits();
+    let mut hx: u32 = (ui >> 32) as u32;
+    let mut k: i32 = 0;
+
+    if (hx < 0x00100000) || ((hx >> 31) != 0) {
+        /* x < 2**-126  */
+        if ui << 1 == 0 {
+            return -1. / (x * x); /* log(+-0)=-inf */
+        }
+        if hx >> 31 != 0 {
+            return (x - x) / 0.0; /* log(-#) = NaN */
+        }
+        /* subnormal number, scale x up */
+        k -= 54;
+        x *= x1p54;
+        ui = x.to_bits();
+        hx = (ui >> 32) as u32;
+    } else if hx >= 0x7ff00000 {
+        return x;
+    } else if hx == 0x3ff00000 && ui << 32 == 0 {
+        return 0.;
+    }
+
+    /* reduce x into [sqrt(2)/2, sqrt(2)] */
+    hx += 0x3ff00000 - 0x3fe6a09e;
+    k += ((hx >> 20) as i32) - 0x3ff;
+    hx = (hx & 0x000fffff) + 0x3fe6a09e;
+    ui = ((hx as u64) << 32) | (ui & 0xffffffff);
+    x = f64::from_bits(ui);
+
+    let f: f64 = x - 1.0;
+    let hfsq: f64 = 0.5 * f * f;
+    let s: f64 = f / (2.0 + f);
+    let z: f64 = s * s;
+    let w: f64 = z * z;
+    let t1: f64 = w * (LG2 + w * (LG4 + w * LG6));
+    let t2: f64 = z * (LG1 + w * (LG3 + w * (LG5 + w * LG7)));
+    let r: f64 = t2 + t1;
+    let dk: f64 = k as f64;
+    return s * (hfsq + r) + dk * LN2_LO - hfsq + f + dk * LN2_HI;
+}

src/math/log10.rs

+/* origin: FreeBSD /usr/src/lib/msun/src/e_log10.c */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+/*
+ * Return the base 10 logarithm of x.  See log.c for most comments.
+ *
+ * Reduce x to 2^k (1+f) and calculate r = log(1+f) - f + f*f/2
+ * as in log.c, then combine and scale in extra precision:
+ *    log10(x) = (f - f*f/2 + r)/log(10) + k*log10(2)
+ */
+
+use core::f64;
+
+const IVLN10HI: f64 = 4.34294481878168880939e-01; /* 0x3fdbcb7b, 0x15200000 */
+const IVLN10LO: f64 = 2.50829467116452752298e-11; /* 0x3dbb9438, 0xca9aadd5 */
+const LOG10_2HI: f64 = 3.01029995663611771306e-01; /* 0x3FD34413, 0x509F6000 */
+const LOG10_2LO: f64 = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */
+const LG1: f64 = 6.666666666666735130e-01; /* 3FE55555 55555593 */
+const LG2: f64 = 3.999999999940941908e-01; /* 3FD99999 9997FA04 */
+const LG3: f64 = 2.857142874366239149e-01; /* 3FD24924 94229359 */
+const LG4: f64 = 2.222219843214978396e-01; /* 3FCC71C5 1D8E78AF */
+const LG5: f64 = 1.818357216161805012e-01; /* 3FC74664 96CB03DE */
+const LG6: f64 = 1.531383769920937332e-01; /* 3FC39A09 D078C69F */
+const LG7: f64 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
+
+#[inline]
+pub fn log10(mut x: f64) -> f64 {
+    let x1p54 = f64::from_bits(0x4350000000000000); // 0x1p54 === 2 ^ 54
+
+    let mut ui: u64 = x.to_bits();
+    let hfsq: f64;
+    let f: f64;
+    let s: f64;
+    let z: f64;
+    let r: f64;
+    let mut w: f64;
+    let t1: f64;
+    let t2: f64;
+    let dk: f64;
+    let y: f64;
+    let mut hi: f64;
+    let lo: f64;
+    let mut val_hi: f64;
+    let mut val_lo: f64;
+    let mut hx: u32;
+    let mut k: i32;
+
+    hx = (ui >> 32) as u32;
+    k = 0;
+    if hx < 0x00100000 || (hx >> 31) > 0 {
+        if ui << 1 == 0 {
+            return -1. / (x * x); /* log(+-0)=-inf */
+        }
+        if (hx >> 31) > 0 {
+            return (x - x) / 0.0; /* log(-#) = NaN */
+        }
+        /* subnormal number, scale x up */
+        k -= 54;
+        x *= x1p54;
+        ui = x.to_bits();
+        hx = (ui >> 32) as u32;
+    } else if hx >= 0x7ff00000 {
+        return x;
+    } else if hx == 0x3ff00000 && ui << 32 == 0 {
+        return 0.;
+    }
+
+    /* reduce x into [sqrt(2)/2, sqrt(2)] */
+    hx += 0x3ff00000 - 0x3fe6a09e;
+    k += (hx >> 20) as i32 - 0x3ff;
+    hx = (hx & 0x000fffff) + 0x3fe6a09e;
+    ui = (hx as u64) << 32 | (ui & 0xffffffff);
+    x = f64::from_bits(ui);
+
+    f = x - 1.0;
+    hfsq = 0.5 * f * f;
+    s = f / (2.0 + f);
+    z = s * s;
+    w = z * z;
+    t1 = w * (LG2 + w * (LG4 + w * LG6));
+    t2 = z * (LG1 + w * (LG3 + w * (LG5 + w * LG7)));
+    r = t2 + t1;
+
+    /* See log2.c for details. */
+    /* hi+lo = f - hfsq + s*(hfsq+R) ~ log(1+f) */
+    hi = f - hfsq;
+    ui = hi.to_bits();
+    ui &= (-1i64 as u64) << 32;
+    hi = f64::from_bits(ui);
+    lo = f - hi - hfsq + s * (hfsq + r);
+
+    /* val_hi+val_lo ~ log10(1+f) + k*log10(2) */
+    val_hi = hi * IVLN10HI;
+    dk = k as f64;
+    y = dk * LOG10_2HI;
+    val_lo = dk * LOG10_2LO + (lo + hi) * IVLN10LO + lo * IVLN10HI;
+
+    /*
+     * Extra precision in for adding y is not strictly needed
+     * since there is no very large cancellation near x = sqrt(2) or
+     * x = 1/sqrt(2), but we do it anyway since it costs little on CPUs
+     * with some parallelism and it reduces the error for many args.
+     */
+    w = y + val_hi;
+    val_lo += (y - w) + val_hi;
+    val_hi = w;
+
+    return val_lo + val_hi;
+}

src/math/log10f.rs

+/* origin: FreeBSD /usr/src/lib/msun/src/e_log10f.c */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+/*
+ * See comments in log10.c.
+ */
+
+use core::f32;
+
+const IVLN10HI: f32 = 4.3432617188e-01; /* 0x3ede6000 */
+const IVLN10LO: f32 = -3.1689971365e-05; /* 0xb804ead9 */
+const LOG10_2HI: f32 = 3.0102920532e-01; /* 0x3e9a2080 */
+const LOG10_2LO: f32 = 7.9034151668e-07; /* 0x355427db */
+/* |(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 */
+
+#[inline]
+pub fn log10f(mut x: f32) -> f32 {
+    let x1p25f = f32::from_bits(0x4c000000); // 0x1p25f === 2 ^ 25
+
+    let mut ui: u32 = x.to_bits();
+    let hfsq: f32;
+    let f: f32;
+    let s: f32;
+    let z: f32;
+    let r: f32;
+    let w: f32;
+    let t1: f32;
+    let t2: f32;
+    let dk: f32;
+    let mut hi: f32;
+    let lo: f32;
+    let mut ix: u32;
+    let mut k: i32;
+
+    ix = ui;
+    k = 0;
+    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.0; /* log(-#) = NaN */
+        }
+        /* subnormal number, scale up x */
+        k -= 25;
+        x *= x1p25f;
+        ui = x.to_bits();
+        ix = ui;
+    } else if ix >= 0x7f800000 {
+        return x;
+    } else if ix == 0x3f800000 {
+        return 0.;
+    }
+
+    /* reduce x into [sqrt(2)/2, sqrt(2)] */
+    ix += 0x3f800000 - 0x3f3504f3;
+    k += (ix >> 23) as i32 - 0x7f;
+    ix = (ix & 0x007fffff) + 0x3f3504f3;
+    ui = ix;
+    x = f32::from_bits(ui);
+
+    f = x - 1.0;
+    s = f / (2.0 + f);
+    z = s * s;
+    w = z * z;
+    t1 = w * (LG2 + w * LG4);
+    t2 = z * (LG1 + w * LG3);
+    r = t2 + t1;
+    hfsq = 0.5 * f * f;
+
+    hi = f - hfsq;
+    ui = hi.to_bits();
+    ui &= 0xfffff000;
+    hi = f32::from_bits(ui);
+    lo = f - hi - hfsq + s * (hfsq + r);
+    dk = k as f32;
+    return dk * LOG10_2LO + (lo + hi) * IVLN10LO + lo * IVLN10HI + hi * IVLN10HI + dk * LOG10_2HI;
+}

src/math/log1p.rs

+/* origin: FreeBSD /usr/src/lib/msun/src/s_log1p.c */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+/* double log1p(double x)
+ * Return the natural logarithm of 1+x.
+ *
+ * Method :
+ *   1. Argument Reduction: find k and f such that
+ *                      1+x = 2^k * (1+f),
+ *         where  sqrt(2)/2 < 1+f < sqrt(2) .
+ *
+ *      Note. If k=0, then f=x is exact. However, if k!=0, then f
+ *      may not be representable exactly. In that case, a correction
+ *      term is need. Let u=1+x rounded. Let c = (1+x)-u, then
+ *      log(1+x) - log(u) ~ c/u. Thus, we proceed to compute log(u),
+ *      and add back the correction term c/u.
+ *      (Note: when x > 2**53, one can simply return log(x))
+ *
+ *   2. Approximation of log(1+f): See log.c
+ *
+ *   3. Finally, log1p(x) = k*ln2 + log(1+f) + c/u. See log.c
+ *
+ * Special cases:
+ *      log1p(x) is NaN with signal if x < -1 (including -INF) ;
+ *      log1p(+INF) is +INF; log1p(-1) is -INF with signal;
+ *      log1p(NaN) is that NaN with no signal.
+ *
+ * Accuracy:
+ *      according to an error analysis, the error is always less than
+ *      1 ulp (unit in the last place).
+ *
+ * Constants:
+ * The hexadecimal values are the intended ones for the following
+ * constants. The decimal values may be used, provided that the
+ * compiler will convert from decimal to binary accurately enough
+ * to produce the hexadecimal values shown.
+ *
+ * Note: Assuming log() return accurate answer, the following
+ *       algorithm can be used to compute log1p(x) to within a few ULP:
+ *
+ *              u = 1+x;
+ *              if(u==1.0) return x ; else
+ *                         return log(u)*(x/(u-1.0));
+ *
+ *       See HP-15C Advanced Functions Handbook, p.193.
+ */
+
+use core::f64;
+
+const LN2_HI: f64 = 6.93147180369123816490e-01; /* 3fe62e42 fee00000 */
+const LN2_LO: f64 = 1.90821492927058770002e-10; /* 3dea39ef 35793c76 */
+const LG1: f64 = 6.666666666666735130e-01; /* 3FE55555 55555593 */
+const LG2: f64 = 3.999999999940941908e-01; /* 3FD99999 9997FA04 */
+const LG3: f64 = 2.857142874366239149e-01; /* 3FD24924 94229359 */
+const LG4: f64 = 2.222219843214978396e-01; /* 3FCC71C5 1D8E78AF */
+const LG5: f64 = 1.818357216161805012e-01; /* 3FC74664 96CB03DE */
+const LG6: f64 = 1.531383769920937332e-01; /* 3FC39A09 D078C69F */
+const LG7: f64 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
+
+#[inline]
+pub fn log1p(x: f64) -> f64 {
+    let mut ui: u64 = x.to_bits();
+    let hfsq: f64;
+    let mut f: f64 = 0.;
+    let mut c: f64 = 0.;
+    let s: f64;
+    let z: f64;
+    let r: f64;
+    let w: f64;
+    let t1: f64;
+    let t2: f64;
+    let dk: f64;
+    let hx: u32;
+    let mut hu: u32;
+    let mut k: i32;
+
+    hx = (ui >> 32) as u32;
+    k = 1;
+    if hx < 0x3fda827a || (hx >> 31) > 0 {
+        /* 1+x < sqrt(2)+ */
+        if hx >= 0xbff00000 {
+            /* x <= -1.0 */
+            if x == -1. {
+                return x / 0.0; /* log1p(-1) = -inf */
+            }
+            return (x - x) / 0.0; /* log1p(x<-1) = NaN */
+        }
+        if hx << 1 < 0x3ca00000 << 1 {
+            /* |x| < 2**-53 */
+            /* underflow if subnormal */
+            if (hx & 0x7ff00000) == 0 {
+                force_eval!(x as f32);
+            }
+            return x;
+        }
+        if hx <= 0xbfd2bec4 {
+            /* sqrt(2)/2- <= 1+x < sqrt(2)+ */
+            k = 0;
+            c = 0.;
+            f = x;
+        }
+    } else if hx >= 0x7ff00000 {
+        return x;
+    }
+    if k > 0 {
+        ui = (1. + x).to_bits();
+        hu = (ui >> 32) as u32;
+        hu += 0x3ff00000 - 0x3fe6a09e;
+        k = (hu >> 20) as i32 - 0x3ff;
+        /* correction term ~ log(1+x)-log(u), avoid underflow in c/u */
+        if k < 54 {
+            c = if k >= 2 {
+                1. - (f64::from_bits(ui) - x)
+            } else {
+                x - (f64::from_bits(ui) - 1.)
+            };
+            c /= f64::from_bits(ui);
+        } else {
+            c = 0.;
+        }
+        /* reduce u into [sqrt(2)/2, sqrt(2)] */
+        hu = (hu & 0x000fffff) + 0x3fe6a09e;
+        ui = (hu as u64) << 32 | (ui & 0xffffffff);
+        f = f64::from_bits(ui) - 1.;
+    }
+    hfsq = 0.5 * f * f;
+    s = f / (2.0 + f);
+    z = s * s;
+    w = z * z;
+    t1 = w * (LG2 + w * (LG4 + w * LG6));
+    t2 = z * (LG1 + w * (LG3 + w * (LG5 + w * LG7)));
+    r = t2 + t1;
+    dk = k as f64;
+    return s * (hfsq + r) + (dk * LN2_LO + c) - hfsq + f + dk * LN2_HI;
+}

src/math/log1pf.rs

+/* origin: FreeBSD /usr/src/lib/msun/src/s_log1pf.c */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
+use core::f32;
+
+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 */
+
+#[inline]
+pub fn log1pf(x: f32) -> f32 {
+    let mut ui: u32 = x.to_bits();
+    let hfsq: f32;
+    let mut f: f32 = 0.;
+    let mut c: f32 = 0.;
+    let s: f32;
+    let z: f32;
+    let r: f32;
+    let w: f32;
+    let t1: f32;
+    let t2: f32;
+    let dk: f32;
+    let ix: u32;
+    let mut iu: u32;
+    let mut k: i32;
+
+    ix = ui;
+    k = 1;
+    if ix < 0x3ed413d0 || (ix >> 31) > 0 {
+        /* 1+x < sqrt(2)+  */
+        if ix >= 0xbf800000 {
+            /* x <= -1.0 */
+            if x == -1. {
+                return x / 0.0; /* log1p(-1)=+inf */
+            }
+            return (x - x) / 0.0; /* log1p(x<-1)=NaN */
+        }
+        if ix << 1 < 0x33800000 << 1 {
+            /* |x| < 2**-24 */
+            /* underflow if subnormal */
+            if (ix & 0x7f800000) == 0 {
+                force_eval!(x * x);
+            }
+            return x;
+        }
+        if ix <= 0xbe95f619 {
+            /* sqrt(2)/2- <= 1+x < sqrt(2)+ */
+            k = 0;
+            c = 0.;
+            f = x;
+        }
+    } else if ix >= 0x7f800000 {
+        return x;
+    }
+    if k > 0 {
+        ui = (1. + x).to_bits();
+        iu = ui;
+        iu += 0x3f800000 - 0x3f3504f3;
+        k = (iu >> 23) as i32 - 0x7f;
+        /* correction term ~ log(1+x)-log(u), avoid underflow in c/u */
+        if k < 25 {
+            c = if k >= 2 {
+                1. - (f32::from_bits(ui) - x)
+            } else {
+                x - (f32::from_bits(ui) - 1.)
+            };
+            c /= f32::from_bits(ui);
+        } else {
+            c = 0.;
+        }
+        /* reduce u into [sqrt(2)/2, sqrt(2)] */
+        iu = (iu & 0x007fffff) + 0x3f3504f3;
+        ui = iu;
+        f = f32::from_bits(ui) - 1.;
+    }
+    s = f / (2.0 + f);
+    z = s * s;
+    w = z * z;
+    t1 = w * (LG2 + w * LG4);
+    t2 = z * (LG1 + w * LG3);
+    r = t2 + t1;
+    hfsq = 0.5 * f * f;
+    dk = k as f32;
+    return s * (hfsq + r) + (dk * LN2_LO + c) - hfsq + f + dk * LN2_HI;
+}

src/math/log2.rs

+/* origin: FreeBSD /usr/src/lib/msun/src/e_log2.c */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+/*
+ * Return the base 2 logarithm of x.  See log.c for most comments.
+ *
+ * Reduce x to 2^k (1+f) and calculate r = log(1+f) - f + f*f/2
+ * as in log.c, then combine and scale in extra precision:
+ *    log2(x) = (f - f*f/2 + r)/log(2) + k
+ */
+
+use core::f64;
+
+const IVLN2HI: f64 = 1.44269504072144627571e+00; /* 0x3ff71547, 0x65200000 */
+const IVLN2LO: f64 = 1.67517131648865118353e-10; /* 0x3de705fc, 0x2eefa200 */
+const LG1: f64 = 6.666666666666735130e-01; /* 3FE55555 55555593 */
+const LG2: f64 = 3.999999999940941908e-01; /* 3FD99999 9997FA04 */
+const LG3: f64 = 2.857142874366239149e-01; /* 3FD24924 94229359 */
+const LG4: f64 = 2.222219843214978396e-01; /* 3FCC71C5 1D8E78AF */
+const LG5: f64 = 1.818357216161805012e-01; /* 3FC74664 96CB03DE */
+const LG6: f64 = 1.531383769920937332e-01; /* 3FC39A09 D078C69F */
+const LG7: f64 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
+
+#[inline]
+pub fn log2(mut x: f64) -> f64 {
+    let x1p54 = f64::from_bits(0x4350000000000000); // 0x1p54 === 2 ^ 54
+
+    let mut ui: u64 = x.to_bits();
+    let hfsq: f64;
+    let f: f64;
+    let s: f64;
+    let z: f64;
+    let r: f64;
+    let mut w: f64;
+    let t1: f64;
+    let t2: f64;
+    let y: f64;
+    let mut hi: f64;
+    let lo: f64;
+    let mut val_hi: f64;
+    let mut val_lo: f64;
+    let mut hx: u32;
+    let mut k: i32;
+
+    hx = (ui >> 32) as u32;
+    k = 0;
+    if hx < 0x00100000 || (hx >> 31) > 0 {
+        if ui << 1 == 0 {
+            return -1. / (x * x); /* log(+-0)=-inf */
+        }
+        if (hx >> 31) > 0 {
+            return (x - x) / 0.0; /* log(-#) = NaN */
+        }
+        /* subnormal number, scale x up */
+        k -= 54;
+        x *= x1p54;
+        ui = x.to_bits();
+        hx = (ui >> 32) as u32;
+    } else if hx >= 0x7ff00000 {
+        return x;
+    } else if hx == 0x3ff00000 && ui << 32 == 0 {
+        return 0.;
+    }
+
+    /* reduce x into [sqrt(2)/2, sqrt(2)] */
+    hx += 0x3ff00000 - 0x3fe6a09e;
+    k += (hx >> 20) as i32 - 0x3ff;
+    hx = (hx & 0x000fffff) + 0x3fe6a09e;
+    ui = (hx as u64) << 32 | (ui & 0xffffffff);
+    x = f64::from_bits(ui);
+
+    f = x - 1.0;
+    hfsq = 0.5 * f * f;
+    s = f / (2.0 + f);
+    z = s * s;
+    w = z * z;
+    t1 = w * (LG2 + w * (LG4 + w * LG6));
+    t2 = z * (LG1 + w * (LG3 + w * (LG5 + w * LG7)));
+    r = t2 + t1;
+
+    /* hi+lo = f - hfsq + s*(hfsq+R) ~ log(1+f) */
+    hi = f - hfsq;
+    ui = hi.to_bits();
+    ui &= (-1i64 as u64) << 32;
+    hi = f64::from_bits(ui);
+    lo = f - hi - hfsq + s * (hfsq + r);
+
+    val_hi = hi * IVLN2HI;
+    val_lo = (lo + hi) * IVLN2LO + lo * IVLN2HI;
+
+    /* spadd(val_hi, val_lo, y), except for not using double_t: */
+    y = k.into();
+    w = y + val_hi;
+    val_lo += (y - w) + val_hi;
+    val_hi = w;
+
+    return val_lo + val_hi;
+}

src/math/log2f.rs

+/* origin: FreeBSD /usr/src/lib/msun/src/e_log2f.c */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+/*
+ * See comments in log2.c.
+ */
+
+use core::f32;
+
+const IVLN2HI: f32 = 1.4428710938e+00; /* 0x3fb8b000 */
+const IVLN2LO: f32 = -1.7605285393e-04; /* 0xb9389ad4 */
+/* |(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 */
+
+#[inline]
+pub fn log2f(mut x: f32) -> f32 {
+    let x1p25f = f32::from_bits(0x4c000000); // 0x1p25f === 2 ^ 25
+
+    let mut ui: u32 = x.to_bits();
+    let hfsq: f32;
+    let f: f32;
+    let s: f32;
+    let z: f32;
+    let r: f32;
+    let w: f32;
+    let t1: f32;
+    let t2: f32;
+    let mut hi: f32;
+    let lo: f32;
+    let mut ix: u32;
+    let mut k: i32;
+
+    ix = ui;
+    k = 0;
+    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.0; /* log(-#) = NaN */
+        }
+        /* subnormal number, scale up x */
+        k -= 25;
+        x *= x1p25f;
+        ui = x.to_bits();
+        ix = ui;
+    } else if ix >= 0x7f800000 {
+        return x;
+    } else if ix == 0x3f800000 {
+        return 0.;
+    }
+
+    /* reduce x into [sqrt(2)/2, sqrt(2)] */
+    ix += 0x3f800000 - 0x3f3504f3;
+    k += (ix >> 23) as i32 - 0x7f;
+    ix = (ix & 0x007fffff) + 0x3f3504f3;
+    ui = ix;
+    x = f32::from_bits(ui);
+
+    f = x - 1.0;
+    s = f / (2.0 + f);
+    z = s * s;
+    w = z * z;
+    t1 = w * (LG2 + w * LG4);
+    t2 = z * (LG1 + w * LG3);
+    r = t2 + t1;
+    hfsq = 0.5 * f * f;
+
+    hi = f - hfsq;
+    ui = hi.to_bits();
+    ui &= 0xfffff000;
+    hi = f32::from_bits(ui);
+    lo = f - hi - hfsq + s * (hfsq + r);
+    return (lo + hi) * IVLN2LO + lo * IVLN2HI + hi * IVLN2HI + k as f32;
+}

src/math/rem_pio2_large.rs

+use super::floor;
+use super::scalbn;
+
+// initial value for jk
+const INIT_JK: [usize; 4] = [3, 4, 4, 6];
+
+// Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
+//
+//              integer array, contains the (24*i)-th to (24*i+23)-th
+//              bit of 2/pi after binary point. The corresponding
+//              floating value is
+//
+//                      ipio2[i] * 2^(-24(i+1)).
+//
+// NB: This table must have at least (e0-3)/24 + jk terms.
+//     For quad precision (e0 <= 16360, jk = 6), this is 686.
+#[cfg(target_pointer_width = "32")]
+const IPIO2: [i32; 66] = [
+    0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 0x95993C, 0x439041, 0xFE5163,
+    0xABDEBB, 0xC561B7, 0x246E3A, 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
+    0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, 0x3991D6, 0x398353, 0x39F49C,
+    0x845F8B, 0xBDF928, 0x3B1FF8, 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
+    0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, 0xF17B3D, 0x0739F7, 0x8A5292,
+    0xEA6BFB, 0x5FB11F, 0x8D5D08, 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
+    0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 0x4D7327, 0x310606, 0x1556CA,
+    0x73A8C9, 0x60E27B, 0xC08C6B,
+];
+
+#[cfg(target_pointer_width = "64")]
+const IPIO2: [i32; 690] = [
+    0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 0x95993C, 0x439041, 0xFE5163,
+    0xABDEBB, 0xC561B7, 0x246E3A, 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129,
+    0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C, 0x7026B4, 0x5F7E41, 0x3991D6, 0x398353, 0x39F49C,
+    0x845F8B, 0xBDF928, 0x3B1FF8, 0x97FFDE, 0x05980F, 0xEF2F11, 0x8B5A0A, 0x6D1F6D, 0x367ECF,
+    0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, 0xF17B3D, 0x0739F7, 0x8A5292,
+    0xEA6BFB, 0x5FB11F, 0x8D5D08, 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3,
+    0x91615E, 0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 0x4D7327, 0x310606, 0x1556CA,
+    0x73A8C9, 0x60E27B, 0xC08C6B, 0x47C419, 0xC367CD, 0xDCE809, 0x2A8359, 0xC4768B, 0x961CA6,
+    0xDDAF44, 0xD15719, 0x053EA5, 0xFF0705, 0x3F7E33, 0xE832C2, 0xDE4F98, 0x327DBB, 0xC33D26,
+    0xEF6B1E, 0x5EF89F, 0x3A1F35, 0xCAF27F, 0x1D87F1, 0x21907C, 0x7C246A, 0xFA6ED5, 0x772D30,
+    0x433B15, 0xC614B5, 0x9D19C3, 0xC2C4AD, 0x414D2C, 0x5D000C, 0x467D86, 0x2D71E3, 0x9AC69B,
+    0x006233, 0x7CD2B4, 0x97A7B4, 0xD55537, 0xF63ED7, 0x1810A3, 0xFC764D, 0x2A9D64, 0xABD770,
+    0xF87C63, 0x57B07A, 0xE71517, 0x5649C0, 0xD9D63B, 0x3884A7, 0xCB2324, 0x778AD6, 0x23545A,
+    0xB91F00, 0x1B0AF1, 0xDFCE19, 0xFF319F, 0x6A1E66, 0x615799, 0x47FBAC, 0xD87F7E, 0xB76522,
+    0x89E832, 0x60BFE6, 0xCDC4EF, 0x09366C, 0xD43F5D, 0xD7DE16, 0xDE3B58, 0x929BDE, 0x2822D2,
+    0xE88628, 0x4D58E2, 0x32CAC6, 0x16E308, 0xCB7DE0, 0x50C017, 0xA71DF3, 0x5BE018, 0x34132E,
+    0x621283, 0x014883, 0x5B8EF5, 0x7FB0AD, 0xF2E91E, 0x434A48, 0xD36710, 0xD8DDAA, 0x425FAE,
+    0xCE616A, 0xA4280A, 0xB499D3, 0xF2A606, 0x7F775C, 0x83C2A3, 0x883C61, 0x78738A, 0x5A8CAF,
+    0xBDD76F, 0x63A62D, 0xCBBFF4, 0xEF818D, 0x67C126, 0x45CA55, 0x36D9CA, 0xD2A828, 0x8D61C2,
+    0x77C912, 0x142604, 0x9B4612, 0xC459C4, 0x44C5C8, 0x91B24D, 0xF31700, 0xAD43D4, 0xE54929,
+    0x10D5FD, 0xFCBE00, 0xCC941E, 0xEECE70, 0xF53E13, 0x80F1EC, 0xC3E7B3, 0x28F8C7, 0x940593,
+    0x3E71C1, 0xB3092E, 0xF3450B, 0x9C1288, 0x7B20AB, 0x9FB52E, 0xC29247, 0x2F327B, 0x6D550C,
+    0x90A772, 0x1FE76B, 0x96CB31, 0x4A1679, 0xE27941, 0x89DFF4, 0x9794E8, 0x84E6E2, 0x973199,
+    0x6BED88, 0x365F5F, 0x0EFDBB, 0xB49A48, 0x6CA467, 0x427271, 0x325D8D, 0xB8159F, 0x09E5BC,
+    0x25318D, 0x3974F7, 0x1C0530, 0x010C0D, 0x68084B, 0x58EE2C, 0x90AA47, 0x02E774, 0x24D6BD,
+    0xA67DF7, 0x72486E, 0xEF169F, 0xA6948E, 0xF691B4, 0x5153D1, 0xF20ACF, 0x339820, 0x7E4BF5,
+    0x6863B2, 0x5F3EDD, 0x035D40, 0x7F8985, 0x295255, 0xC06437, 0x10D86D, 0x324832, 0x754C5B,
+    0xD4714E, 0x6E5445, 0xC1090B, 0x69F52A, 0xD56614, 0x9D0727, 0x50045D, 0xDB3BB4, 0xC576EA,
+    0x17F987, 0x7D6B49, 0xBA271D, 0x296996, 0xACCCC6, 0x5414AD, 0x6AE290, 0x89D988, 0x50722C,
+    0xBEA404, 0x940777, 0x7030F3, 0x27FC00, 0xA871EA, 0x49C266, 0x3DE064, 0x83DD97, 0x973FA3,
+    0xFD9443, 0x8C860D, 0xDE4131, 0x9D3992, 0x8C70DD, 0xE7B717, 0x3BDF08, 0x2B3715, 0xA0805C,
+    0x93805A, 0x921110, 0xD8E80F, 0xAF806C, 0x4BFFDB, 0x0F9038, 0x761859, 0x15A562, 0xBBCB61,
+    0xB989C7, 0xBD4010, 0x04F2D2, 0x277549, 0xF6B6EB, 0xBB22DB, 0xAA140A, 0x2F2689, 0x768364,
+    0x333B09, 0x1A940E, 0xAA3A51, 0xC2A31D, 0xAEEDAF, 0x12265C, 0x4DC26D, 0x9C7A2D, 0x9756C0,
+    0x833F03, 0xF6F009, 0x8C402B, 0x99316D, 0x07B439, 0x15200C, 0x5BC3D8, 0xC492F5, 0x4BADC6,
+    0xA5CA4E, 0xCD37A7, 0x36A9E6, 0x9492AB, 0x6842DD, 0xDE6319, 0xEF8C76, 0x528B68, 0x37DBFC,
+    0xABA1AE, 0x3115DF, 0xA1AE00, 0xDAFB0C, 0x664D64, 0xB705ED, 0x306529, 0xBF5657, 0x3AFF47,
+    0xB9F96A, 0xF3BE75, 0xDF9328, 0x3080AB, 0xF68C66, 0x15CB04, 0x0622FA, 0x1DE4D9, 0xA4B33D,
+    0x8F1B57, 0x09CD36, 0xE9424E, 0xA4BE13, 0xB52333, 0x1AAAF0, 0xA8654F, 0xA5C1D2, 0x0F3F0B,
+    0xCD785B, 0x76F923, 0x048B7B, 0x721789, 0x53A6C6, 0xE26E6F, 0x00EBEF, 0x584A9B, 0xB7DAC4,
+    0xBA66AA, 0xCFCF76, 0x1D02D1, 0x2DF1B1, 0xC1998C, 0x77ADC3, 0xDA4886, 0xA05DF7, 0xF480C6,
+    0x2FF0AC, 0x9AECDD, 0xBC5C3F, 0x6DDED0, 0x1FC790, 0xB6DB2A, 0x3A25A3, 0x9AAF00, 0x9353AD,
+    0x0457B6, 0xB42D29, 0x7E804B, 0xA707DA, 0x0EAA76, 0xA1597B, 0x2A1216, 0x2DB7DC, 0xFDE5FA,
+    0xFEDB89, 0xFDBE89, 0x6C76E4, 0xFCA906, 0x70803E, 0x156E85, 0xFF87FD, 0x073E28, 0x336761,
+    0x86182A, 0xEABD4D, 0xAFE7B3, 0x6E6D8F, 0x396795, 0x5BBF31, 0x48D784, 0x16DF30, 0x432DC7,
+    0x356125, 0xCE70C9, 0xB8CB30, 0xFD6CBF, 0xA200A4, 0xE46C05, 0xA0DD5A, 0x476F21, 0xD21262,
+    0x845CB9, 0x496170, 0xE0566B, 0x015299, 0x375550, 0xB7D51E, 0xC4F133, 0x5F6E13, 0xE4305D,
+    0xA92E85, 0xC3B21D, 0x3632A1, 0xA4B708, 0xD4B1EA, 0x21F716, 0xE4698F, 0x77FF27, 0x80030C,
+    0x2D408D, 0xA0CD4F, 0x99A520, 0xD3A2B3, 0x0A5D2F, 0x42F9B4, 0xCBDA11, 0xD0BE7D, 0xC1DB9B,
+    0xBD17AB, 0x81A2CA, 0x5C6A08, 0x17552E, 0x550027, 0xF0147F, 0x8607E1, 0x640B14, 0x8D4196,
+    0xDEBE87, 0x2AFDDA, 0xB6256B, 0x34897B, 0xFEF305, 0x9EBFB9, 0x4F6A68, 0xA82A4A, 0x5AC44F,
+    0xBCF82D, 0x985AD7, 0x95C7F4, 0x8D4D0D, 0xA63A20, 0x5F57A4, 0xB13F14, 0x953880, 0x0120CC,
+    0x86DD71, 0xB6DEC9, 0xF560BF, 0x11654D, 0x6B0701, 0xACB08C, 0xD0C0B2, 0x485551, 0x0EFB1E,
+    0xC37295, 0x3B06A3, 0x3540C0, 0x7BDC06, 0xCC45E0, 0xFA294E, 0xC8CAD6, 0x41F3E8, 0xDE647C,
+    0xD8649B, 0x31BED9, 0xC397A4, 0xD45877, 0xC5E369, 0x13DAF0, 0x3C3ABA, 0x461846, 0x5F7555,
+    0xF5BDD2, 0xC6926E, 0x5D2EAC, 0xED440E, 0x423E1C, 0x87C461, 0xE9FD29, 0xF3D6E7, 0xCA7C22,
+    0x35916F, 0xC5E008, 0x8DD7FF, 0xE26A6E, 0xC6FDB0, 0xC10893, 0x745D7C, 0xB2AD6B, 0x9D6ECD,
+    0x7B723E, 0x6A11C6, 0xA9CFF7, 0xDF7329, 0xBAC9B5, 0x5100B7, 0x0DB2E2, 0x24BA74, 0x607DE5,
+    0x8AD874, 0x2C150D, 0x0C1881, 0x94667E, 0x162901, 0x767A9F, 0xBEFDFD, 0xEF4556, 0x367ED9,
+    0x13D9EC, 0xB9BA8B, 0xFC97C4, 0x27A831, 0xC36EF1, 0x36C594, 0x56A8D8, 0xB5A8B4, 0x0ECCCF,
+    0x2D8912, 0x34576F, 0x89562C, 0xE3CE99, 0xB920D6, 0xAA5E6B, 0x9C2A3E, 0xCC5F11, 0x4A0BFD,
+    0xFBF4E1, 0x6D3B8E, 0x2C86E2, 0x84D4E9, 0xA9B4FC, 0xD1EEEF, 0xC9352E, 0x61392F, 0x442138,
+    0xC8D91B, 0x0AFC81, 0x6A4AFB, 0xD81C2F, 0x84B453, 0x8C994E, 0xCC2254, 0xDC552A, 0xD6C6C0,
+    0x96190B, 0xB8701A, 0x649569, 0x605A26, 0xEE523F, 0x0F117F, 0x11B5F4, 0xF5CBFC, 0x2DBC34,
+    0xEEBC34, 0xCC5DE8, 0x605EDD, 0x9B8E67, 0xEF3392, 0xB817C9, 0x9B5861, 0xBC57E1, 0xC68351,
+    0x103ED8, 0x4871DD, 0xDD1C2D, 0xA118AF, 0x462C21, 0xD7F359, 0x987AD9, 0xC0549E, 0xFA864F,
+    0xFC0656, 0xAE79E5, 0x362289, 0x22AD38, 0xDC9367, 0xAAE855, 0x382682, 0x9BE7CA, 0xA40D51,
+    0xB13399, 0x0ED7A9, 0x480569, 0xF0B265, 0xA7887F, 0x974C88, 0x36D1F9, 0xB39221, 0x4A827B,
+    0x21CF98, 0xDC9F40, 0x5547DC, 0x3A74E1, 0x42EB67, 0xDF9DFE, 0x5FD45E, 0xA4677B, 0x7AACBA,
+    0xA2F655, 0x23882B, 0x55BA41, 0x086E59, 0x862A21, 0x834739, 0xE6E389, 0xD49EE5, 0x40FB49,
+    0xE956FF, 0xCA0F1C, 0x8A59C5, 0x2BFA94, 0xC5C1D3, 0xCFC50F, 0xAE5ADB, 0x86C547, 0x624385,
+    0x3B8621, 0x94792C, 0x876110, 0x7B4C2A, 0x1A2C80, 0x12BF43, 0x902688, 0x893C78, 0xE4C4A8,
+    0x7BDBE5, 0xC23AC4, 0xEAF426, 0x8A67F7, 0xBF920D, 0x2BA365, 0xB1933D, 0x0B7CBD, 0xDC51A4,
+    0x63DD27, 0xDDE169, 0x19949A, 0x9529A8, 0x28CE68, 0xB4ED09, 0x209F44, 0xCA984E, 0x638270,
+    0x237C7E, 0x32B90F, 0x8EF5A7, 0xE75614, 0x08F121, 0x2A9DB5, 0x4D7E6F, 0x5119A5, 0xABF9B5,
+    0xD6DF82, 0x61DD96, 0x023616, 0x9F3AC4, 0xA1A283, 0x6DED72, 0x7A8D39, 0xA9B882, 0x5C326B,
+    0x5B2746, 0xED3400, 0x7700D2, 0x55F4FC, 0x4D5901, 0x8071E0,
+];
+
+const PIO2: [f64; 8] = [
+    1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
+    7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
+    5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */
+    3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */
+    1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */
+    1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */
+    2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */
+    2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
+];
+
+// fn rem_pio2_large(x : &[f64], y : &mut [f64], e0 : i32, prec : usize) -> i32
+//
+// Input parameters:
+//      x[]     The input value (must be positive) is broken into nx
+//              pieces of 24-bit integers in double precision format.
+//              x[i] will be the i-th 24 bit of x. The scaled exponent
+//              of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
+//              match x's up to 24 bits.
+//
+//              Example of breaking a double positive z into x[0]+x[1]+x[2]:
+//                      e0 = ilogb(z)-23
+//                      z  = scalbn(z,-e0)
+//              for i = 0,1,2
+//                      x[i] = floor(z)
+//                      z    = (z-x[i])*2**24
+//
+//      y[]     ouput result in an array of double precision numbers.
+//              The dimension of y[] is:
+//                      24-bit  precision       1
+//                      53-bit  precision       2
+//                      64-bit  precision       2
+//                      113-bit precision       3
+//              The actual value is the sum of them. Thus for 113-bit
+//              precison, one may have to do something like:
+//
+//              long double t,w,r_head, r_tail;
+//              t = (long double)y[2] + (long double)y[1];
+//              w = (long double)y[0];
+//              r_head = t+w;
+//              r_tail = w - (r_head - t);
+//
+//      e0      The exponent of x[0]. Must be <= 16360 or you need to
+//              expand the ipio2 table.
+//
+//      prec    an integer indicating the precision:
+//                      0       24  bits (single)
+//                      1       53  bits (double)
+//                      2       64  bits (extended)
+//                      3       113 bits (quad)
+//
+// Here is the description of some local variables:
+//
+//      jk      jk+1 is the initial number of terms of ipio2[] needed
+//              in the computation. The minimum and recommended value
+//              for jk is 3,4,4,6 for single, double, extended, and quad.
+//              jk+1 must be 2 larger than you might expect so that our
+//              recomputation test works. (Up to 24 bits in the integer
+//              part (the 24 bits of it that we compute) and 23 bits in
+//              the fraction part may be lost to cancelation before we
+//              recompute.)
+//
+//      jz      local integer variable indicating the number of
+//              terms of ipio2[] used.
+//
+//      jx      nx - 1
+//
+//      jv      index for pointing to the suitable ipio2[] for the
+//              computation. In general, we want
+//                      ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8
+//              is an integer. Thus
+//                      e0-3-24*jv >= 0 or (e0-3)/24 >= jv
+//              Hence jv = max(0,(e0-3)/24).
+//
+//      jp      jp+1 is the number of terms in PIo2[] needed, jp = jk.
+//
+//      q[]     double array with integral value, representing the
+//              24-bits chunk of the product of x and 2/pi.
+//
+//      q0      the corresponding exponent of q[0]. Note that the
+//              exponent for q[i] would be q0-24*i.
+//
+//      PIo2[]  double precision array, obtained by cutting pi/2
+//              into 24 bits chunks.
+//
+//      f[]     ipio2[] in floating point
+//
+//      iq[]    integer array by breaking up q[] in 24-bits chunk.
+//
+//      fq[]    final product of x*(2/pi) in fq[0],..,fq[jk]
+//
+//      ih      integer. If >0 it indicates q[] is >= 0.5, hence
+//              it also indicates the *sign* of the result.
+
+/// Return the last three digits of N with y = x - N*pi/2
+/// so that |y| < pi/2.
+///
+/// The method is to compute the integer (mod 8) and fraction parts of
+/// (2/pi)*x without doing the full multiplication. In general we
+/// skip the part of the product that are known to be a huge integer (
+/// more accurately, = 0 mod 8 ). Thus the number of operations are
+/// independent of the exponent of the input.
+#[inline]
+pub(crate) fn rem_pio2_large(x: &[f64], y: &mut [f64], e0: i32, prec: usize) -> i32 {
+    let x1p24 = f64::from_bits(0x4170000000000000); // 0x1p24 === 2 ^ 24
+    let x1p_24 = f64::from_bits(0x3e70000000000000); // 0x1p_24 === 2 ^ (-24)
+
+    #[cfg(target_pointer_width = "64")]
+    assert!(e0 <= 16360);
+
+    let nx = x.len();
+
+    let mut fw: f64;
+    let mut n: i32;
+    let mut ih: i32;
+    let mut z: f64;
+    let mut f: [f64; 20] = [0.; 20];
+    let mut fq: [f64; 20] = [0.; 20];
+    let mut q: [f64; 20] = [0.; 20];
+    let mut iq: [i32; 20] = [0; 20];
+
+    /* initialize jk*/
+    let jk = INIT_JK[prec];
+    let jp = jk;
+
+    /* determine jx,jv,q0, note that 3>q0 */
+    let jx = nx - 1;
+    let mut jv = (e0 - 3) / 24;
+    if jv < 0 {
+        jv = 0;
+    }
+    let mut q0 = e0 - 24 * (jv + 1);
+    let jv = jv as usize;
+
+    /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
+    let mut j = (jv - jx) as i32;
+    let m = jx + jk;
+    for i in 0..=m {
+        f[i] = if j < 0 { 0. } else { IPIO2[j as usize] as f64 };
+        j += 1
+    }
+
+    /* compute q[0],q[1],...q[jk] */
+    for i in 0..=jk {
+        fw = 0f64;
+        for j in 0..=jx {
+            fw += x[j] * f[jx + i - j];
+        }
+        q[i] = fw;
+    }
+
+    let mut jz = jk;
+
+    'recompute: loop {
+        /* distill q[] into iq[] reversingly */
+        let mut i = 0i32;
+        z = q[jz];
+        for j in (1..=jz).rev() {
+            fw = (x1p_24 * z) as i32 as f64;
+            iq[i as usize] = (z - x1p24 * fw) as i32;
+            z = q[j - 1] + fw;
+            i += 1;
+        }
+
+        /* compute n */
+        z = scalbn(z, q0); /* actual value of z */
+        z -= 8.0 * floor(z * 0.125); /* trim off integer >= 8 */
+        n = z as i32;
+        z -= n as f64;
+        ih = 0;
+        if q0 > 0 {
+            /* need iq[jz-1] to determine n */
+            i = iq[jz - 1] >> (24 - q0);
+            n += i;
+            iq[jz - 1] -= i << (24 - q0);
+            ih = iq[jz - 1] >> (23 - q0);
+        } else if q0 == 0 {
+            ih = iq[jz - 1] >> 23;
+        } else if z >= 0.5 {
+            ih = 2;
+        }
+
+        if ih > 0 {
+            /* q > 0.5 */
+            n += 1;
+            let mut carry = 0i32;
+            for i in 0..jz {
+                /* compute 1-q */
+                let j = iq[i];
+                if carry == 0 {
+                    if j != 0 {
+                        carry = 1;
+                        iq[i] = 0x1000000 - j;
+                    }
+                } else {
+                    iq[i] = 0xffffff - j;
+                }
+            }
+            if q0 > 0 {
+                /* rare case: chance is 1 in 12 */
+                match q0 {
+                    1 => {
+                        iq[jz - 1] &= 0x7fffff;
+                    }
+                    2 => {
+                        iq[jz - 1] &= 0x3fffff;
+                    }
+                    _ => {}
+                }
+            }
+            if ih == 2 {
+                z = 1. - z;
+                if carry != 0 {
+                    z -= scalbn(1., q0);
+                }
+            }
+        }
+
+        /* check if recomputation is needed */
+        if z == 0. {
+            let mut j = 0;
+            for i in (jk..=jz - 1).rev() {
+                j |= iq[i];
+            }
+            if j == 0 {
+                /* need recomputation */
+                let mut k = 1;
+                while iq[jk - k] == 0 {
+                    k += 1; /* k = no. of terms needed */
+                }
+
+                for i in (jz + 1)..=(jz + k) {
+                    /* add q[jz+1] to q[jz+k] */
+                    f[jx + i] = IPIO2[jv + i] as f64;
+                    fw = 0f64;
+                    for j in 0..=jx {
+                        fw += x[j] * f[jx + i - j];
+                    }
+                    q[i] = fw;
+                }
+                jz += k;
+                continue 'recompute;
+            }
+        }
+
+        break;
+    }
+
+    /* chop off zero terms */
+    if z == 0. {
+        jz -= 1;
+        q0 -= 24;
+        while iq[jz] == 0 {
+            jz -= 1;
+            q0 -= 24;
+        }
+    } else {
+        /* break z into 24-bit if necessary */
+        z = scalbn(z, -q0);
+        if z >= x1p24 {
+            fw = (x1p_24 * z) as i32 as f64;
+            iq[jz] = (z - x1p24 * fw) as i32;
+            jz += 1;
+            q0 += 24;
+            iq[jz] = fw as i32;
+        } else {
+            iq[jz] = z as i32;
+        }
+    }
+
+    /* convert integer "bit" chunk to floating-point value */
+    fw = scalbn(1., q0);
+    for i in (0..=jz).rev() {
+        q[i] = fw * (iq[i] as f64);
+        fw *= x1p_24;
+    }
+
+    /* compute PIo2[0,...,jp]*q[jz,...,0] */
+    for i in (0..=jz).rev() {
+        fw = 0f64;
+        let mut k = 0;
+        while (k <= jp) && (k <= jz - i) {
+            fw += PIO2[k] * q[i + k];
+            k += 1;
+        }
+        fq[jz - i] = fw;
+    }
+
+    /* compress fq[] into y[] */
+    match prec {
+        0 => {
+            fw = 0f64;
+            for i in (0..=jz).rev() {
+                fw += fq[i];
+            }
+            y[0] = if ih == 0 { fw } else { -fw };
+        }
+        1 | 2 => {
+            fw = 0f64;
+            for i in (0..=jz).rev() {
+                fw += fq[i];
+            }
+            // TODO: drop excess precision here once double_t is used
+            fw = fw as f64;
+            y[0] = if ih == 0 { fw } else { -fw };
+            fw = fq[0] - fw;
+            for i in 1..=jz {
+                fw += fq[i];
+            }
+            y[1] = if ih == 0 { fw } else { -fw };
+        }
+        3 => {
+            /* painful */
+            for i in (1..=jz).rev() {
+                fw = fq[i - 1] + fq[i];
+                fq[i] += fq[i - 1] - fw;
+                fq[i - 1] = fw;
+            }
+            for i in (2..=jz).rev() {
+                fw = fq[i - 1] + fq[i];
+                fq[i] += fq[i - 1] - fw;
+                fq[i - 1] = fw;
+            }
+            fw = 0f64;
+            for i in (2..=jz).rev() {
+                fw += fq[i];
+            }
+            if ih == 0 {
+                y[0] = fq[0];
+                y[1] = fq[1];
+                y[2] = fw;
+            } else {
+                y[0] = -fq[0];
+                y[1] = -fq[1];
+                y[2] = -fw;
+            }
+        }
+        _ => unreachable!(),
+    }
+    n & 7
+}

src/math/rem_pio2f.rs

+use super::rem_pio2_large;
+
+use core::f64;
+
+const TOINT: f64 = 1.5 / f64::EPSILON;
+
+/// 53 bits of 2/pi
+const INV_PIO2: f64 = 6.36619772367581382433e-01; /* 0x3FE45F30, 0x6DC9C883 */
+/// first 25 bits of pi/2
+const PIO2_1: f64 = 1.57079631090164184570e+00; /* 0x3FF921FB, 0x50000000 */
+/// pi/2 - pio2_1
+const PIO2_1T: f64 = 1.58932547735281966916e-08; /* 0x3E5110b4, 0x611A6263 */
+
+/// Return the remainder of x rem pi/2 in *y
+///
+/// use double precision for everything except passing x
+/// use __rem_pio2_large() for large x
+#[inline]
+pub(crate) fn rem_pio2f(x: f32) -> (i32, f64) {
+    let x64 = x as f64;
+
+    let mut tx: [f64; 1] = [0.];
+    let mut ty: [f64; 1] = [0.];
+
+    let ix = x.to_bits() & 0x7fffffff;
+    /* 25+53 bit pi is good enough for medium size */
+    if ix < 0x4dc90fdb {
+        /* |x| ~< 2^28*(pi/2), medium size */
+        /* Use a specialized rint() to get fn.  Assume round-to-nearest. */
+        let f_n = x64 * INV_PIO2 + TOINT - TOINT;
+        return (f_n as i32, x64 - f_n * PIO2_1 - f_n * PIO2_1T);
+    }
+    if ix >= 0x7f800000 {
+        /* x is inf or NaN */
+        return (0, x64 - x64);
+    }
+    /* scale x into [2^23, 2^24-1] */
+    let sign = (x.to_bits() >> 31) != 0;
+    let e0 = ((ix >> 23) - (0x7f + 23)) as i32; /* e0 = ilogb(|x|)-23, positive */
+    tx[0] = f32::from_bits(ix - (e0 << 23) as u32) as f64;
+    let n = rem_pio2_large(&tx, &mut ty, e0, 0);
+    if sign {
+        return (-n, -ty[0]);
+    }
+    (n, ty[0])
+}

src/math/roundf.rs

+use core::f32;
+
+const TOINT: f32 = 1.0 / f32::EPSILON;
+
+#[inline]
+pub fn roundf(mut x: f32) -> f32 {
+    let i = x.to_bits();
+    let e: u32 = i >> 23 & 0xff;
+    let mut y: f32;
+
+    if e >= 0x7f + 23 {
+        return x;
+    }
+    if i >> 31 != 0 {
+        x = -x;
+    }
+    if e < 0x7f - 1 {
+        force_eval!(x + TOINT);
+        return 0.0 * x;
+    }
+    y = x + TOINT - TOINT - x;
+    if y > 0.5f32 {
+        y = y + x - 1.0;
+    } else if y <= -0.5f32 {
+        y = y + x + 1.0;
+    } else {
+        y = y + x;
+    }
+    if i >> 31 != 0 {
+        -y
+    } else {
+        y
+    }
+}

src/math/sinf.rs

+use super::{k_cosf, k_sinf, rem_pio2f};
+
+use core::f64::consts::FRAC_PI_2;
+
+/* Small multiples of pi/2 rounded to double precision. */
+const S1_PIO2: f64 = 1. * FRAC_PI_2; /* 0x3FF921FB, 0x54442D18 */
+const S2_PIO2: f64 = 2. * FRAC_PI_2; /* 0x400921FB, 0x54442D18 */
+const S3_PIO2: f64 = 3. * FRAC_PI_2; /* 0x4012D97C, 0x7F3321D2 */
+const S4_PIO2: f64 = 4. * FRAC_PI_2; /* 0x401921FB, 0x54442D18 */
+
+#[inline]
+pub fn sinf(x: f32) -> f32 {
+    let x64 = x as f64;
+
+    let x1p120 = f32::from_bits(0x7b800000); // 0x1p120f === 2 ^ 120
+
+    let mut ix = x.to_bits();
+    let sign = (ix >> 31) != 0;
+    ix &= 0x7fffffff;
+
+    if ix <= 0x3f490fda {
+        /* |x| ~<= pi/4 */
+        if ix < 0x39800000 {
+            /* |x| < 2**-12 */
+            /* raise inexact if x!=0 and underflow if subnormal */
+            force_eval!(if ix < 0x00800000 {
+                x / x1p120
+            } else {
+                x + x1p120
+            });
+            return x;
+        }
+        return k_sinf(x64);
+    }
+    if ix <= 0x407b53d1 {
+        /* |x| ~<= 5*pi/4 */
+        if ix <= 0x4016cbe3 {
+            /* |x| ~<= 3pi/4 */
+            if sign {
+                return -k_cosf(x64 + S1_PIO2);
+            } else {
+                return k_cosf(x64 - S1_PIO2);
+            }
+        }
+        return k_sinf(if sign {
+            -(x64 + S2_PIO2)
+        } else {
+            -(x64 - S2_PIO2)
+        });
+    }
+    if ix <= 0x40e231d5 {
+        /* |x| ~<= 9*pi/4 */
+        if ix <= 0x40afeddf {
+            /* |x| ~<= 7*pi/4 */
+            if sign {
+                return k_cosf(x64 + S3_PIO2);
+            } else {
+                return -k_cosf(x64 - S3_PIO2);
+            }
+        }
+        return k_sinf(if sign { x64 + S4_PIO2 } else { x64 - S4_PIO2 });
+    }
+
+    /* sin(Inf or NaN) is NaN */
+    if ix >= 0x7f800000 {
+        return x - x;
+    }
+
+    /* general argument reduction needed */
+    let (n, y) = rem_pio2f(x);
+    match n & 3 {
+        0 => k_sinf(y),
+        1 => k_cosf(y),
+        2 => return k_sinf(-y),
+        _ => -k_cosf(y),
+    }
+}

src/math/tanf.rs

+use super::{k_tanf, rem_pio2f};
+
+use core::f64::consts::FRAC_PI_2;
+
+/* Small multiples of pi/2 rounded to double precision. */
+const T1_PIO2: f64 = 1. * FRAC_PI_2; /* 0x3FF921FB, 0x54442D18 */
+const T2_PIO2: f64 = 2. * FRAC_PI_2; /* 0x400921FB, 0x54442D18 */
+const T3_PIO2: f64 = 3. * FRAC_PI_2; /* 0x4012D97C, 0x7F3321D2 */
+const T4_PIO2: f64 = 4. * FRAC_PI_2; /* 0x401921FB, 0x54442D18 */
+
+#[inline]
+pub fn tanf(x: f32) -> f32 {
+    let x64 = x as f64;
+
+    let x1p120 = f32::from_bits(0x7b800000); // 0x1p120f === 2 ^ 120
+
+    let mut ix = x.to_bits();
+    let sign = (ix >> 31) != 0;
+    ix &= 0x7fffffff;
+
+    if ix <= 0x3f490fda {
+        /* |x| ~<= pi/4 */
+        if ix < 0x39800000 {
+            /* |x| < 2**-12 */
+            /* raise inexact if x!=0 and underflow if subnormal */
+            force_eval!(if ix < 0x00800000 {
+                x / x1p120
+            } else {
+                x + x1p120
+            });
+            return x;
+        }
+        return k_tanf(x64, false);
+    }
+    if ix <= 0x407b53d1 {
+        /* |x| ~<= 5*pi/4 */
+        if ix <= 0x4016cbe3 {
+            /* |x| ~<= 3pi/4 */
+            return k_tanf(if sign { x64 + T1_PIO2 } else { x64 - T1_PIO2 }, true);
+        } else {
+            return k_tanf(if sign { x64 + T2_PIO2 } else { x64 - T2_PIO2 }, false);
+        }
+    }
+    if ix <= 0x40e231d5 {
+        /* |x| ~<= 9*pi/4 */
+        if ix <= 0x40afeddf {
+            /* |x| ~<= 7*pi/4 */
+            return k_tanf(if sign { x64 + T3_PIO2 } else { x64 - T3_PIO2 }, true);
+        } else {
+            return k_tanf(if sign { x64 + T4_PIO2 } else { x64 - T4_PIO2 }, false);
+        }
+    }
+
+    /* tan(Inf or NaN) is NaN */
+    if ix >= 0x7f800000 {
+        return x - x;
+    }
+
+    /* argument reduction */
+    let (n, y) = rem_pio2f(x);
+    k_tanf(y, n & 1 != 0)
+}

.travis.yml

@@ -18,2 +18,4 @@ language: rust
     - env: TARGET=x86_64-unknown-linux-gnu
+    - env: TARGET=cargo-fmt
+      rust: beta
 
@@ -37,4 +39,3 @@ before_install: set -e
   only:
-    - master
     - staging
     - trying

Cargo.toml

@@ -15,3 +15,3 @@ # THIS FILE IS AUTOMATICALLY GENERATED BY CARGO
 name = "libm"
-version = "0.1.0"
+version = "0.1.1"
 authors = ["Jorge Aparicio <[e-mail]>"]
@@ -18,0 +18,0 @@ description = "libm in pure Rust"

CHANGELOG.md

@@ -8,2 +8,37 @@ # Change Log
 
+## [v0.1.1] - 2018-07-14
+
+### Added
+
+- acos
+- acosf
+- asin
+- asinf
+- atanf
+- cbrt
+- cbrtf
+- ceil
+- ceilf
+- cosf
+- exp
+- exp2
+- exp2f
+- expm1
+- expm1f
+- fdim
+- fdimf
+- floorf
+- fma
+- fmod
+- log
+- log2
+- log10
+- log10f
+- log1p
+- log1pf
+- log2f
+- roundf
+- sinf
+- tanf
+
 ## v0.1.0 - 2018-07-13
@@ -13,2 +48,3 @@
 
-[Unreleased]: https://github.com/japaric/libm/compare/v0.1.0...HEAD
+[Unreleased]: https://github.com/japaric/libm/compare/v0.1.1...HEAD
+[v0.1.1]: https://github.com/japaric/libm/compare/v0.1.0...v0.1.1

ci/install.sh

 set -euxo pipefail
 
 main() {
+    if [ $TARGET = cargo-fmt ]; then
+        rustup component add rustfmt-preview
+        return
+    fi
+
     if ! hash cross >/dev/null 2>&1; then
@@ -5,0 +10,0 @@ cargo install cross

ci/script.sh

 set -euxo pipefail
 
 main() {
+    if [ $TARGET = cargo-fmt ]; then
+        cargo fmt -- --check
+        return
+    fi
+
     # quick check
@@ -5,0 +10,0 @@ cargo check

CONTRIBUTING.md

@@ -13,3 +13,4 @@ # How to contribute
   you can run unit tests with the `cargo test --lib` command.
-- Send us a pull request!
+- Send us a pull request! Make sure to run `cargo fmt` on your code before sending the PR. Also
+  include "closes #42" in the PR description to close the corresponding issue.
 - :tada:
@@ -31,3 +32,3 @@
 - Only use relative imports within the `math` directory / module, e.g. `use self::fabs::fabs` or
-`use super::isnanf`. Absolute imports from core are OK, e.g. `use core::u64`.
+`use super::k_cos`. Absolute imports from core are OK, e.g. `use core::u64`.
 
@@ -40,2 +41,9 @@ - To reinterpret a float as an integer use the `to_bits` method. The MUSL code uses the
 
+- You may use other methods from core like `f64::is_nan`, etc. as appropriate.
+
+- If you're implementing one of the private double-underscore functions, take a look at the
+  "source" name in the comment at the top for an idea for alternate naming. For example, `__sin`
+  was renamed to `k_sin` after the FreeBSD source code naming. Do `use` these private functions in
+  `mod.rs`.
+
 - You may encounter weird literals like `0x1p127f` in the MUSL code. These are hexadecimal floating
@@ -42,0 +50,0 @@ point literals. Rust (the language) doesn't support these kind of literals. The best way I have

README.md

@@ -17,2 +17,26 @@ # `libm`
 
+## Already usable
+
+This crate is [on crates.io] and can be used today in stable `#![no_std]` programs like this:
+
+[on crates.io]: https://crates.io/crates/libm
+
+``` rust
+#![no_std]
+
+extern crate libm;
+
+use libm::F32Ext; // adds methods to `f32`
+
+fn foo(x: f32) {
+    let y = x.sqrt();
+    let z = libm::truncf(x);
+}
+```
+
+Not all the math functions are available at the moment. Check the [API docs] to learn what's
+currently supported.
+
+[API docs]: https://docs.rs/libm
+
 ## Contributing
@@ -19,0 +43,0 @@

src/lib.rs

@@ -37,9 +37,6 @@ //! Port of MUSL's libm to Rust
 pub trait F32Ext: private::Sealed {
-    #[cfg(todo)]
     fn floor(self) -> Self;
 
-    #[cfg(todo)]
     fn ceil(self) -> Self;
 
-    #[cfg(todo)]
     fn round(self) -> Self;
@@ -49,2 +46,4 @@
 
+    fn fdim(self, rhs: Self) -> Self;
+
     #[cfg(todo)]
@@ -76,3 +75,2 @@ fn fract(self) -> Self;
 
-    #[cfg(todo)]
     fn exp2(self) -> Self;
@@ -84,9 +82,6 @@
 
-    #[cfg(todo)]
     fn log2(self) -> Self;
 
-    #[cfg(todo)]
     fn log10(self) -> Self;
 
-    #[cfg(todo)]
     fn cbrt(self) -> Self;
@@ -96,18 +91,12 @@
 
-    #[cfg(todo)]
     fn sin(self) -> Self;
 
-    #[cfg(todo)]
     fn cos(self) -> Self;
 
-    #[cfg(todo)]
     fn tan(self) -> Self;
 
-    #[cfg(todo)]
     fn asin(self) -> Self;
 
-    #[cfg(todo)]
     fn acos(self) -> Self;
 
-    #[cfg(todo)]
     fn atan(self) -> Self;
@@ -124,6 +113,4 @@
 
-    #[cfg(todo)]
     fn exp_m1(self) -> Self;
 
-    #[cfg(todo)]
     fn ln_1p(self) -> Self;
@@ -151,3 +138,2 @@
 impl F32Ext for f32 {
-    #[cfg(todo)]
     #[inline]
@@ -158,3 +144,2 @@ fn floor(self) -> Self {
 
-    #[cfg(todo)]
     #[inline]
@@ -165,3 +150,2 @@ fn ceil(self) -> Self {
 
-    #[cfg(todo)]
     #[inline]
@@ -177,2 +161,7 @@ fn round(self) -> Self {
 
+    #[inline]
+    fn fdim(self, rhs: Self) -> Self {
+        fdimf(self, rhs)
+    }
+
     #[cfg(todo)]
@@ -231,3 +220,2 @@ #[inline]
 
-    #[cfg(todo)]
     #[inline]
@@ -248,3 +236,2 @@ fn exp2(self) -> Self {
 
-    #[cfg(todo)]
     #[inline]
@@ -255,3 +242,2 @@ fn log2(self) -> Self {
 
-    #[cfg(todo)]
     #[inline]
@@ -262,3 +248,2 @@ fn log10(self) -> Self {
 
-    #[cfg(todo)]
     #[inline]
@@ -274,3 +259,2 @@ fn cbrt(self) -> Self {
 
-    #[cfg(todo)]
     #[inline]
@@ -281,3 +265,2 @@ fn sin(self) -> Self {
 
-    #[cfg(todo)]
     #[inline]
@@ -288,3 +271,2 @@ fn cos(self) -> Self {
 
-    #[cfg(todo)]
     #[inline]
@@ -295,3 +277,2 @@ fn tan(self) -> Self {
 
-    #[cfg(todo)]
     #[inline]
@@ -302,3 +283,2 @@ fn asin(self) -> Self {
 
-    #[cfg(todo)]
     #[inline]
@@ -309,3 +289,2 @@ fn acos(self) -> Self {
 
-    #[cfg(todo)]
     #[inline]
@@ -322,3 +301,2 @@ fn atan(self) -> Self {
 
-    #[cfg(todo)]
     #[inline]
@@ -329,3 +307,2 @@ fn exp_m1(self) -> Self {
 
-    #[cfg(todo)]
     #[inline]
@@ -386,3 +363,2 @@ fn ln_1p(self) -> Self {
 
-    #[cfg(todo)]
     fn ceil(self) -> Self;
@@ -394,2 +370,4 @@
 
+    fn fdim(self, rhs: Self) -> Self;
+
     #[cfg(todo)]
@@ -403,3 +381,2 @@ fn fract(self) -> Self;
 
-    #[cfg(todo)]
     fn mul_add(self, a: Self, b: Self) -> Self;
@@ -421,21 +398,14 @@
 
-    #[cfg(todo)]
     fn exp(self) -> Self;
 
-    #[cfg(todo)]
     fn exp2(self) -> Self;
 
-    #[cfg(todo)]
     fn ln(self) -> Self;
 
-    #[cfg(todo)]
     fn log(self, base: Self) -> Self;
 
-    #[cfg(todo)]
     fn log2(self) -> Self;
 
-    #[cfg(todo)]
     fn log10(self) -> Self;
 
-    #[cfg(todo)]
     fn cbrt(self) -> Self;
@@ -457,3 +427,2 @@
 
-    #[cfg(todo)]
     fn acos(self) -> Self;
@@ -473,6 +442,4 @@
 
-    #[cfg(todo)]
     fn exp_m1(self) -> Self;
 
-    #[cfg(todo)]
     fn ln_1p(self) -> Self;
@@ -505,3 +472,2 @@
 
-    #[cfg(todo)]
     #[inline]
@@ -522,2 +488,6 @@ fn ceil(self) -> Self {
 
+    #[inline]
+    fn fdim(self, rhs: Self) -> Self {
+        fdim(self, rhs)
+    }
     #[cfg(todo)]
@@ -534,3 +504,2 @@ #[inline]
 
-    #[cfg(todo)]
     #[inline]
@@ -573,3 +542,2 @@ fn mul_add(self, a: Self, b: Self) -> Self {
 
-    #[cfg(todo)]
     #[inline]
@@ -580,3 +548,2 @@ fn exp(self) -> Self {
 
-    #[cfg(todo)]
     #[inline]
@@ -587,3 +554,2 @@ fn exp2(self) -> Self {
 
-    #[cfg(todo)]
     #[inline]
@@ -594,3 +560,2 @@ fn ln(self) -> Self {
 
-    #[cfg(todo)]
     #[inline]
@@ -601,3 +566,2 @@ fn log(self, base: Self) -> Self {
 
-    #[cfg(todo)]
     #[inline]
@@ -608,3 +572,2 @@ fn log2(self) -> Self {
 
-    #[cfg(todo)]
     #[inline]
@@ -615,3 +578,2 @@ fn log10(self) -> Self {
 
-    #[cfg(todo)]
     #[inline]
@@ -651,3 +613,2 @@ fn cbrt(self) -> Self {
 
-    #[cfg(todo)]
     #[inline]
@@ -670,3 +631,2 @@ fn acos(self) -> Self {
 
-    #[cfg(todo)]
     #[inline]
@@ -677,3 +637,2 @@ fn exp_m1(self) -> Self {
 
-    #[cfg(todo)]
     #[inline]
@@ -680,0 +639,0 @@ fn ln_1p(self) -> Self {

src/math/expf.rs

@@ -0,7 +1,22 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/e_expf.c */
+/*
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, [e-mail].
+ */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
 use super::scalbnf;
 
-const HALF    : [f32; 2] = [0.5,-0.5];
-const LN2_HI  : f32      = 6.9314575195e-01;  /* 0x3f317200 */
-const LN2_LO  : f32      = 1.4286067653e-06;  /* 0x35bfbe8e */
-const INV_LN2 : f32      = 1.4426950216e+00;  /* 0x3fb8aa3b */
+const HALF: [f32; 2] = [0.5, -0.5];
+const LN2_HI: f32 = 6.9314575195e-01; /* 0x3f317200 */
+const LN2_LO: f32 = 1.4286067653e-06; /* 0x35bfbe8e */
+const INV_LN2: f32 = 1.4426950216e+00; /* 0x3fb8aa3b */
 /*
@@ -11,21 +26,24 @@ * Domain [-0.34568, 0.34568], range ~[-4.278e-9, 4.447e-9]:
  */
-const P1 : f32 =  1.6666625440e-1;    /*  0xaaaa8f.0p-26 */
-const P2 : f32 = -2.7667332906e-3;    /* -0xb55215.0p-32 */
+const P1: f32 = 1.6666625440e-1; /*  0xaaaa8f.0p-26 */
+const P2: f32 = -2.7667332906e-3; /* -0xb55215.0p-32 */
 
 #[inline]
 pub fn expf(mut x: f32) -> f32 {
-    let x1p127  = f32::from_bits(0x7f000000); // 0x1p127f === 2 ^ 127
-    let x1p_126 = f32::from_bits(0x800000);   // 0x1p-126f === 2 ^ -126  /*original 0x1p-149f    ??????????? */
-    
+    let x1p127 = f32::from_bits(0x7f000000); // 0x1p127f === 2 ^ 127
+    let x1p_126 = f32::from_bits(0x800000); // 0x1p-126f === 2 ^ -126  /*original 0x1p-149f    ??????????? */
+
     let mut hx = x.to_bits();
-    let sign = (hx >> 31) as i32;   /* sign bit of x */
-    let signb : bool = sign != 0;
-    hx &= 0x7fffffff;  /* high word of |x| */
-    
+    let sign = (hx >> 31) as i32; /* sign bit of x */
+    let signb: bool = sign != 0;
+    hx &= 0x7fffffff; /* high word of |x| */
+
     /* special cases */
-    if hx >= 0x42aeac50 {   /* if |x| >= -87.33655f or NaN */
-        if hx > 0x7f800000 {/* NaN */
+    if hx >= 0x42aeac50 {
+        /* if |x| >= -87.33655f or NaN */
+        if hx > 0x7f800000 {
+            /* NaN */
             return x;
         }
-        if (hx >= 0x42b17218) && (!signb) {  /* x >= 88.722839f */
+        if (hx >= 0x42b17218) && (!signb) {
+            /* x >= 88.722839f */
             /* overflow */
@@ -37,16 +55,19 @@ x *= x1p127;
             /* underflow */
-            force_eval!(-x1p_126/x);
-            if hx >= 0x42cff1b5 { /* x <= -103.972084f */
-                return 0.
+            force_eval!(-x1p_126 / x);
+            if hx >= 0x42cff1b5 {
+                /* x <= -103.972084f */
+                return 0.;
             }
         }
     }
-    
+
     /* argument reduction */
-    let k  : i32;
-    let hi : f32;
-    let lo : f32;
-    if hx > 0x3eb17218 {     /* if |x| > 0.5 ln2 */
-        if hx > 0x3f851592 { /* if |x| > 1.5 ln2 */
-            k = (INV_LN2*x + HALF[sign as usize]) as i32;
+    let k: i32;
+    let hi: f32;
+    let lo: f32;
+    if hx > 0x3eb17218 {
+        /* if |x| > 0.5 ln2 */
+        if hx > 0x3f851592 {
+            /* if |x| > 1.5 ln2 */
+            k = (INV_LN2 * x + HALF[sign as usize]) as i32;
         } else {
@@ -56,6 +77,7 @@ k = 1 - sign - sign;
         let kf = k as f32;
-        hi = x - kf*LN2_HI;  /* k*ln2hi is exact here */
-        lo = kf*LN2_LO;
+        hi = x - kf * LN2_HI; /* k*ln2hi is exact here */
+        lo = kf * LN2_LO;
         x = hi - lo;
-    } else if hx > 0x39000000 {  /* |x| > 2**-14 */
+    } else if hx > 0x39000000 {
+        /* |x| > 2**-14 */
         k = 0;
@@ -69,7 +91,7 @@ hi = x;
     }
-    
+
     /* x is now in primary range */
-    let xx = x*x;
-    let c = x - xx*(P1+xx*P2);
-    let y = 1. + (x*c/(2.-c) - lo + hi);
+    let xx = x * x;
+    let c = x - xx * (P1 + xx * P2);
+    let y = 1. + (x * c / (2. - c) - lo + hi);
     if k == 0 {
@@ -76,0 +98,0 @@ y

src/math/floor.rs

 use core::f64;
 
-const TOINT    : f64 = 1. / f64::EPSILON;
+const TOINT: f64 = 1. / f64::EPSILON;
 
 #[inline]
-pub fn floor(x : f64) -> f64 {
+pub fn floor(x: f64) -> f64 {
     let ui = x.to_bits();
-	let e = ((ui >> 52) & 0x7ff) as i32;
+    let e = ((ui >> 52) & 0x7ff) as i32;
 
-	if (e >= 0x3ff+52) || (x == 0.) {
-		return x;
+    if (e >= 0x3ff + 52) || (x == 0.) {
+        return x;
     }
-	/* y = int(x) - x, where int(x) is an integer neighbor of x */
-	let y = if (ui >> 63) != 0 {
-		x - TOINT + TOINT - x
-	} else {
-		x + TOINT - TOINT - x
+    /* y = int(x) - x, where int(x) is an integer neighbor of x */
+    let y = if (ui >> 63) != 0 {
+        x - TOINT + TOINT - x
+    } else {
+        x + TOINT - TOINT - x
     };
-	/* special case because of non-nearest rounding modes */
-	if e <= 0x3ff-1 {
-		force_eval!(y);
-		return if (ui >> 63) != 0 { -1. } else { 0. };
-	}
-	if y > 0. {
+    /* special case because of non-nearest rounding modes */
+    if e <= 0x3ff - 1 {
+        force_eval!(y);
+        return if (ui >> 63) != 0 { -1. } else { 0. };
+    }
+    if y > 0. {
         x + y - 1.
@@ -26,0 +26,0 @@ } else {

src/math/fmodf.rs

@@ -0,5 +1,4 @@
+use core::f32;
 use core::u32;
 
-use super::isnanf;
-
 #[inline]
@@ -14,3 +13,3 @@ pub fn fmodf(x: f32, y: f32) -> f32 {
 
-    if uyi << 1 == 0 || isnanf(y) || ex == 0xff {
+    if uyi << 1 == 0 || y.is_nan() || ex == 0xff {
         return (x * y) / (x * y);
@@ -17,0 +16,0 @@ }

src/math/hypot.rs

@@ -7,2 +7,3 @@ use core::f64;
 
+#[inline]
 fn sq(x: f64) -> (f64, f64) {
@@ -16,4 +17,4 @@ let xh: f64;
     xl = x - xh;
-    let hi = x*x;
-    let lo = xh*xh - hi + 2.*xh*xl + xl*xl;
+    let hi = x * x;
+    let lo = xh * xh - hi + 2. * xh * xl + xl * xl;
     (hi, lo)
@@ -44,4 +45,4 @@ }
     /* special cases */
-    ex = (uxi>>52) as i64;
-    ey = (uyi>>52) as i64;
+    ex = (uxi >> 52) as i64;
+    ey = (uyi >> 52) as i64;
     x = f64::from_bits(uxi);
@@ -65,7 +66,7 @@ y = f64::from_bits(uyi);
     z = 1.;
-    if ex > 0x3ff+510 {
+    if ex > 0x3ff + 510 {
         z = x1p700;
         x *= x1p_700;
         y *= x1p_700;
-    } else if ey < 0x3ff-450 {
+    } else if ey < 0x3ff - 450 {
         z = x1p_700;
@@ -77,3 +78,3 @@ x *= x1p700;
     let (hy, ly) = sq(y);
-    return z*sqrt(ly+lx+hy+hx);
+    return z * sqrt(ly + lx + hy + hx);
 }

src/math/hypotf.rs

@@ -25,6 +25,6 @@ use core::f32;
     y = f32::from_bits(uyi);
-    if uyi == 0xff<<23 {
+    if uyi == 0xff << 23 {
         return y;
     }
-    if uxi >= 0xff<<23 || uyi == 0 || uxi - uyi >= 25<<23 {
+    if uxi >= 0xff << 23 || uyi == 0 || uxi - uyi >= 25 << 23 {
         return x + y;
@@ -34,7 +34,7 @@ }
     z = 1.;
-    if uxi >= (0x7f+60)<<23 {
+    if uxi >= (0x7f + 60) << 23 {
         z = x1p90;
         x *= x1p_90;
         y *= x1p_90;
-    } else if uyi < (0x7f-60)<<23 {
+    } else if uyi < (0x7f - 60) << 23 {
         z = x1p_90;
@@ -44,3 +44,3 @@ x *= x1p90;
     }
-    z*sqrtf((x as f64 * x as f64 + y as f64 * y as f64) as f32)
+    z * sqrtf((x as f64 * x as f64 + y as f64 * y as f64) as f32)
 }

src/math/logf.rs

@@ -1,8 +0,23 @@
-const LN2_HI : f32 = 6.9313812256e-01;  /* 0x3f317180 */
-const LN2_LO : f32 = 9.0580006145e-06; /* 0x3717f7d1 */
+/* origin: FreeBSD /usr/src/lib/msun/src/e_logf.c */
+/*
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, [e-mail].
+ */
+/*
+ * ====================================================
+ * 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 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 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 */
 
@@ -12,12 +27,13 @@ #[inline]
     let x1p25 = f32::from_bits(0x4c000000); // 0x1p25f === 2 ^ 25
-    
+
     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 < 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 */
+        if (ix >> 31) != 0 {
+            return (x - x) / 0.; /* log(-#) = NaN */
         }
@@ -36,16 +52,16 @@ /* subnormal number, scale up x */
     ix += 0x3f800000 - 0x3f3504f3;
-    k += ((ix>>23) as i32) - 0x7f;
+    k += ((ix >> 23) as i32) - 0x7f;
     ix = (ix & 0x007fffff) + 0x3f3504f3;
-    x =  f32::from_bits(ix);
+    x = f32::from_bits(ix);
 
     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 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 hfsq = 0.5 * f * f;
     let dk = k as f32;
-    s*(hfsq+r) + dk*LN2_LO - hfsq + f + dk*LN2_HI
+    s * (hfsq + r) + dk * LN2_LO - hfsq + f + dk * LN2_HI
 }

src/math/mod.rs

 macro_rules! force_eval {
     ($e:expr) => {
-        unsafe { ::core::ptr::read_volatile(&$e); }
-    }
+        unsafe {
+            ::core::ptr::read_volatile(&$e);
+        }
+    };
 }
 
+mod acos;
+mod acosf;
+mod asin;
+mod asinf;
+mod atanf;
+mod cbrt;
+mod cbrtf;
+mod ceil;
+mod ceilf;
+mod cosf;
+mod exp;
+mod exp2;
+mod exp2f;
+mod expf;
+mod expm1;
+mod expm1f;
 mod fabs;
 mod fabsf;
+mod fdim;
+mod fdimf;
+mod floor;
+mod floorf;
+mod fma;
+mod fmod;
 mod fmodf;
+mod hypot;
+mod hypotf;
+mod log;
+mod log10;
+mod log10f;
+mod log1p;
+mod log1pf;
+mod log2;
+mod log2f;
+mod logf;
 mod powf;
 mod round;
+mod roundf;
 mod scalbn;
 mod scalbnf;
+mod sinf;
 mod sqrt;
 mod sqrtf;
-mod logf;
-mod expf;
-mod floor;
+mod tanf;
 mod trunc;
 mod truncf;
-mod hypot;
-mod hypotf;
 
-//mod service;
+// Use separated imports instead of {}-grouped imports for easier merging.
+pub use self::acos::acos;
+pub use self::acosf::acosf;
+pub use self::asin::asin;
+pub use self::asinf::asinf;
+pub use self::atanf::atanf;
+pub use self::cbrt::cbrt;
+pub use self::cbrtf::cbrtf;
+pub use self::ceil::ceil;
+pub use self::ceilf::ceilf;
+pub use self::cosf::cosf;
+pub use self::exp::exp;
+pub use self::exp2::exp2;
+pub use self::exp2f::exp2f;
+pub use self::expf::expf;
+pub use self::expm1::expm1;
+pub use self::expm1f::expm1f;
+pub use self::fabs::fabs;
+pub use self::fabsf::fabsf;
+pub use self::fdim::fdim;
+pub use self::fdimf::fdimf;
+pub use self::floor::floor;
+pub use self::floorf::floorf;
+pub use self::fma::fma;
+pub use self::fmod::fmod;
+pub use self::fmodf::fmodf;
+pub use self::hypot::hypot;
+pub use self::hypotf::hypotf;
+pub use self::log::log;
+pub use self::log10::log10;
+pub use self::log10f::log10f;
+pub use self::log1p::log1p;
+pub use self::log1pf::log1pf;
+pub use self::log2::log2;
+pub use self::log2f::log2f;
+pub use self::logf::logf;
+pub use self::powf::powf;
+pub use self::round::round;
+pub use self::roundf::roundf;
+pub use self::scalbn::scalbn;
+pub use self::scalbnf::scalbnf;
+pub use self::sinf::sinf;
+pub use self::sqrt::sqrt;
+pub use self::sqrtf::sqrtf;
+pub use self::tanf::tanf;
+pub use self::trunc::trunc;
+pub use self::truncf::truncf;
 
-pub use self::{
-    fabs::fabs,
-    fabsf::fabsf,
-    fmodf::fmodf,
-    powf::powf,
-    round::round,
-    scalbn::scalbn,
-    scalbnf::scalbnf,
-    sqrt::sqrt,
-    sqrtf::sqrtf,
-    logf::logf,
-    expf::expf,
-    floor::floor,
-    trunc::trunc,
-    truncf::truncf,
-    hypot::hypot,
-    hypotf::hypotf,
+mod k_cosf;
+mod k_sinf;
+mod k_tanf;
+mod rem_pio2_large;
+mod rem_pio2f;
+
+use self::{
+    k_cosf::k_cosf, k_sinf::k_sinf, k_tanf::k_tanf, rem_pio2_large::rem_pio2_large,
+    rem_pio2f::rem_pio2f,
 };
 
-fn isnanf(x: f32) -> bool {
-    x.to_bits() & 0x7fffffff > 0x7f800000
+#[inline]
+pub fn get_high_word(x: f64) -> u32 {
+    (x.to_bits() >> 32) as u32
 }
+
+#[inline]
+pub fn get_low_word(x: f64) -> u32 {
+    x.to_bits() as u32
+}
+
+#[inline]
+pub fn with_set_high_word(f: f64, hi: u32) -> f64 {
+    let mut tmp = f.to_bits();
+    tmp &= 0x00000000_ffffffff;
+    tmp |= (hi as u64) << 32;
+    f64::from_bits(tmp)
+}
+
+#[inline]
+pub fn with_set_low_word(f: f64, lo: u32) -> f64 {
+    let mut tmp = f.to_bits();
+    tmp &= 0xffffffff_00000000;
+    tmp |= lo as u64;
+    f64::from_bits(tmp)
+}

src/math/powf.rs

@@ -0,1 +1,16 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/e_powf.c */
+/*
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, [e-mail].
+ */
+/*
+ * ====================================================
+ * 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.
+ * ====================================================
+ */
+
 use super::{fabsf, scalbnf, sqrtf};
@@ -2,0 +17,0 @@

src/math/round.rs

@@ -5,2 +5,3 @@ use core::f64;
 
+#[inline]
 pub fn round(mut x: f64) -> f64 {
@@ -30,3 +31,3 @@ let (f, i) = (x, x.to_bits());
     }
-    
+
     if i >> 63 != 0 {
@@ -33,0 +34,0 @@ -y

src/math/scalbn.rs

 #[inline]
-pub fn scalbn(x : f64, mut n: i32) -> f64 {
-    let x1p1023  = f64::from_bits(0x7fe0000000000000); // 0x1p1023 === 2 ^ 1023
-    let x1p53    = f64::from_bits(0x4340000000000000); // 0x1p53 === 2 ^ 53
+pub fn scalbn(x: f64, mut n: i32) -> f64 {
+    let x1p1023 = f64::from_bits(0x7fe0000000000000); // 0x1p1023 === 2 ^ 1023
+    let x1p53 = f64::from_bits(0x4340000000000000); // 0x1p53 === 2 ^ 53
     let x1p_1022 = f64::from_bits(0x0010000000000000); // 0x1p-1022 === 2 ^ (-1022)
-    
+
     let mut y = x;
-    
+
     if n > 1023 {
@@ -32,3 +32,3 @@ y *= x1p1023;
     }
-    y*f64::from_bits(((0x3ff+n) as u64)<<52)
+    y * f64::from_bits(((0x3ff + n) as u64) << 52)
 }

src/math/scalbnf.rs

 #[inline]
-pub fn scalbnf(mut x: f32, mut n : i32) -> f32 {
-    let x1p127  = f32::from_bits(0x7f000000); // 0x1p127f === 2 ^ 127
-    let x1p_126 = f32::from_bits(0x800000);   // 0x1p-126f === 2 ^ -126
-    let x1p24   = f32::from_bits(0x4b800000); // 0x1p24f === 2 ^ 24
-    
+pub fn scalbnf(mut x: f32, mut n: i32) -> f32 {
+    let x1p127 = f32::from_bits(0x7f000000); // 0x1p127f === 2 ^ 127
+    let x1p_126 = f32::from_bits(0x800000); // 0x1p-126f === 2 ^ -126
+    let x1p24 = f32::from_bits(0x4b800000); // 0x1p24f === 2 ^ 24
+
     if n > 127 {
@@ -28,3 +28,3 @@ x *= x1p127;
     }
-    x * f32::from_bits(((0x7f+n) as u32)<<23)
+    x * f32::from_bits(((0x7f + n) as u32) << 23)
 }

src/math/sqrt.rs

@@ -0,1 +1,79 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/e_sqrt.c */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+/* sqrt(x)
+ * Return correctly rounded sqrt.
+ *           ------------------------------------------
+ *           |  Use the hardware sqrt if you have one |
+ *           ------------------------------------------
+ * Method:
+ *   Bit by bit method using integer arithmetic. (Slow, but portable)
+ *   1. Normalization
+ *      Scale x to y in [1,4) with even powers of 2:
+ *      find an integer k such that  1 <= (y=x*2^(2k)) < 4, then
+ *              sqrt(x) = 2^k * sqrt(y)
+ *   2. Bit by bit computation
+ *      Let q  = sqrt(y) truncated to i bit after binary point (q = 1),
+ *           i                                                   0
+ *                                     i+1         2
+ *          s  = 2*q , and      y  =  2   * ( y - q  ).         (1)
+ *           i      i            i                 i
+ *
+ *      To compute q    from q , one checks whether
+ *                  i+1       i
+ *
+ *                            -(i+1) 2
+ *                      (q + 2      ) <= y.                     (2)
+ *                        i
+ *                                                            -(i+1)
+ *      If (2) is false, then q   = q ; otherwise q   = q  + 2      .
+ *                             i+1   i             i+1   i
+ *
+ *      With some algebric manipulation, it is not difficult to see
+ *      that (2) is equivalent to
+ *                             -(i+1)
+ *                      s  +  2       <= y                      (3)
+ *                       i                i
+ *
+ *      The advantage of (3) is that s  and y  can be computed by
+ *                                    i      i
+ *      the following recurrence formula:
+ *          if (3) is false
+ *
+ *          s     =  s  ,       y    = y   ;                    (4)
+ *           i+1      i          i+1    i
+ *
+ *          otherwise,
+ *                         -i                     -(i+1)
+ *          s     =  s  + 2  ,  y    = y  -  s  - 2             (5)
+ *           i+1      i          i+1    i     i
+ *
+ *      One may easily use induction to prove (4) and (5).
+ *      Note. Since the left hand side of (3) contain only i+2 bits,
+ *            it does not necessary to do a full (53-bit) comparison
+ *            in (3).
+ *   3. Final rounding
+ *      After generating the 53 bits result, we compute one more bit.
+ *      Together with the remainder, we can decide whether the
+ *      result is exact, bigger than 1/2ulp, or less than 1/2ulp
+ *      (it will never equal to 1/2ulp).
+ *      The rounding mode can be detected by checking whether
+ *      huge + tiny is equal to huge, and whether huge - tiny is
+ *      equal to huge for some floating point number "huge" and "tiny".
+ *
+ * Special cases:
+ *      sqrt(+-0) = +-0         ... exact
+ *      sqrt(inf) = inf
+ *      sqrt(-ve) = NaN         ... with invalid signal
+ *      sqrt(NaN) = NaN         ... with invalid signal for signaling NaN
+ */
+
 use core::f64;
@@ -25,24 +103,25 @@
     /* take care of Inf and NaN */
-    if (ix0&0x7ff00000) == 0x7ff00000 {
-        return x*x + x;  /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */
+    if (ix0 & 0x7ff00000) == 0x7ff00000 {
+        return x * x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */
     }
     /* take care of zero */
     if ix0 <= 0 {
-        if ((ix0&!(sign as i32))|ix1 as i32) == 0 {
-            return x;  /* sqrt(+-0) = +-0 */
+        if ((ix0 & !(sign as i32)) | ix1 as i32) == 0 {
+            return x; /* sqrt(+-0) = +-0 */
         }
         if ix0 < 0 {
-            return (x - x) / (x - x);  /* sqrt(-ve) = sNaN */
+            return (x - x) / (x - x); /* sqrt(-ve) = sNaN */
         }
     }
     /* normalize x */
-    m = ix0>>20;
-    if m == 0 {  /* subnormal x */
+    m = ix0 >> 20;
+    if m == 0 {
+        /* subnormal x */
         while ix0 == 0 {
             m -= 21;
-            ix0 |= (ix1>>11) as i32;
+            ix0 |= (ix1 >> 11) as i32;
             ix1 <<= 21;
         }
-        i=0;
-        while (ix0&0x00100000) == 0 {
+        i = 0;
+        while (ix0 & 0x00100000) == 0 {
             i += 1;
@@ -52,15 +131,16 @@ ix0 <<= 1;
         m -= i - 1;
-        ix0 |= (ix1>>(32-i)) as i32;
+        ix0 |= (ix1 >> (32 - i)) as i32;
         ix1 <<= i;
     }
-    m -= 1023;    /* unbias exponent */
-    ix0 = (ix0&0x000fffff)|0x00100000;
-    if (m & 1) == 1 {  /* odd m, double x to make it even */
-        ix0 += ix0 + ((ix1&sign)>>31) as i32;
+    m -= 1023; /* unbias exponent */
+    ix0 = (ix0 & 0x000fffff) | 0x00100000;
+    if (m & 1) == 1 {
+        /* odd m, double x to make it even */
+        ix0 += ix0 + ((ix1 & sign) >> 31) as i32;
         ix1 += ix1;
     }
-    m >>= 1;      /* m = [m/2] */
+    m >>= 1; /* m = [m/2] */
 
     /* generate sqrt(x) bit by bit */
-    ix0 += ix0 + ((ix1&sign)>>31) as i32;
+    ix0 += ix0 + ((ix1 & sign) >> 31) as i32;
     ix1 += ix1;
@@ -71,3 +151,3 @@ q = 0; /* [q,q1] = sqrt(x) */
     s1 = 0;
-    r = 0x00200000;        /* r = moving bit from right to left */
+    r = 0x00200000; /* r = moving bit from right to left */
 
@@ -77,7 +157,7 @@ while r != 0 {
         if t <= ix0 {
-            s0   = t + r as i32;
+            s0 = t + r as i32;
             ix0 -= t;
-            q   += r as i32;
+            q += r as i32;
         }
-        ix0 += ix0 + ((ix1&sign)>>31) as i32;
+        ix0 += ix0 + ((ix1 & sign) >> 31) as i32;
         ix1 += ix1;
@@ -90,6 +170,6 @@ r >>= 1;
         t1 = s1 + r;
-        t  = s0;
+        t = s0;
         if t < ix0 || (t == ix0 && t1 <= ix1) {
             s1 = t1 + r;
-            if (t1&sign) == sign && (s1&sign) == 0 {
+            if (t1 & sign) == sign && (s1 & sign) == 0 {
                 s0 += 1;
@@ -104,3 +184,3 @@ }
         }
-        ix0 += ix0 + ((ix1&sign)>>31) as i32;
+        ix0 += ix0 + ((ix1 & sign) >> 31) as i32;
         ix1 += ix1;
@@ -111,3 +191,3 @@ r >>= 1;
     /* use floating add to find out rounding direction */
-    if (ix0 as u32|ix1) != 0 {
+    if (ix0 as u32 | ix1) != 0 {
         z = 1.0 - TINY; /* raise inexact flag */
@@ -118,3 +198,3 @@ if z >= 1.0 {
                 q1 = 0;
-                q+=1;
+                q += 1;
             } else if z > 1.0 {
@@ -130,5 +210,5 @@ if q1 == 0xfffffffe {
     }
-    ix0 = (q>>1) + 0x3fe00000;
-    ix1 = q1>>1;
-    if (q&1) == 1 {
+    ix0 = (q >> 1) + 0x3fe00000;
+    ix1 = q1 >> 1;
+    if (q & 1) == 1 {
         ix1 |= sign;
@@ -135,0 +215,0 @@ }

src/math/sqrtf.rs

@@ -0,1 +1,16 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/e_sqrtf.c */
+/*
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, [e-mail].
+ */
+/*
+ * ====================================================
+ * 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 TINY: f32 = 1.0e-30;
@@ -2,0 +17,0 @@

No alert changes

No metric changes