diff --git a/src/math/acosh.rs b/src/math/acosh.rs
new file mode 100644
index 0000000..3494e34
--- /dev/null
+++ b/src/math/acosh.rs
@@ -0,0 +1,22 @@
+use super::{log, log1p, sqrt};
+
+const LN2: f64 = 0.693147180559945309417232121458176568; /* 0x3fe62e42,  0xfefa39ef*/
+
+/* acosh(x) = log(x + sqrt(x*x-1)) */
+pub fn acosh(x: f64) -> f64 {
+    let u = x.to_bits();
+    let e = ((u >> 52) as usize) & 0x7ff;
+
+    /* x < 1 domain error is handled in the called functions */
+
+    if e < 0x3ff + 1 {
+        /* |x| < 2, up to 2ulp error in [1,1.125] */
+        return log1p(x-1.0+sqrt((x-1.0)*(x-1.0)+2.0*(x-1.0)));
+    }
+    if e < 0x3ff + 26 {
+        /* |x| < 0x1p26 */
+        return log(2.0*x-1.0/(x+sqrt(x*x-1.0)));
+    }
+    /* |x| >= 0x1p26 or nan */
+    return log(x) + LN2;
+}
diff --git a/src/math/acoshf.rs b/src/math/acoshf.rs
new file mode 100644
index 0000000..1e298a9
--- /dev/null
+++ b/src/math/acoshf.rs
@@ -0,0 +1,21 @@
+use super::{log1pf, logf, sqrtf};
+
+const LN2: f32 = 0.693147180559945309417232121458176568;
+
+/* acosh(x) = log(x + sqrt(x*x-1)) */
+pub fn acoshf(x: f32) -> f32 {
+    let u = x.to_bits();
+    let a = u & 0x7fffffff;
+
+    if a < 0x3f800000+(1<<23) {
+        /* |x| < 2, invalid if x < 1 or nan */
+        /* up to 2ulp error in [1,1.125] */
+        return log1pf(x-1.0 + sqrtf((x-1.0)*(x-1.0)+2.0*(x-1.0)));
+    }
+    if a < 0x3f800000+(12<<23) {
+        /* |x| < 0x1p12 */
+        return logf(2.0*x - 1.0/(x+sqrtf(x*x-1.0)));
+    }
+    /* x >= 0x1p12 */
+    return logf(x) + LN2;
+}
diff --git a/src/math/asinef.rs b/src/math/asinef.rs
new file mode 100644
index 0000000..d2cd826
--- /dev/null
+++ b/src/math/asinef.rs
@@ -0,0 +1,95 @@
+/* @(#)z_asinef.c 1.0 98/08/13 */
+/******************************************************************
+ * The following routines are coded directly from the algorithms
+ * and coefficients given in "Software Manual for the Elementary
+ * Functions" by William J. Cody, Jr. and William Waite, Prentice
+ * Hall, 1980.
+ ******************************************************************/
+/******************************************************************
+ * Arcsine
+ *
+ * Input:
+ *   x - floating point value
+ *   acosine - indicates acos calculation
+ *
+ * Output:
+ *   Arcsine of x.
+ *
+ * Description:
+ *   This routine calculates arcsine / arccosine.
+ *
+ *****************************************************************/
+
+use super::{fabsf, sqrtf};
+
+const P: [f32; 2] = [ 0.933935835, -0.504400557 ];
+const Q: [f32; 2] = [ 0.560363004e+1, -0.554846723e+1 ];
+const A: [f32; 2] = [ 0.0, 0.785398163 ];
+const B: [f32; 2] = [ 1.570796326, 0.785398163 ];
+const Z_ROOTEPS_F: f32 = 1.7263349182589107e-4;
+
+pub fn asinef(x: f32, acosine: usize) -> f32
+{
+    let flag: usize;
+    let i: usize;
+    let mut branch: bool = false;
+    let g: f32;
+    let mut res: f32 = 0.0;
+    let mut y: f32;
+
+    /* Check for special values. */
+    //i = numtestf (x);
+    if x.is_nan() || x.is_infinite() {
+        force_eval!(x);
+        return x;
+    }
+
+    y = fabsf(x);
+    flag = acosine;
+
+    if y > 0.5 {
+        i = 1 - flag;
+
+        /* Check for range error. */
+        if y > 1.0 {
+            return 0.0 / 0.0;
+        }
+
+        g = (1.0 - y) / 2.0;
+        y = -2.0 * sqrtf(g);
+        branch = true;
+    } else {
+        i = flag;
+        if y < Z_ROOTEPS_F {
+            res = y;
+            g = 0.0; // pleasing the uninitialized variable
+        } else {
+            g = y * y;
+        }
+    }
+
+    if y >= Z_ROOTEPS_F || branch {
+        /* Calculate the Taylor series. */
+        let p = (P[1] * g + P[0]) * g;
+        let q = (g + Q[1]) * g + Q[0];
+        let r = p / q;
+
+        res = y + y * r;
+    }
+
+    /* Calculate asine or acose. */
+    if flag == 0 {
+        res = (A[i] + res) + A[i];
+        if x < 0.0 {
+            res = -res;
+        }
+    } else {
+        if x < 0.0 {
+            res = (B[i] + res) + B[i];
+        } else {
+            res = (A[i] - res) + A[i];
+        }
+    }
+
+    return res;
+}
diff --git a/src/math/asinh.rs b/src/math/asinh.rs
new file mode 100644
index 0000000..09e8945
--- /dev/null
+++ b/src/math/asinh.rs
@@ -0,0 +1,35 @@
+use super::{log, log1p, sqrt};
+
+const LN2: f64 = 0.693147180559945309417232121458176568; /* 0x3fe62e42,  0xfefa39ef*/
+
+/* asinh(x) = sign(x)*log(|x|+sqrt(x*x+1)) ~= x - x^3/6 + o(x^5) */
+pub fn asinh(mut x: f64) -> f64 {
+    let mut u = x.to_bits();
+    let e = ((u >> 52) as usize) & 0x7ff;
+    let sign = (u >> 63) != 0;
+
+    /* |x| */
+    u &= (!0) >> 1;
+    x = f64::from_bits(u);
+
+    if e >= 0x3ff + 26 {
+        /* |x| >= 0x1p26 or inf or nan */
+        x = log(x) + LN2;
+    } else if e >= 0x3ff + 1 {
+        /* |x| >= 2 */
+        x = log(2.0*x + 1.0/(sqrt(x*x+1.0)+x));
+    } else if e >= 0x3ff - 26 {
+        /* |x| >= 0x1p-26, up to 1.6ulp error in [0.125,0.5] */
+        x = log1p(x + x*x/(sqrt(x*x+1.0)+1.0));
+    } else {
+        /* |x| < 0x1p-26, raise inexact if x != 0 */
+        let x1p120 = f64::from_bits(0x4770000000000000);
+        force_eval!(x + x1p120);
+    }
+
+    if sign {
+        -x
+    } else {
+        x
+    }
+}
diff --git a/src/math/asinhf.rs b/src/math/asinhf.rs
new file mode 100644
index 0000000..236916d
--- /dev/null
+++ b/src/math/asinhf.rs
@@ -0,0 +1,34 @@
+use super::{logf, log1pf, sqrtf};
+
+const LN2: f32 = 0.693147180559945309417232121458176568;
+
+/* asinh(x) = sign(x)*log(|x|+sqrt(x*x+1)) ~= x - x^3/6 + o(x^5) */
+pub fn asinhf(mut x: f32) -> f32 {
+    let u = x.to_bits();
+    let i = u & 0x7fffffff;
+    let sign = (u >> 31) != 0;
+
+    /* |x| */
+    x = f32::from_bits(i);
+
+    if i >= 0x3f800000 + (12<<23) {
+        /* |x| >= 0x1p12 or inf or nan */
+        x = logf(x) + LN2;
+    } else if i >= 0x3f800000 + (1<<23) {
+        /* |x| >= 2 */
+        x = logf(2.0*x + 1.0/(sqrtf(x*x+1.0)+x));
+    } else if i >= 0x3f800000 - (12<<23) {
+        /* |x| >= 0x1p-12, up to 1.6ulp error in [0.125,0.5] */
+        x = log1pf(x + x*x/(sqrtf(x*x+1.0)+1.0));
+    } else {
+        /* |x| < 0x1p-12, raise inexact if x!=0 */
+        let x1p120 = f32::from_bits(0x7b800000);
+        force_eval!(x + x1p120);
+    }
+
+    if sign {
+        -x
+    } else {
+        x
+    }
+}
diff --git a/src/math/atanh.rs b/src/math/atanh.rs
new file mode 100644
index 0000000..ea44480
--- /dev/null
+++ b/src/math/atanh.rs
@@ -0,0 +1,33 @@
+use super::{log1p};
+
+/* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */
+pub fn atanh(mut x: f64) -> f64 {
+    let mut u = x.to_bits();
+    let e = ((u >> 52) as usize) & 0x7ff;
+    let sign = (u >> 63) != 0;
+
+    /* |x| */
+    u &= 0x7fffffff;
+    x = f64::from_bits(u);
+
+    if e < 0x3ff - 1 {
+        if e < 0x3ff - 32 {
+            /* handle underflow */
+            if e == 0 {
+                force_eval!(x as f32);
+            }
+        } else {
+            /* |x| < 0.5, up to 1.7ulp error */
+            x = 0.5*log1p(2.0*x + 2.0*x*x/(1.0-x));
+        }
+    } else {
+        /* avoid overflow */
+        x = 0.5*log1p(2.0*(x/(1.0-x)));
+    }
+
+    if sign {
+        -x
+    } else {
+        x
+    }
+}
diff --git a/src/math/atanhf.rs b/src/math/atanhf.rs
new file mode 100644
index 0000000..77d451b
--- /dev/null
+++ b/src/math/atanhf.rs
@@ -0,0 +1,32 @@
+use super::{log1pf};
+
+/* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */
+pub fn atanhf(mut x: f32) -> f32 {
+    let mut u = x.to_bits();
+    let sign = (u >> 31) != 0;
+
+    /* |x| */
+    u &= 0x7fffffff;
+    x = f32::from_bits(u);
+
+    if u < 0x3f800000 - (1<<23) {
+        if u < 0x3f800000 - (32<<23) {
+            /* handle underflow */
+            if u < (1<<23) {
+                force_eval!((x*x) as f32);
+            }
+        } else {
+            /* |x| < 0.5, up to 1.7ulp error */
+            x = 0.5*log1pf(2.0*x + 2.0*x*x/(1.0-x));
+        }
+    } else {
+        /* avoid overflow */
+        x = 0.5*log1pf(2.0*(x/(1.0-x)));
+    }
+
+    if sign {
+        -x
+    } else {
+        x
+    }
+}
diff --git a/src/math/copysign.rs b/src/math/copysign.rs
new file mode 100644
index 0000000..74b761e
--- /dev/null
+++ b/src/math/copysign.rs
@@ -0,0 +1,7 @@
+pub fn copysign(x: f64, y: f64) -> f64 {
+    let mut ux = x.to_bits();
+    let uy = y.to_bits();
+    ux &= (!0) >> 1;
+    ux |= uy & (1<<63);
+    f64::from_bits(ux)
+}
diff --git a/src/math/copysignf.rs b/src/math/copysignf.rs
new file mode 100644
index 0000000..a0a814b
--- /dev/null
+++ b/src/math/copysignf.rs
@@ -0,0 +1,7 @@
+pub fn copysignf(x: f32, y: f32) -> f32 {
+    let mut ux = x.to_bits();
+    let uy = y.to_bits();
+    ux &= 0x7fffffff;
+    ux |= uy & 0x80000000;
+    f32::from_bits(ux)
+}
diff --git a/src/math/erf.rs b/src/math/erf.rs
new file mode 100644
index 0000000..b3ad2ce
--- /dev/null
+++ b/src/math/erf.rs
@@ -0,0 +1,297 @@
+use super::{exp, fabs, get_high_word, with_set_low_word};
+/* origin: FreeBSD /usr/src/lib/msun/src/s_erf.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 erf(double x)
+ * double erfc(double x)
+ *                           x
+ *                    2      |\
+ *     erf(x)  =  ---------  | exp(-t*t)dt
+ *                 sqrt(pi) \|
+ *                           0
+ *
+ *     erfc(x) =  1-erf(x)
+ *  Note that
+ *              erf(-x) = -erf(x)
+ *              erfc(-x) = 2 - erfc(x)
+ *
+ * Method:
+ *      1. For |x| in [0, 0.84375]
+ *          erf(x)  = x + x*R(x^2)
+ *          erfc(x) = 1 - erf(x)           if x in [-.84375,0.25]
+ *                  = 0.5 + ((0.5-x)-x*R)  if x in [0.25,0.84375]
+ *         where R = P/Q where P is an odd poly of degree 8 and
+ *         Q is an odd poly of degree 10.
+ *                                               -57.90
+ *                      | R - (erf(x)-x)/x | <= 2
+ *
+ *
+ *         Remark. The formula is derived by noting
+ *          erf(x) = (2/sqrt(pi))*(x - x^3/3 + x^5/10 - x^7/42 + ....)
+ *         and that
+ *          2/sqrt(pi) = 1.128379167095512573896158903121545171688
+ *         is close to one. The interval is chosen because the fix
+ *         point of erf(x) is near 0.6174 (i.e., erf(x)=x when x is
+ *         near 0.6174), and by some experiment, 0.84375 is chosen to
+ *         guarantee the error is less than one ulp for erf.
+ *
+ *      2. For |x| in [0.84375,1.25], let s = |x| - 1, and
+ *         c = 0.84506291151 rounded to single (24 bits)
+ *              erf(x)  = sign(x) * (c  + P1(s)/Q1(s))
+ *              erfc(x) = (1-c)  - P1(s)/Q1(s) if x > 0
+ *                        1+(c+P1(s)/Q1(s))    if x < 0
+ *              |P1/Q1 - (erf(|x|)-c)| <= 2**-59.06
+ *         Remark: here we use the taylor series expansion at x=1.
+ *              erf(1+s) = erf(1) + s*Poly(s)
+ *                       = 0.845.. + P1(s)/Q1(s)
+ *         That is, we use rational approximation to approximate
+ *                      erf(1+s) - (c = (single)0.84506291151)
+ *         Note that |P1/Q1|< 0.078 for x in [0.84375,1.25]
+ *         where
+ *              P1(s) = degree 6 poly in s
+ *              Q1(s) = degree 6 poly in s
+ *
+ *      3. For x in [1.25,1/0.35(~2.857143)],
+ *              erfc(x) = (1/x)*exp(-x*x-0.5625+R1/S1)
+ *              erf(x)  = 1 - erfc(x)
+ *         where
+ *              R1(z) = degree 7 poly in z, (z=1/x^2)
+ *              S1(z) = degree 8 poly in z
+ *
+ *      4. For x in [1/0.35,28]
+ *              erfc(x) = (1/x)*exp(-x*x-0.5625+R2/S2) if x > 0
+ *                      = 2.0 - (1/x)*exp(-x*x-0.5625+R2/S2) if -6<x<0
+ *                      = 2.0 - tiny            (if x <= -6)
+ *              erf(x)  = sign(x)*(1.0 - erfc(x)) if x < 6, else
+ *              erf(x)  = sign(x)*(1.0 - tiny)
+ *         where
+ *              R2(z) = degree 6 poly in z, (z=1/x^2)
+ *              S2(z) = degree 7 poly in z
+ *
+ *      Note1:
+ *         To compute exp(-x*x-0.5625+R/S), let s be a single
+ *         precision number and s := x; then
+ *              -x*x = -s*s + (s-x)*(s+x)
+ *              exp(-x*x-0.5626+R/S) =
+ *                      exp(-s*s-0.5625)*exp((s-x)*(s+x)+R/S);
+ *      Note2:
+ *         Here 4 and 5 make use of the asymptotic series
+ *                        exp(-x*x)
+ *              erfc(x) ~ ---------- * ( 1 + Poly(1/x^2) )
+ *                        x*sqrt(pi)
+ *         We use rational approximation to approximate
+ *              g(s)=f(1/x^2) = log(erfc(x)*x) - x*x + 0.5625
+ *         Here is the error bound for R1/S1 and R2/S2
+ *              |R1/S1 - f(x)|  < 2**(-62.57)
+ *              |R2/S2 - f(x)|  < 2**(-61.52)
+ *
+ *      5. For inf > x >= 28
+ *              erf(x)  = sign(x) *(1 - tiny)  (raise inexact)
+ *              erfc(x) = tiny*tiny (raise underflow) if x > 0
+ *                      = 2 - tiny if x<0
+ *
+ *      7. Special case:
+ *              erf(0)  = 0, erf(inf)  = 1, erf(-inf) = -1,
+ *              erfc(0) = 1, erfc(inf) = 0, erfc(-inf) = 2,
+ *              erfc/erf(NaN) is NaN
+ */
+
+const ERX: f64  = 8.45062911510467529297e-01; /* 0x3FEB0AC1, 0x60000000 */
+/*
+ * Coefficients for approximation to  erf on [0,0.84375]
+ */
+const EFX8: f64 =  1.02703333676410069053e+00; /* 0x3FF06EBA, 0x8214DB69 */
+const PP0: f64  =  1.28379167095512558561e-01; /* 0x3FC06EBA, 0x8214DB68 */
+const PP1: f64  = -3.25042107247001499370e-01; /* 0xBFD4CD7D, 0x691CB913 */
+const PP2: f64  = -2.84817495755985104766e-02; /* 0xBF9D2A51, 0xDBD7194F */
+const PP3: f64  = -5.77027029648944159157e-03; /* 0xBF77A291, 0x236668E4 */
+const PP4: f64  = -2.37630166566501626084e-05; /* 0xBEF8EAD6, 0x120016AC */
+const QQ1: f64  =  3.97917223959155352819e-01; /* 0x3FD97779, 0xCDDADC09 */
+const QQ2: f64  =  6.50222499887672944485e-02; /* 0x3FB0A54C, 0x5536CEBA */
+const QQ3: f64  =  5.08130628187576562776e-03; /* 0x3F74D022, 0xC4D36B0F */
+const QQ4: f64  =  1.32494738004321644526e-04; /* 0x3F215DC9, 0x221C1A10 */
+const QQ5: f64  = -3.96022827877536812320e-06; /* 0xBED09C43, 0x42A26120 */
+/*
+ * Coefficients for approximation to  erf  in [0.84375,1.25]
+ */
+const PA0: f64  = -2.36211856075265944077e-03; /* 0xBF6359B8, 0xBEF77538 */
+const PA1: f64  =  4.14856118683748331666e-01; /* 0x3FDA8D00, 0xAD92B34D */
+const PA2: f64  = -3.72207876035701323847e-01; /* 0xBFD7D240, 0xFBB8C3F1 */
+const PA3: f64  =  3.18346619901161753674e-01; /* 0x3FD45FCA, 0x805120E4 */
+const PA4: f64  = -1.10894694282396677476e-01; /* 0xBFBC6398, 0x3D3E28EC */
+const PA5: f64  =  3.54783043256182359371e-02; /* 0x3FA22A36, 0x599795EB */
+const PA6: f64  = -2.16637559486879084300e-03; /* 0xBF61BF38, 0x0A96073F */
+const QA1: f64  =  1.06420880400844228286e-01; /* 0x3FBB3E66, 0x18EEE323 */
+const QA2: f64  =  5.40397917702171048937e-01; /* 0x3FE14AF0, 0x92EB6F33 */
+const QA3: f64  =  7.18286544141962662868e-02; /* 0x3FB2635C, 0xD99FE9A7 */
+const QA4: f64  =  1.26171219808761642112e-01; /* 0x3FC02660, 0xE763351F */
+const QA5: f64  =  1.36370839120290507362e-02; /* 0x3F8BEDC2, 0x6B51DD1C */
+const QA6: f64  =  1.19844998467991074170e-02; /* 0x3F888B54, 0x5735151D */
+/*
+ * Coefficients for approximation to  erfc in [1.25,1/0.35]
+ */
+const RA0: f64  = -9.86494403484714822705e-03; /* 0xBF843412, 0x600D6435 */
+const RA1: f64  = -6.93858572707181764372e-01; /* 0xBFE63416, 0xE4BA7360 */
+const RA2: f64  = -1.05586262253232909814e+01; /* 0xC0251E04, 0x41B0E726 */
+const RA3: f64  = -6.23753324503260060396e+01; /* 0xC04F300A, 0xE4CBA38D */
+const RA4: f64  = -1.62396669462573470355e+02; /* 0xC0644CB1, 0x84282266 */
+const RA5: f64  = -1.84605092906711035994e+02; /* 0xC067135C, 0xEBCCABB2 */
+const RA6: f64  = -8.12874355063065934246e+01; /* 0xC0545265, 0x57E4D2F2 */
+const RA7: f64  = -9.81432934416914548592e+00; /* 0xC023A0EF, 0xC69AC25C */
+const SA1: f64  =  1.96512716674392571292e+01; /* 0x4033A6B9, 0xBD707687 */
+const SA2: f64  =  1.37657754143519042600e+02; /* 0x4061350C, 0x526AE721 */
+const SA3: f64  =  4.34565877475229228821e+02; /* 0x407B290D, 0xD58A1A71 */
+const SA4: f64  =  6.45387271733267880336e+02; /* 0x40842B19, 0x21EC2868 */
+const SA5: f64  =  4.29008140027567833386e+02; /* 0x407AD021, 0x57700314 */
+const SA6: f64  =  1.08635005541779435134e+02; /* 0x405B28A3, 0xEE48AE2C */
+const SA7: f64  =  6.57024977031928170135e+00; /* 0x401A47EF, 0x8E484A93 */
+const SA8: f64  = -6.04244152148580987438e-02; /* 0xBFAEEFF2, 0xEE749A62 */
+/*
+ * Coefficients for approximation to  erfc in [1/.35,28]
+ */
+const RB0: f64  = -9.86494292470009928597e-03; /* 0xBF843412, 0x39E86F4A */
+const RB1: f64  = -7.99283237680523006574e-01; /* 0xBFE993BA, 0x70C285DE */
+const RB2: f64  = -1.77579549177547519889e+01; /* 0xC031C209, 0x555F995A */
+const RB3: f64  = -1.60636384855821916062e+02; /* 0xC064145D, 0x43C5ED98 */
+const RB4: f64  = -6.37566443368389627722e+02; /* 0xC083EC88, 0x1375F228 */
+const RB5: f64  = -1.02509513161107724954e+03; /* 0xC0900461, 0x6A2E5992 */
+const RB6: f64  = -4.83519191608651397019e+02; /* 0xC07E384E, 0x9BDC383F */
+const SB1: f64  =  3.03380607434824582924e+01; /* 0x403E568B, 0x261D5190 */
+const SB2: f64  =  3.25792512996573918826e+02; /* 0x40745CAE, 0x221B9F0A */
+const SB3: f64  =  1.53672958608443695994e+03; /* 0x409802EB, 0x189D5118 */
+const SB4: f64  =  3.19985821950859553908e+03; /* 0x40A8FFB7, 0x688C246A */
+const SB5: f64  =  2.55305040643316442583e+03; /* 0x40A3F219, 0xCEDF3BE6 */
+const SB6: f64  =  4.74528541206955367215e+02; /* 0x407DA874, 0xE79FE763 */
+const SB7: f64  = -2.24409524465858183362e+01; /* 0xC03670E2, 0x42712D62 */
+
+fn erfc1(x: f64) -> f64 {
+    let s: f64;
+    let p: f64;
+    let q: f64;
+
+    s = fabs(x) - 1.0;
+    p = PA0+s*(PA1+s*(PA2+s*(PA3+s*(PA4+s*(PA5+s*PA6)))));
+    q = 1.0+s*(QA1+s*(QA2+s*(QA3+s*(QA4+s*(QA5+s*QA6)))));
+
+    1.0 - ERX - p/q
+}
+
+fn erfc2(ix: u32, mut x: f64) -> f64 {
+    let s: f64;
+    let r: f64;
+    let big_s: f64;
+    let z: f64;
+
+    if ix < 0x3ff40000 {  /* |x| < 1.25 */
+        return erfc1(x);
+    }
+
+    x = fabs(x);
+    s = 1.0/(x*x);
+    if ix < 0x4006db6d {  /* |x| < 1/.35 ~ 2.85714 */
+        r = RA0+s*(RA1+s*(RA2+s*(RA3+s*(RA4+s*(
+             RA5+s*(RA6+s*RA7))))));
+        big_s = 1.0+s*(SA1+s*(SA2+s*(SA3+s*(SA4+s*(
+             SA5+s*(SA6+s*(SA7+s*SA8)))))));
+    } else {                /* |x| > 1/.35 */
+        r = RB0+s*(RB1+s*(RB2+s*(RB3+s*(RB4+s*(
+             RB5+s*RB6)))));
+        big_s = 1.0+s*(SB1+s*(SB2+s*(SB3+s*(SB4+s*(
+             SB5+s*(SB6+s*SB7))))));
+    }
+    z = with_set_low_word(x, 0);
+
+    exp(-z*z-0.5625)*exp((z-x)*(z+x)+r/big_s)/x
+}
+
+pub fn erf(x: f64) -> f64 {
+    let r: f64;
+    let s: f64;
+    let z: f64;
+    let y: f64;
+    let mut ix: u32;
+    let sign: usize;
+
+    ix = get_high_word(x);
+    sign = (ix>>31) as usize;
+    ix &= 0x7fffffff;
+    if ix >= 0x7ff00000 {
+        /* erf(nan)=nan, erf(+-inf)=+-1 */
+        return 1.0-2.0*(sign as f64) + 1.0/x;
+    }
+    if ix < 0x3feb0000 {  /* |x| < 0.84375 */
+        if ix < 0x3e300000 {  /* |x| < 2**-28 */
+            /* avoid underflow */
+            return 0.125*(8.0*x + EFX8*x);
+        }
+        z = x*x;
+        r = PP0+z*(PP1+z*(PP2+z*(PP3+z*PP4)));
+        s = 1.0+z*(QQ1+z*(QQ2+z*(QQ3+z*(QQ4+z*QQ5))));
+        y = r/s;
+        return x + x*y;
+    }
+    if ix < 0x40180000 { /* 0.84375 <= |x| < 6 */
+        y = 1.0 - erfc2(ix,x);
+    } else {
+        let x1p_1022 = f64::from_bits(0x0010000000000000);
+        y = 1.0 - x1p_1022;
+    }
+
+    if sign != 0 {
+        -y
+    } else {
+        y
+    }
+}
+
+pub fn erfc(x: f64) -> f64 {
+    let r: f64;
+    let s: f64;
+    let z: f64;
+    let y: f64;
+    let mut ix: u32;
+    let sign: usize;
+
+    ix = get_high_word(x);
+    sign = (ix>>31) as usize;
+    ix &= 0x7fffffff;
+    if ix >= 0x7ff00000 {
+        /* erfc(nan)=nan, erfc(+-inf)=0,2 */
+        return 2.0*(sign as f64) + 1.0/x;
+    }
+    if ix < 0x3feb0000 {  /* |x| < 0.84375 */
+        if ix < 0x3c700000 { /* |x| < 2**-56 */
+            return 1.0 - x;
+        }
+        z = x*x;
+        r = PP0+z*(PP1+z*(PP2+z*(PP3+z*PP4)));
+        s = 1.0+z*(QQ1+z*(QQ2+z*(QQ3+z*(QQ4+z*QQ5))));
+        y = r/s;
+        if sign != 0 || ix < 0x3fd00000 {  /* x < 1/4 */
+            return 1.0 - (x+x*y);
+        }
+        return 0.5 - (x - 0.5 + x*y);
+    }
+    if ix < 0x403c0000 {  /* 0.84375 <= |x| < 28 */
+        if sign != 0 {
+            return 2.0 - erfc2(ix,x);
+        } else {
+            return erfc2(ix,x);
+        }
+    }
+
+    let x1p_1022 = f64::from_bits(0x0010000000000000);
+    if sign != 0 {
+        2.0 - x1p_1022
+    } else {
+        x1p_1022*x1p_1022
+    }
+}
diff --git a/src/math/erff.rs b/src/math/erff.rs
new file mode 100644
index 0000000..0aaa897
--- /dev/null
+++ b/src/math/erff.rs
@@ -0,0 +1,210 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_erff.c */
+/*
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+/*
+ * ====================================================
+ * 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::{expf, fabsf};
+
+const ERX: f32  =  8.4506291151e-01; /* 0x3f58560b */
+/*
+ * Coefficients for approximation to  erf on [0,0.84375]
+ */
+const EFX8: f32 =  1.0270333290e+00; /* 0x3f8375d4 */
+const PP0: f32  =  1.2837916613e-01; /* 0x3e0375d4 */
+const PP1: f32  = -3.2504209876e-01; /* 0xbea66beb */
+const PP2: f32  = -2.8481749818e-02; /* 0xbce9528f */
+const PP3: f32  = -5.7702702470e-03; /* 0xbbbd1489 */
+const PP4: f32  = -2.3763017452e-05; /* 0xb7c756b1 */
+const QQ1: f32  =  3.9791721106e-01; /* 0x3ecbbbce */
+const QQ2: f32  =  6.5022252500e-02; /* 0x3d852a63 */
+const QQ3: f32  =  5.0813062117e-03; /* 0x3ba68116 */
+const QQ4: f32  =  1.3249473704e-04; /* 0x390aee49 */
+const QQ5: f32  = -3.9602282413e-06; /* 0xb684e21a */
+/*
+ * Coefficients for approximation to  erf  in [0.84375,1.25]
+ */
+const PA0: f32  = -2.3621185683e-03; /* 0xbb1acdc6 */
+const PA1: f32  =  4.1485610604e-01; /* 0x3ed46805 */
+const PA2: f32  = -3.7220788002e-01; /* 0xbebe9208 */
+const PA3: f32  =  3.1834661961e-01; /* 0x3ea2fe54 */
+const PA4: f32  = -1.1089469492e-01; /* 0xbde31cc2 */
+const PA5: f32  =  3.5478305072e-02; /* 0x3d1151b3 */
+const PA6: f32  = -2.1663755178e-03; /* 0xbb0df9c0 */
+const QA1: f32  =  1.0642088205e-01; /* 0x3dd9f331 */
+const QA2: f32  =  5.4039794207e-01; /* 0x3f0a5785 */
+const QA3: f32  =  7.1828655899e-02; /* 0x3d931ae7 */
+const QA4: f32  =  1.2617121637e-01; /* 0x3e013307 */
+const QA5: f32  =  1.3637083583e-02; /* 0x3c5f6e13 */
+const QA6: f32  =  1.1984500103e-02; /* 0x3c445aa3 */
+/*
+ * Coefficients for approximation to  erfc in [1.25,1/0.35]
+ */
+const RA0: f32  = -9.8649440333e-03; /* 0xbc21a093 */
+const RA1: f32  = -6.9385856390e-01; /* 0xbf31a0b7 */
+const RA2: f32  = -1.0558626175e+01; /* 0xc128f022 */
+const RA3: f32  = -6.2375331879e+01; /* 0xc2798057 */
+const RA4: f32  = -1.6239666748e+02; /* 0xc322658c */
+const RA5: f32  = -1.8460508728e+02; /* 0xc3389ae7 */
+const RA6: f32  = -8.1287437439e+01; /* 0xc2a2932b */
+const RA7: f32  = -9.8143291473e+00; /* 0xc11d077e */
+const SA1: f32  =  1.9651271820e+01; /* 0x419d35ce */
+const SA2: f32  =  1.3765776062e+02; /* 0x4309a863 */
+const SA3: f32  =  4.3456588745e+02; /* 0x43d9486f */
+const SA4: f32  =  6.4538726807e+02; /* 0x442158c9 */
+const SA5: f32  =  4.2900814819e+02; /* 0x43d6810b */
+const SA6: f32  =  1.0863500214e+02; /* 0x42d9451f */
+const SA7: f32  =  6.5702495575e+00; /* 0x40d23f7c */
+const SA8: f32  = -6.0424413532e-02; /* 0xbd777f97 */
+/*
+ * Coefficients for approximation to  erfc in [1/.35,28]
+ */
+const RB0: f32  = -9.8649431020e-03; /* 0xbc21a092 */
+const RB1: f32  = -7.9928326607e-01; /* 0xbf4c9dd4 */
+const RB2: f32  = -1.7757955551e+01; /* 0xc18e104b */
+const RB3: f32  = -1.6063638306e+02; /* 0xc320a2ea */
+const RB4: f32  = -6.3756646729e+02; /* 0xc41f6441 */
+const RB5: f32  = -1.0250950928e+03; /* 0xc480230b */
+const RB6: f32  = -4.8351919556e+02; /* 0xc3f1c275 */
+const SB1: f32  =  3.0338060379e+01; /* 0x41f2b459 */
+const SB2: f32  =  3.2579251099e+02; /* 0x43a2e571 */
+const SB3: f32  =  1.5367296143e+03; /* 0x44c01759 */
+const SB4: f32  =  3.1998581543e+03; /* 0x4547fdbb */
+const SB5: f32  =  2.5530502930e+03; /* 0x451f90ce */
+const SB6: f32  =  4.7452853394e+02; /* 0x43ed43a7 */
+const SB7: f32  = -2.2440952301e+01; /* 0xc1b38712 */
+
+fn erfc1(x: f32) -> f32 {
+    let s: f32;
+    let p: f32;
+    let q: f32;
+
+    s = fabsf(x) - 1.0;
+    p = PA0+s*(PA1+s*(PA2+s*(PA3+s*(PA4+s*(PA5+s*PA6)))));
+    q = 1.0+s*(QA1+s*(QA2+s*(QA3+s*(QA4+s*(QA5+s*QA6)))));
+    return 1.0 - ERX - p/q;
+}
+
+fn erfc2(mut ix: u32, mut x: f32) -> f32 {
+    let s: f32;
+    let r: f32;
+    let big_s: f32;
+    let z: f32;
+
+    if ix < 0x3fa00000 { /* |x| < 1.25 */
+        return erfc1(x);
+    }
+
+    x = fabsf(x);
+    s = 1.0/(x*x);
+    if ix < 0x4036db6d {   /* |x| < 1/0.35 */
+        r = RA0+s*(RA1+s*(RA2+s*(RA3+s*(RA4+s*(
+             RA5+s*(RA6+s*RA7))))));
+        big_s = 1.0+s*(SA1+s*(SA2+s*(SA3+s*(SA4+s*(
+             SA5+s*(SA6+s*(SA7+s*SA8)))))));
+    } else {                 /* |x| >= 1/0.35 */
+        r = RB0+s*(RB1+s*(RB2+s*(RB3+s*(RB4+s*(
+             RB5+s*RB6)))));
+        big_s = 1.0+s*(SB1+s*(SB2+s*(SB3+s*(SB4+s*(
+             SB5+s*(SB6+s*SB7))))));
+    }
+    ix = x.to_bits();
+    z = f32::from_bits(ix&0xffffe000);
+
+    expf(-z*z - 0.5625) * expf((z-x)*(z+x) + r/big_s)/x
+}
+
+pub fn erff(x: f32) -> f32
+{
+    let r: f32;
+    let s: f32;
+    let z: f32;
+    let y: f32;
+    let mut ix: u32;
+    let sign: usize;
+
+    ix = x.to_bits();
+    sign = (ix>>31) as usize;
+    ix &= 0x7fffffff;
+    if ix >= 0x7f800000 {
+        /* erf(nan)=nan, erf(+-inf)=+-1 */
+        return 1.0-2.0*(sign as f32) + 1.0/x;
+    }
+    if ix < 0x3f580000 {  /* |x| < 0.84375 */
+        if ix < 0x31800000 {  /* |x| < 2**-28 */
+            /*avoid underflow */
+            return 0.125*(8.0*x + EFX8*x);
+        }
+        z = x*x;
+        r = PP0+z*(PP1+z*(PP2+z*(PP3+z*PP4)));
+        s = 1.0+z*(QQ1+z*(QQ2+z*(QQ3+z*(QQ4+z*QQ5))));
+        y = r/s;
+        return x + x*y;
+    }
+    if ix < 0x40c00000 {  /* |x| < 6 */
+        y = 1.0 - erfc2(ix,x);
+    } else {
+        let x1p_120 = f32::from_bits(0x03800000);
+        y = 1.0 - x1p_120;
+    }
+
+    if sign != 0 {
+        -y
+    } else {
+        y
+    }
+}
+
+pub fn erfcf(x: f32) -> f32 {
+    let r: f32;
+    let s: f32;
+    let z: f32;
+    let y: f32;
+    let mut ix: u32;
+    let sign: usize;
+
+    ix = x.to_bits();
+    sign = (ix>>31) as usize;
+    ix &= 0x7fffffff;
+    if ix >= 0x7f800000 {
+        /* erfc(nan)=nan, erfc(+-inf)=0,2 */
+        return 2.0*(sign as f32) + 1.0/x;
+    }
+
+    if ix < 0x3f580000 {  /* |x| < 0.84375 */
+        if ix < 0x23800000 { /* |x| < 2**-56 */
+            return 1.0 - x;
+        }
+        z = x*x;
+        r = PP0+z*(PP1+z*(PP2+z*(PP3+z*PP4)));
+        s = 1.0+z*(QQ1+z*(QQ2+z*(QQ3+z*(QQ4+z*QQ5))));
+        y = r/s;
+        if sign != 0 || ix < 0x3e800000 {  /* x < 1/4 */
+            return 1.0 - (x+x*y);
+        }
+        return 0.5 - (x - 0.5 + x*y);
+    }
+    if ix < 0x41e00000 {  /* |x| < 28 */
+        if sign != 0 {
+            return 2.0 - erfc2(ix, x);
+        } else {
+            return erfc2(ix, x);
+        }
+    }
+
+    let x1p_120 = f32::from_bits(0x03800000);
+    if sign != 0 {
+        2.0 - x1p_120
+    } else {
+        x1p_120*x1p_120
+    }
+}
diff --git a/src/math/exp10.rs b/src/math/exp10.rs
new file mode 100644
index 0000000..d12fa0b
--- /dev/null
+++ b/src/math/exp10.rs
@@ -0,0 +1,24 @@
+use super::{exp2, modf, pow};
+
+const LN10: f64 = 3.32192809488736234787031942948939;
+const P10: &[f64] = &[
+    1e-15, 1e-14, 1e-13, 1e-12, 1e-11, 1e-10,
+    1e-9, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1,
+    1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,
+    1e10, 1e11, 1e12, 1e13, 1e14, 1e15
+];
+
+pub fn exp10(x: f64) -> f64
+{
+    let (mut y, n) = modf(x);
+    let u: u64 = n.to_bits();
+    /* fabs(n) < 16 without raising invalid on nan */
+    if (u>>52 & 0x7ff) < 0x3ff+4 {
+        if y == 0.0 {
+            return P10[((n as isize) + 15) as usize];
+        }
+        y = exp2(LN10 * y);
+        return y * P10[((n as isize) + 15) as usize];
+    }
+    return pow(10.0, x);
+}
diff --git a/src/math/exp10f.rs b/src/math/exp10f.rs
new file mode 100644
index 0000000..8fb88a5
--- /dev/null
+++ b/src/math/exp10f.rs
@@ -0,0 +1,22 @@
+use super::{exp2, exp2f, modff};
+
+const LN10_F32: f32 = 3.32192809488736234787031942948939;
+const LN10_F64: f64 = 3.32192809488736234787031942948939;
+const P10: &[f32] = &[
+    1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1,
+    1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7
+];
+
+pub fn exp10f(x: f32) -> f32 {
+    let (mut y, n) = modff(x);
+    let u = n.to_bits();
+    /* fabsf(n) < 8 without raising invalid on nan */
+    if (u>>23 & 0xff) < 0x7f+3 {
+        if y == 0.0 {
+            return P10[((n as isize) + 7) as usize]
+        }
+        y = exp2f(LN10_F32 * y);
+        return y * P10[((n as isize) + 7) as usize];
+    }
+    return exp2(LN10_F64 * (x as f64)) as f32;
+}
diff --git a/src/math/frexp.rs b/src/math/frexp.rs
new file mode 100644
index 0000000..45733a3
--- /dev/null
+++ b/src/math/frexp.rs
@@ -0,0 +1,20 @@
+pub fn frexp(x: f64) -> (f64, isize) {
+    let mut y = x.to_bits();
+    let ee = ((y>>52) & 0x7ff) as isize;
+
+    if ee == 0 {
+        if x != 0.0 {
+            let x1p64 = f64::from_bits(0x43f0000000000000);
+            let (x, e) = frexp(x*x1p64);
+            return (x, e - 64);
+        }
+        return (x, 0);
+    } else if ee == 0x7ff {
+        return (x, 0);
+    }
+
+    let e = ee - 0x3fe;
+    y &= 0x800fffffffffffff;
+    y |= 0x3fe0000000000000;
+    return (f64::from_bits(y), e);
+}
diff --git a/src/math/frexpf.rs b/src/math/frexpf.rs
new file mode 100644
index 0000000..1c9dae0
--- /dev/null
+++ b/src/math/frexpf.rs
@@ -0,0 +1,21 @@
+pub fn frexpf(x: f32) -> (f32, isize) {
+    let mut y = x.to_bits();
+    let ee: isize = ((y>>23) & 0xff) as isize;
+
+    if ee == 0 {
+        if x != 0.0 {
+            let x1p64 = f32::from_bits(0x5f800000);
+            let (x, e) = frexpf(x*x1p64);
+            return (x, e - 64);
+        } else {
+            return (x, 0);
+        }
+    } else if ee == 0xff {
+        return (x, 0);
+    }
+
+    let e = ee - 0x7e;
+    y &= 0x807fffff;
+    y |= 0x3f000000;
+    return (f32::from_bits(y), e);
+} 
diff --git a/src/math/ilogb.rs b/src/math/ilogb.rs
new file mode 100644
index 0000000..78fe030
--- /dev/null
+++ b/src/math/ilogb.rs
@@ -0,0 +1,31 @@
+const FP_ILOGBNAN: isize = -1 - (((!0) >> 1));
+const FP_ILOGB0: isize = FP_ILOGBNAN;
+
+pub fn ilogb(x: f64) -> isize {
+    let mut i: u64 = x.to_bits();
+    let e = ((i>>52) & 0x7ff) as isize;
+
+    if e == 0 {
+        i <<= 12;
+        if i == 0 {
+            force_eval!(0.0/0.0);
+            return FP_ILOGB0;
+        }
+        /* subnormal x */
+        let mut e = -0x3ff;
+        while (i>>63) == 0 {
+            e -= 1;
+            i <<= 1;
+        }
+        return e;
+    }
+    if e == 0x7ff {
+        force_eval!(0.0/0.0);
+        if (i<<12) != 0 {
+            return FP_ILOGBNAN;
+        } else {
+            return isize::max_value();
+        }
+    }
+    return e - 0x3ff;
+}
diff --git a/src/math/ilogbf.rs b/src/math/ilogbf.rs
new file mode 100644
index 0000000..9ca1c36
--- /dev/null
+++ b/src/math/ilogbf.rs
@@ -0,0 +1,31 @@
+const FP_ILOGBNAN: isize = -1 - (((!0) >> 1));
+const FP_ILOGB0: isize = FP_ILOGBNAN;
+
+pub fn ilogbf(x: f32) -> isize {
+    let mut i = x.to_bits();
+    let e = ((i>>23) & 0xff) as isize;
+
+    if e == 0 {
+        i <<= 9;
+        if i == 0 {
+            force_eval!(0.0/0.0);
+            return FP_ILOGB0;
+        }
+        /* subnormal x */
+        let mut e = -0x7f;
+        while (i>>31) == 0 {
+            e -= 1;
+            i <<= 1;
+        }
+        return e;
+    }
+    if e == 0xff {
+        force_eval!(0.0/0.0);
+        if (i<<9) != 0 {
+            return FP_ILOGBNAN;
+        } else {
+            return isize::max_value();
+        }
+    }
+    return e - 0x7f;
+}
diff --git a/src/math/j0.rs b/src/math/j0.rs
new file mode 100644
index 0000000..02625b0
--- /dev/null
+++ b/src/math/j0.rs
@@ -0,0 +1,392 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/e_j0.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.
+ * ====================================================
+ */
+/* j0(x), y0(x)
+ * Bessel function of the first and second kinds of order zero.
+ * Method -- j0(x):
+ *      1. For tiny x, we use j0(x) = 1 - x^2/4 + x^4/64 - ...
+ *      2. Reduce x to |x| since j0(x)=j0(-x),  and
+ *         for x in (0,2)
+ *              j0(x) = 1-z/4+ z^2*R0/S0,  where z = x*x;
+ *         (precision:  |j0-1+z/4-z^2R0/S0 |<2**-63.67 )
+ *         for x in (2,inf)
+ *              j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)-q0(x)*sin(x0))
+ *         where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
+ *         as follow:
+ *              cos(x0) = cos(x)cos(pi/4)+sin(x)sin(pi/4)
+ *                      = 1/sqrt(2) * (cos(x) + sin(x))
+ *              sin(x0) = sin(x)cos(pi/4)-cos(x)sin(pi/4)
+ *                      = 1/sqrt(2) * (sin(x) - cos(x))
+ *         (To avoid cancellation, use
+ *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
+ *          to compute the worse one.)
+ *
+ *      3 Special cases
+ *              j0(nan)= nan
+ *              j0(0) = 1
+ *              j0(inf) = 0
+ *
+ * Method -- y0(x):
+ *      1. For x<2.
+ *         Since
+ *              y0(x) = 2/pi*(j0(x)*(ln(x/2)+Euler) + x^2/4 - ...)
+ *         therefore y0(x)-2/pi*j0(x)*ln(x) is an even function.
+ *         We use the following function to approximate y0,
+ *              y0(x) = U(z)/V(z) + (2/pi)*(j0(x)*ln(x)), z= x^2
+ *         where
+ *              U(z) = u00 + u01*z + ... + u06*z^6
+ *              V(z) = 1  + v01*z + ... + v04*z^4
+ *         with absolute approximation error bounded by 2**-72.
+ *         Note: For tiny x, U/V = u0 and j0(x)~1, hence
+ *              y0(tiny) = u0 + (2/pi)*ln(tiny), (choose tiny<2**-27)
+ *      2. For x>=2.
+ *              y0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x0)+q0(x)*sin(x0))
+ *         where x0 = x-pi/4. It is better to compute sin(x0),cos(x0)
+ *         by the method mentioned above.
+ *      3. Special cases: y0(0)=-inf, y0(x<0)=NaN, y0(inf)=0.
+ */
+
+use super::{cos, get_low_word, get_high_word, fabs, log, sin, sqrt};
+const INVSQRTPI: f64 = 5.64189583547756279280e-01; /* 0x3FE20DD7, 0x50429B6D */
+const TPI: f64       = 6.36619772367581382433e-01; /* 0x3FE45F30, 0x6DC9C883 */
+
+/* common method when |x|>=2 */
+fn common(ix: u32, x: f64, y0: bool) -> f64 {
+    let s: f64;
+    let mut c: f64;
+    let mut ss: f64;
+    let mut cc: f64;
+    let z: f64;
+
+    /*
+     * j0(x) = sqrt(2/(pi*x))*(p0(x)*cos(x-pi/4)-q0(x)*sin(x-pi/4))
+     * y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x-pi/4)+q0(x)*cos(x-pi/4))
+     *
+     * sin(x-pi/4) = (sin(x) - cos(x))/sqrt(2)
+     * cos(x-pi/4) = (sin(x) + cos(x))/sqrt(2)
+     * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
+     */
+    s = sin(x);
+    c = cos(x);
+    if y0 {
+        c = -c;
+    }
+    cc = s+c;
+    /* avoid overflow in 2*x, big ulp error when x>=0x1p1023 */
+    if ix < 0x7fe00000 {
+        ss = s-c;
+        z = -cos(2.0*x);
+        if s*c < 0.0 {
+            cc = z/ss;
+        } else {
+            ss = z/cc;
+        }
+        if ix < 0x48000000 {
+            if y0 {
+                ss = -ss;
+            }
+            cc = pzero(x)*cc-qzero(x)*ss;
+        }
+    }
+    return INVSQRTPI*cc/sqrt(x);
+}
+
+/* R0/S0 on [0, 2.00] */
+const R02: f64 =  1.56249999999999947958e-02; /* 0x3F8FFFFF, 0xFFFFFFFD */
+const R03: f64 = -1.89979294238854721751e-04; /* 0xBF28E6A5, 0xB61AC6E9 */
+const R04: f64 =  1.82954049532700665670e-06; /* 0x3EBEB1D1, 0x0C503919 */
+const R05: f64 = -4.61832688532103189199e-09; /* 0xBE33D5E7, 0x73D63FCE */
+const S01: f64 =  1.56191029464890010492e-02; /* 0x3F8FFCE8, 0x82C8C2A4 */
+const S02: f64 =  1.16926784663337450260e-04; /* 0x3F1EA6D2, 0xDD57DBF4 */
+const S03: f64 =  5.13546550207318111446e-07; /* 0x3EA13B54, 0xCE84D5A9 */
+const S04: f64 =  1.16614003333790000205e-09; /* 0x3E1408BC, 0xF4745D8F */
+
+pub fn j0(mut x: f64) -> f64
+{
+    let z: f64;
+    let r: f64;
+    let s: f64;
+    let mut ix: u32;
+
+    ix = get_high_word(x);
+    ix &= 0x7fffffff;
+
+    /* j0(+-inf)=0, j0(nan)=nan */
+    if ix >= 0x7ff00000 {
+        return 1.0/(x*x);
+    }
+    x = fabs(x);
+
+    if ix >= 0x40000000 {  /* |x| >= 2 */
+        /* large ulp error near zeros: 2.4, 5.52, 8.6537,.. */
+        return common(ix,x,false);
+    }
+
+    /* 1 - x*x/4 + x*x*R(x^2)/S(x^2) */
+    if ix >= 0x3f200000 {  /* |x| >= 2**-13 */
+        /* up to 4ulp error close to 2 */
+        z = x*x;
+        r = z*(R02+z*(R03+z*(R04+z*R05)));
+        s = 1.0+z*(S01+z*(S02+z*(S03+z*S04)));
+        return (1.0+x/2.0)*(1.0-x/2.0) + z*(r/s);
+    }
+
+    /* 1 - x*x/4 */
+    /* prevent underflow */
+    /* inexact should be raised when x!=0, this is not done correctly */
+    if ix >= 0x38000000 { /* |x| >= 2**-127 */
+        x = 0.25*x*x;
+    }
+    return 1.0 - x;
+}
+
+const U00: f64  = -7.38042951086872317523e-02; /* 0xBFB2E4D6, 0x99CBD01F */
+const U01: f64  =  1.76666452509181115538e-01; /* 0x3FC69D01, 0x9DE9E3FC */
+const U02: f64  = -1.38185671945596898896e-02; /* 0xBF8C4CE8, 0xB16CFA97 */
+const U03: f64  =  3.47453432093683650238e-04; /* 0x3F36C54D, 0x20B29B6B */
+const U04: f64  = -3.81407053724364161125e-06; /* 0xBECFFEA7, 0x73D25CAD */
+const U05: f64  =  1.95590137035022920206e-08; /* 0x3E550057, 0x3B4EABD4 */
+const U06: f64  = -3.98205194132103398453e-11; /* 0xBDC5E43D, 0x693FB3C8 */
+const V01: f64  =  1.27304834834123699328e-02; /* 0x3F8A1270, 0x91C9C71A */
+const V02: f64  =  7.60068627350353253702e-05; /* 0x3F13ECBB, 0xF578C6C1 */
+const V03: f64  =  2.59150851840457805467e-07; /* 0x3E91642D, 0x7FF202FD */
+const V04: f64  =  4.41110311332675467403e-10; /* 0x3DFE5018, 0x3BD6D9EF */
+
+pub fn y0(x: f64) -> f64
+{
+    let z: f64;
+    let u: f64;
+    let v: f64;
+    let ix: u32;
+    let lx: u32;
+
+    ix = get_high_word(x);
+    lx = get_low_word(x);
+
+    /* y0(nan)=nan, y0(<0)=nan, y0(0)=-inf, y0(inf)=0 */
+    if ((ix<<1) | lx) == 0 {
+        return -1.0/0.0;
+    }
+    if (ix>>31) != 0 {
+        return 0.0/0.0;
+    }
+    if ix >= 0x7ff00000 {
+        return 1.0/x;
+    }
+
+    if ix >= 0x40000000 {  /* x >= 2 */
+        /* large ulp errors near zeros: 3.958, 7.086,.. */
+        return common(ix,x,true);
+    }
+
+    /* U(x^2)/V(x^2) + (2/pi)*j0(x)*log(x) */
+    if ix >= 0x3e400000 {  /* x >= 2**-27 */
+        /* large ulp error near the first zero, x ~= 0.89 */
+        z = x*x;
+        u = U00+z*(U01+z*(U02+z*(U03+z*(U04+z*(U05+z*U06)))));
+        v = 1.0+z*(V01+z*(V02+z*(V03+z*V04)));
+        return u/v + TPI*(j0(x)*log(x));
+    }
+    return U00 + TPI*log(x);
+}
+
+/* The asymptotic expansions of pzero is
+ *      1 - 9/128 s^2 + 11025/98304 s^4 - ...,  where s = 1/x.
+ * For x >= 2, We approximate pzero by
+ *      pzero(x) = 1 + (R/S)
+ * where  R = pR0 + pR1*s^2 + pR2*s^4 + ... + pR5*s^10
+ *        S = 1 + pS0*s^2 + ... + pS4*s^10
+ * and
+ *      | pzero(x)-1-R/S | <= 2  ** ( -60.26)
+ */
+const PR8: [f64; 6] = [ /* for x in [inf, 8]=1/[0,0.125] */
+  0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
+ -7.03124999999900357484e-02, /* 0xBFB1FFFF, 0xFFFFFD32 */
+ -8.08167041275349795626e+00, /* 0xC02029D0, 0xB44FA779 */
+ -2.57063105679704847262e+02, /* 0xC0701102, 0x7B19E863 */
+ -2.48521641009428822144e+03, /* 0xC0A36A6E, 0xCD4DCAFC */
+ -5.25304380490729545272e+03, /* 0xC0B4850B, 0x36CC643D */
+];
+const PS8: [f64; 5] = [
+  1.16534364619668181717e+02, /* 0x405D2233, 0x07A96751 */
+  3.83374475364121826715e+03, /* 0x40ADF37D, 0x50596938 */
+  4.05978572648472545552e+04, /* 0x40E3D2BB, 0x6EB6B05F */
+  1.16752972564375915681e+05, /* 0x40FC810F, 0x8F9FA9BD */
+  4.76277284146730962675e+04, /* 0x40E74177, 0x4F2C49DC */
+];
+
+const PR5: [f64; 6] = [ /* for x in [8,4.5454]=1/[0.125,0.22001] */
+ -1.14125464691894502584e-11, /* 0xBDA918B1, 0x47E495CC */
+ -7.03124940873599280078e-02, /* 0xBFB1FFFF, 0xE69AFBC6 */
+ -4.15961064470587782438e+00, /* 0xC010A370, 0xF90C6BBF */
+ -6.76747652265167261021e+01, /* 0xC050EB2F, 0x5A7D1783 */
+ -3.31231299649172967747e+02, /* 0xC074B3B3, 0x6742CC63 */
+ -3.46433388365604912451e+02, /* 0xC075A6EF, 0x28A38BD7 */
+];
+const PS5: [f64; 5] = [
+  6.07539382692300335975e+01, /* 0x404E6081, 0x0C98C5DE */
+  1.05125230595704579173e+03, /* 0x40906D02, 0x5C7E2864 */
+  5.97897094333855784498e+03, /* 0x40B75AF8, 0x8FBE1D60 */
+  9.62544514357774460223e+03, /* 0x40C2CCB8, 0xFA76FA38 */
+  2.40605815922939109441e+03, /* 0x40A2CC1D, 0xC70BE864 */
+];
+
+const PR3: [f64; 6] = [/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
+ -2.54704601771951915620e-09, /* 0xBE25E103, 0x6FE1AA86 */
+ -7.03119616381481654654e-02, /* 0xBFB1FFF6, 0xF7C0E24B */
+ -2.40903221549529611423e+00, /* 0xC00345B2, 0xAEA48074 */
+ -2.19659774734883086467e+01, /* 0xC035F74A, 0x4CB94E14 */
+ -5.80791704701737572236e+01, /* 0xC04D0A22, 0x420A1A45 */
+ -3.14479470594888503854e+01, /* 0xC03F72AC, 0xA892D80F */
+];
+const PS3: [f64; 5] = [
+  3.58560338055209726349e+01, /* 0x4041ED92, 0x84077DD3 */
+  3.61513983050303863820e+02, /* 0x40769839, 0x464A7C0E */
+  1.19360783792111533330e+03, /* 0x4092A66E, 0x6D1061D6 */
+  1.12799679856907414432e+03, /* 0x40919FFC, 0xB8C39B7E */
+  1.73580930813335754692e+02, /* 0x4065B296, 0xFC379081 */
+];
+
+const PR2: [f64; 6] = [/* for x in [2.8570,2]=1/[0.3499,0.5] */
+ -8.87534333032526411254e-08, /* 0xBE77D316, 0xE927026D */
+ -7.03030995483624743247e-02, /* 0xBFB1FF62, 0x495E1E42 */
+ -1.45073846780952986357e+00, /* 0xBFF73639, 0x8A24A843 */
+ -7.63569613823527770791e+00, /* 0xC01E8AF3, 0xEDAFA7F3 */
+ -1.11931668860356747786e+01, /* 0xC02662E6, 0xC5246303 */
+ -3.23364579351335335033e+00, /* 0xC009DE81, 0xAF8FE70F */
+];
+const PS2: [f64; 5] = [
+  2.22202997532088808441e+01, /* 0x40363865, 0x908B5959 */
+  1.36206794218215208048e+02, /* 0x4061069E, 0x0EE8878F */
+  2.70470278658083486789e+02, /* 0x4070E786, 0x42EA079B */
+  1.53875394208320329881e+02, /* 0x40633C03, 0x3AB6FAFF */
+  1.46576176948256193810e+01, /* 0x402D50B3, 0x44391809 */
+];
+
+fn pzero(x: f64) -> f64
+{
+    let p: &[f64; 6];
+    let q: &[f64; 5];
+    let z: f64;
+    let r: f64;
+    let s: f64;
+    let mut ix: u32;
+
+    ix = get_high_word(x);
+    ix &= 0x7fffffff;
+    if      ix >= 0x40200000 {p = &PR8; q = &PS8;}
+    else if ix >= 0x40122E8B {p = &PR5; q = &PS5;}
+    else if ix >= 0x4006DB6D {p = &PR3; q = &PS3;}
+    else /*ix >= 0x40000000*/{p = &PR2; q = &PS2;}
+    z = 1.0/(x*x);
+    r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
+    s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
+    return 1.0 + r/s;
+}
+
+
+/* For x >= 8, the asymptotic expansions of qzero is
+ *      -1/8 s + 75/1024 s^3 - ..., where s = 1/x.
+ * We approximate pzero by
+ *      qzero(x) = s*(-1.25 + (R/S))
+ * where  R = qR0 + qR1*s^2 + qR2*s^4 + ... + qR5*s^10
+ *        S = 1 + qS0*s^2 + ... + qS5*s^12
+ * and
+ *      | qzero(x)/s +1.25-R/S | <= 2  ** ( -61.22)
+ */
+const QR8: [f64; 6] = [ /* for x in [inf, 8]=1/[0,0.125] */
+  0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
+  7.32421874999935051953e-02, /* 0x3FB2BFFF, 0xFFFFFE2C */
+  1.17682064682252693899e+01, /* 0x40278952, 0x5BB334D6 */
+  5.57673380256401856059e+02, /* 0x40816D63, 0x15301825 */
+  8.85919720756468632317e+03, /* 0x40C14D99, 0x3E18F46D */
+  3.70146267776887834771e+04, /* 0x40E212D4, 0x0E901566 */
+];
+const QS8: [f64; 6] = [
+  1.63776026895689824414e+02, /* 0x406478D5, 0x365B39BC */
+  8.09834494656449805916e+03, /* 0x40BFA258, 0x4E6B0563 */
+  1.42538291419120476348e+05, /* 0x41016652, 0x54D38C3F */
+  8.03309257119514397345e+05, /* 0x412883DA, 0x83A52B43 */
+  8.40501579819060512818e+05, /* 0x4129A66B, 0x28DE0B3D */
+ -3.43899293537866615225e+05, /* 0xC114FD6D, 0x2C9530C5 */
+];
+
+const QR5: [f64; 6] = [ /* for x in [8,4.5454]=1/[0.125,0.22001] */
+  1.84085963594515531381e-11, /* 0x3DB43D8F, 0x29CC8CD9 */
+  7.32421766612684765896e-02, /* 0x3FB2BFFF, 0xD172B04C */
+  5.83563508962056953777e+00, /* 0x401757B0, 0xB9953DD3 */
+  1.35111577286449829671e+02, /* 0x4060E392, 0x0A8788E9 */
+  1.02724376596164097464e+03, /* 0x40900CF9, 0x9DC8C481 */
+  1.98997785864605384631e+03, /* 0x409F17E9, 0x53C6E3A6 */
+];
+const QS5: [f64; 6] = [
+  8.27766102236537761883e+01, /* 0x4054B1B3, 0xFB5E1543 */
+  2.07781416421392987104e+03, /* 0x40A03BA0, 0xDA21C0CE */
+  1.88472887785718085070e+04, /* 0x40D267D2, 0x7B591E6D */
+  5.67511122894947329769e+04, /* 0x40EBB5E3, 0x97E02372 */
+  3.59767538425114471465e+04, /* 0x40E19118, 0x1F7A54A0 */
+ -5.35434275601944773371e+03, /* 0xC0B4EA57, 0xBEDBC609 */
+];
+
+const QR3: [f64; 6] = [/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
+  4.37741014089738620906e-09, /* 0x3E32CD03, 0x6ADECB82 */
+  7.32411180042911447163e-02, /* 0x3FB2BFEE, 0x0E8D0842 */
+  3.34423137516170720929e+00, /* 0x400AC0FC, 0x61149CF5 */
+  4.26218440745412650017e+01, /* 0x40454F98, 0x962DAEDD */
+  1.70808091340565596283e+02, /* 0x406559DB, 0xE25EFD1F */
+  1.66733948696651168575e+02, /* 0x4064D77C, 0x81FA21E0 */
+];
+const QS3: [f64; 6] = [
+  4.87588729724587182091e+01, /* 0x40486122, 0xBFE343A6 */
+  7.09689221056606015736e+02, /* 0x40862D83, 0x86544EB3 */
+  3.70414822620111362994e+03, /* 0x40ACF04B, 0xE44DFC63 */
+  6.46042516752568917582e+03, /* 0x40B93C6C, 0xD7C76A28 */
+  2.51633368920368957333e+03, /* 0x40A3A8AA, 0xD94FB1C0 */
+ -1.49247451836156386662e+02, /* 0xC062A7EB, 0x201CF40F */
+];
+
+const QR2: [f64; 6] = [/* for x in [2.8570,2]=1/[0.3499,0.5] */
+  1.50444444886983272379e-07, /* 0x3E84313B, 0x54F76BDB */
+  7.32234265963079278272e-02, /* 0x3FB2BEC5, 0x3E883E34 */
+  1.99819174093815998816e+00, /* 0x3FFFF897, 0xE727779C */
+  1.44956029347885735348e+01, /* 0x402CFDBF, 0xAAF96FE5 */
+  3.16662317504781540833e+01, /* 0x403FAA8E, 0x29FBDC4A */
+  1.62527075710929267416e+01, /* 0x403040B1, 0x71814BB4 */
+];
+const QS2: [f64; 6] = [
+  3.03655848355219184498e+01, /* 0x403E5D96, 0xF7C07AED */
+  2.69348118608049844624e+02, /* 0x4070D591, 0xE4D14B40 */
+  8.44783757595320139444e+02, /* 0x408A6645, 0x22B3BF22 */
+  8.82935845112488550512e+02, /* 0x408B977C, 0x9C5CC214 */
+  2.12666388511798828631e+02, /* 0x406A9553, 0x0E001365 */
+ -5.31095493882666946917e+00, /* 0xC0153E6A, 0xF8B32931 */
+];
+
+fn qzero(x: f64) -> f64
+{
+    let p: &[f64; 6];
+    let q: &[f64; 6];
+    let s: f64;
+    let r: f64;
+    let z: f64;
+    let mut ix: u32;
+
+    ix = get_high_word(x);
+    ix &= 0x7fffffff;
+    if      ix >= 0x40200000 {p = &QR8; q = &QS8;}
+    else if ix >= 0x40122E8B {p = &QR5; q = &QS5;}
+    else if ix >= 0x4006DB6D {p = &QR3; q = &QS3;}
+    else /*ix >= 0x40000000*/{p = &QR2; q = &QS2;}
+    z = 1.0/(x*x);
+    r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
+    s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
+    return (-0.125 + r/s)/x;
+} 
diff --git a/src/math/j0f.rs b/src/math/j0f.rs
new file mode 100644
index 0000000..e2faed0
--- /dev/null
+++ b/src/math/j0f.rs
@@ -0,0 +1,330 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/e_j0f.c */
+/*
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+/*
+ * ====================================================
+ * 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::{cosf, fabsf, logf, sinf, sqrtf};
+
+const INVSQRTPI: f32 = 5.6418961287e-01; /* 0x3f106ebb */
+const TPI: f32       = 6.3661974669e-01; /* 0x3f22f983 */
+
+fn common(ix: u32, x: f32, y0: bool) -> f32
+{
+    let z: f32;
+    let s: f32;
+    let mut c: f32;
+    let mut ss: f32;
+    let mut cc: f32;
+    /*
+     * j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
+     * y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
+     */
+    s = sinf(x);
+    c = cosf(x);
+    if y0 {
+        c = -c;
+    }
+    cc = s+c;
+    if ix < 0x7f000000 {
+        ss = s-c;
+        z = -cosf(2.0*x);
+        if s*c < 0.0 {
+            cc = z/ss;
+        } else {
+            ss = z/cc;
+        }
+        if ix < 0x58800000 {
+            if y0 {
+                ss = -ss;
+            }
+            cc = pzerof(x)*cc-qzerof(x)*ss;
+        }
+    }
+    return INVSQRTPI*cc/sqrtf(x);
+}
+
+/* R0/S0 on [0, 2.00] */
+const R02: f32 =  1.5625000000e-02; /* 0x3c800000 */
+const R03: f32 = -1.8997929874e-04; /* 0xb947352e */
+const R04: f32 =  1.8295404516e-06; /* 0x35f58e88 */
+const R05: f32 = -4.6183270541e-09; /* 0xb19eaf3c */
+const S01: f32 =  1.5619102865e-02; /* 0x3c7fe744 */
+const S02: f32 =  1.1692678527e-04; /* 0x38f53697 */
+const S03: f32 =  5.1354652442e-07; /* 0x3509daa6 */
+const S04: f32 =  1.1661400734e-09; /* 0x30a045e8 */
+
+pub fn j0f(mut x: f32) -> f32
+{
+    let z: f32;
+    let r: f32;
+    let s: f32;
+    let mut ix: u32;
+
+    ix = x.to_bits();
+    ix &= 0x7fffffff;
+    if ix >= 0x7f800000 {
+        return 1.0/(x*x);
+    }
+    x = fabsf(x);
+
+    if ix >= 0x40000000 {  /* |x| >= 2 */
+        /* large ulp error near zeros */
+        return common(ix, x, false);
+    }
+    if ix >= 0x3a000000 {  /* |x| >= 2**-11 */
+        /* up to 4ulp error near 2 */
+        z = x*x;
+        r = z*(R02+z*(R03+z*(R04+z*R05)));
+        s = 1.0+z*(S01+z*(S02+z*(S03+z*S04)));
+        return (1.0+x/2.0)*(1.0-x/2.0) + z*(r/s);
+    }
+    if ix >= 0x21800000 {  /* |x| >= 2**-60 */
+        x = 0.25*x*x;
+    }
+    return 1.0 - x;
+}
+
+const U00: f32  = -7.3804296553e-02; /* 0xbd9726b5 */
+const U01: f32  =  1.7666645348e-01; /* 0x3e34e80d */
+const U02: f32  = -1.3818567619e-02; /* 0xbc626746 */
+const U03: f32  =  3.4745343146e-04; /* 0x39b62a69 */
+const U04: f32  = -3.8140706238e-06; /* 0xb67ff53c */
+const U05: f32  =  1.9559013964e-08; /* 0x32a802ba */
+const U06: f32  = -3.9820518410e-11; /* 0xae2f21eb */
+const V01: f32  =  1.2730483897e-02; /* 0x3c509385 */
+const V02: f32  =  7.6006865129e-05; /* 0x389f65e0 */
+const V03: f32  =  2.5915085189e-07; /* 0x348b216c */
+const V04: f32  =  4.4111031494e-10; /* 0x2ff280c2 */
+
+pub fn y0f(x: f32) -> f32
+{
+    let z: f32;
+    let u: f32;
+    let v: f32;
+    let ix: u32;
+
+    ix = x.to_bits();
+    if (ix & 0x7fffffff) == 0 {
+        return -1.0/0.0;
+    }
+    if (ix>>31) !=0 {
+        return 0.0/0.0;
+    }
+    if ix >= 0x7f800000 {
+        return 1.0/x;
+    }
+    if ix >= 0x40000000 {  /* |x| >= 2.0 */
+        /* large ulp error near zeros */
+        return common(ix,x,true);
+    }
+    if ix >= 0x39000000 {  /* x >= 2**-13 */
+        /* large ulp error at x ~= 0.89 */
+        z = x*x;
+        u = U00+z*(U01+z*(U02+z*(U03+z*(U04+z*(U05+z*U06)))));
+        v = 1.0+z*(V01+z*(V02+z*(V03+z*V04)));
+        return u/v + TPI*(j0f(x)*logf(x));
+    }
+    return U00 + TPI*logf(x);
+}
+
+/* The asymptotic expansions of pzero is
+ *      1 - 9/128 s^2 + 11025/98304 s^4 - ...,  where s = 1/x.
+ * For x >= 2, We approximate pzero by
+ *      pzero(x) = 1 + (R/S)
+ * where  R = pR0 + pR1*s^2 + pR2*s^4 + ... + pR5*s^10
+ *        S = 1 + pS0*s^2 + ... + pS4*s^10
+ * and
+ *      | pzero(x)-1-R/S | <= 2  ** ( -60.26)
+ */
+const PR8: [f32; 6] = [ /* for x in [inf, 8]=1/[0,0.125] */
+  0.0000000000e+00, /* 0x00000000 */
+ -7.0312500000e-02, /* 0xbd900000 */
+ -8.0816707611e+00, /* 0xc1014e86 */
+ -2.5706311035e+02, /* 0xc3808814 */
+ -2.4852163086e+03, /* 0xc51b5376 */
+ -5.2530439453e+03, /* 0xc5a4285a */
+];
+const PS8: [f32; 5] = [
+  1.1653436279e+02, /* 0x42e91198 */
+  3.8337448730e+03, /* 0x456f9beb */
+  4.0597855469e+04, /* 0x471e95db */
+  1.1675296875e+05, /* 0x47e4087c */
+  4.7627726562e+04, /* 0x473a0bba */
+];
+const PR5: [f32; 6] = [ /* for x in [8,4.5454]=1/[0.125,0.22001] */
+ -1.1412546255e-11, /* 0xad48c58a */
+ -7.0312492549e-02, /* 0xbd8fffff */
+ -4.1596107483e+00, /* 0xc0851b88 */
+ -6.7674766541e+01, /* 0xc287597b */
+ -3.3123129272e+02, /* 0xc3a59d9b */
+ -3.4643338013e+02, /* 0xc3ad3779 */
+];
+const PS5: [f32; 5] = [
+  6.0753936768e+01, /* 0x42730408 */
+  1.0512523193e+03, /* 0x44836813 */
+  5.9789707031e+03, /* 0x45bad7c4 */
+  9.6254453125e+03, /* 0x461665c8 */
+  2.4060581055e+03, /* 0x451660ee */
+];
+
+const PR3: [f32; 6] = [/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
+ -2.5470459075e-09, /* 0xb12f081b */
+ -7.0311963558e-02, /* 0xbd8fffb8 */
+ -2.4090321064e+00, /* 0xc01a2d95 */
+ -2.1965976715e+01, /* 0xc1afba52 */
+ -5.8079170227e+01, /* 0xc2685112 */
+ -3.1447946548e+01, /* 0xc1fb9565 */
+];
+const PS3: [f32; 5] = [
+  3.5856033325e+01, /* 0x420f6c94 */
+  3.6151397705e+02, /* 0x43b4c1ca */
+  1.1936077881e+03, /* 0x44953373 */
+  1.1279968262e+03, /* 0x448cffe6 */
+  1.7358093262e+02, /* 0x432d94b8 */
+];
+
+const PR2: [f32; 6] = [/* for x in [2.8570,2]=1/[0.3499,0.5] */
+ -8.8753431271e-08, /* 0xb3be98b7 */
+ -7.0303097367e-02, /* 0xbd8ffb12 */
+ -1.4507384300e+00, /* 0xbfb9b1cc */
+ -7.6356959343e+00, /* 0xc0f4579f */
+ -1.1193166733e+01, /* 0xc1331736 */
+ -3.2336456776e+00, /* 0xc04ef40d */
+];
+const PS2: [f32; 5] = [
+  2.2220300674e+01, /* 0x41b1c32d */
+  1.3620678711e+02, /* 0x430834f0 */
+  2.7047027588e+02, /* 0x43873c32 */
+  1.5387539673e+02, /* 0x4319e01a */
+  1.4657617569e+01, /* 0x416a859a */
+];
+
+fn pzerof(x: f32) -> f32
+{
+    let p: &[f32; 6];
+    let q: &[f32; 5];
+    let z: f32;
+    let r: f32;
+    let s: f32;
+    let mut ix: u32;
+
+    ix = x.to_bits();
+    ix &= 0x7fffffff;
+    if      ix >= 0x41000000 {p = &PR8; q = &PS8;}
+    else if ix >= 0x409173eb {p = &PR5; q = &PS5;}
+    else if ix >= 0x4036d917 {p = &PR3; q = &PS3;}
+    else /*ix >= 0x40000000*/{p = &PR2; q = &PS2;}
+    z = 1.0/(x*x);
+    r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
+    s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
+    return 1.0 + r/s;
+}
+
+
+/* For x >= 8, the asymptotic expansions of qzero is
+ *      -1/8 s + 75/1024 s^3 - ..., where s = 1/x.
+ * We approximate pzero by
+ *      qzero(x) = s*(-1.25 + (R/S))
+ * where  R = qR0 + qR1*s^2 + qR2*s^4 + ... + qR5*s^10
+ *        S = 1 + qS0*s^2 + ... + qS5*s^12
+ * and
+ *      | qzero(x)/s +1.25-R/S | <= 2  ** ( -61.22)
+ */
+const QR8: [f32; 6] = [ /* for x in [inf, 8]=1/[0,0.125] */
+  0.0000000000e+00, /* 0x00000000 */
+  7.3242187500e-02, /* 0x3d960000 */
+  1.1768206596e+01, /* 0x413c4a93 */
+  5.5767340088e+02, /* 0x440b6b19 */
+  8.8591972656e+03, /* 0x460a6cca */
+  3.7014625000e+04, /* 0x471096a0 */
+];
+const QS8: [f32; 6] = [
+  1.6377603149e+02, /* 0x4323c6aa */
+  8.0983447266e+03, /* 0x45fd12c2 */
+  1.4253829688e+05, /* 0x480b3293 */
+  8.0330925000e+05, /* 0x49441ed4 */
+  8.4050156250e+05, /* 0x494d3359 */
+ -3.4389928125e+05, /* 0xc8a7eb69 */
+];
+
+const QR5: [f32; 6] = [ /* for x in [8,4.5454]=1/[0.125,0.22001] */
+  1.8408595828e-11, /* 0x2da1ec79 */
+  7.3242180049e-02, /* 0x3d95ffff */
+  5.8356351852e+00, /* 0x40babd86 */
+  1.3511157227e+02, /* 0x43071c90 */
+  1.0272437744e+03, /* 0x448067cd */
+  1.9899779053e+03, /* 0x44f8bf4b */
+];
+const QS5: [f32; 6] = [
+  8.2776611328e+01, /* 0x42a58da0 */
+  2.0778142090e+03, /* 0x4501dd07 */
+  1.8847289062e+04, /* 0x46933e94 */
+  5.6751113281e+04, /* 0x475daf1d */
+  3.5976753906e+04, /* 0x470c88c1 */
+ -5.3543427734e+03, /* 0xc5a752be */
+];
+
+const QR3: [f32; 6] = [/* for x in [4.547,2.8571]=1/[0.2199,0.35001] */
+  4.3774099900e-09, /* 0x3196681b */
+  7.3241114616e-02, /* 0x3d95ff70 */
+  3.3442313671e+00, /* 0x405607e3 */
+  4.2621845245e+01, /* 0x422a7cc5 */
+  1.7080809021e+02, /* 0x432acedf */
+  1.6673394775e+02, /* 0x4326bbe4 */
+];
+const QS3: [f32; 6] = [
+  4.8758872986e+01, /* 0x42430916 */
+  7.0968920898e+02, /* 0x44316c1c */
+  3.7041481934e+03, /* 0x4567825f */
+  6.4604252930e+03, /* 0x45c9e367 */
+  2.5163337402e+03, /* 0x451d4557 */
+ -1.4924745178e+02, /* 0xc3153f59 */
+];
+
+const QR2: [f32; 6] = [/* for x in [2.8570,2]=1/[0.3499,0.5] */
+  1.5044444979e-07, /* 0x342189db */
+  7.3223426938e-02, /* 0x3d95f62a */
+  1.9981917143e+00, /* 0x3fffc4bf */
+  1.4495602608e+01, /* 0x4167edfd */
+  3.1666231155e+01, /* 0x41fd5471 */
+  1.6252708435e+01, /* 0x4182058c */
+];
+const QS2: [f32; 6] = [
+  3.0365585327e+01, /* 0x41f2ecb8 */
+  2.6934811401e+02, /* 0x4386ac8f */
+  8.4478375244e+02, /* 0x44533229 */
+  8.8293585205e+02, /* 0x445cbbe5 */
+  2.1266638184e+02, /* 0x4354aa98 */
+ -5.3109550476e+00, /* 0xc0a9f358 */
+];
+
+fn qzerof(x: f32) -> f32
+{
+    let p: &[f32; 6];
+    let q: &[f32; 6];
+    let s: f32;
+    let r: f32;
+    let z: f32;
+    let mut ix: u32;
+
+    ix = x.to_bits();
+    ix &= 0x7fffffff;
+    if      ix >= 0x41000000 {p = &QR8; q = &QS8;}
+    else if ix >= 0x409173eb {p = &QR5; q = &QS5;}
+    else if ix >= 0x4036d917 {p = &QR3; q = &QS3;}
+    else /*ix >= 0x40000000*/{p = &QR2; q = &QS2;}
+    z = 1.0/(x*x);
+    r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
+    s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
+    return (-0.125 + r/s)/x;
+} 
diff --git a/src/math/j1.rs b/src/math/j1.rs
new file mode 100644
index 0000000..92289a6
--- /dev/null
+++ b/src/math/j1.rs
@@ -0,0 +1,387 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/e_j1.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.
+ * ====================================================
+ */
+/* j1(x), y1(x)
+ * Bessel function of the first and second kinds of order zero.
+ * Method -- j1(x):
+ *      1. For tiny x, we use j1(x) = x/2 - x^3/16 + x^5/384 - ...
+ *      2. Reduce x to |x| since j1(x)=-j1(-x),  and
+ *         for x in (0,2)
+ *              j1(x) = x/2 + x*z*R0/S0,  where z = x*x;
+ *         (precision:  |j1/x - 1/2 - R0/S0 |<2**-61.51 )
+ *         for x in (2,inf)
+ *              j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x1)-q1(x)*sin(x1))
+ *              y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
+ *         where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
+ *         as follow:
+ *              cos(x1) =  cos(x)cos(3pi/4)+sin(x)sin(3pi/4)
+ *                      =  1/sqrt(2) * (sin(x) - cos(x))
+ *              sin(x1) =  sin(x)cos(3pi/4)-cos(x)sin(3pi/4)
+ *                      = -1/sqrt(2) * (sin(x) + cos(x))
+ *         (To avoid cancellation, use
+ *              sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
+ *          to compute the worse one.)
+ *
+ *      3 Special cases
+ *              j1(nan)= nan
+ *              j1(0) = 0
+ *              j1(inf) = 0
+ *
+ * Method -- y1(x):
+ *      1. screen out x<=0 cases: y1(0)=-inf, y1(x<0)=NaN
+ *      2. For x<2.
+ *         Since
+ *              y1(x) = 2/pi*(j1(x)*(ln(x/2)+Euler)-1/x-x/2+5/64*x^3-...)
+ *         therefore y1(x)-2/pi*j1(x)*ln(x)-1/x is an odd function.
+ *         We use the following function to approximate y1,
+ *              y1(x) = x*U(z)/V(z) + (2/pi)*(j1(x)*ln(x)-1/x), z= x^2
+ *         where for x in [0,2] (abs err less than 2**-65.89)
+ *              U(z) = U0[0] + U0[1]*z + ... + U0[4]*z^4
+ *              V(z) = 1  + v0[0]*z + ... + v0[4]*z^5
+ *         Note: For tiny x, 1/x dominate y1 and hence
+ *              y1(tiny) = -2/pi/tiny, (choose tiny<2**-54)
+ *      3. For x>=2.
+ *              y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x1)+q1(x)*cos(x1))
+ *         where x1 = x-3*pi/4. It is better to compute sin(x1),cos(x1)
+ *         by method mentioned above.
+ */
+
+use super::{cos, get_high_word, get_low_word, fabs, log, sin, sqrt};
+
+const INVSQRTPI: f64 = 5.64189583547756279280e-01; /* 0x3FE20DD7, 0x50429B6D */
+const TPI: f64       = 6.36619772367581382433e-01; /* 0x3FE45F30, 0x6DC9C883 */
+
+fn common(ix: u32, x: f64, y1: bool, sign: bool) -> f64
+{
+    let z: f64;
+    let mut s: f64;
+    let c: f64;
+    let mut ss: f64;
+    let mut cc: f64;
+
+    /*
+     * j1(x) = sqrt(2/(pi*x))*(p1(x)*cos(x-3pi/4)-q1(x)*sin(x-3pi/4))
+     * y1(x) = sqrt(2/(pi*x))*(p1(x)*sin(x-3pi/4)+q1(x)*cos(x-3pi/4))
+     *
+     * sin(x-3pi/4) = -(sin(x) + cos(x))/sqrt(2)
+     * cos(x-3pi/4) = (sin(x) - cos(x))/sqrt(2)
+     * sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
+     */
+    s = sin(x);
+    if y1 {
+        s = -s;
+    }
+    c = cos(x);
+    cc = s-c;
+    if ix < 0x7fe00000 {
+        /* avoid overflow in 2*x */
+        ss = -s-c;
+        z = cos(2.0*x);
+        if s*c > 0.0 {
+            cc = z/ss;
+        } else {
+            ss = z/cc;
+        }
+        if ix < 0x48000000 {
+            if y1 {
+                ss = -ss;
+            }
+            cc = pone(x)*cc-qone(x)*ss;
+        }
+    }
+    if sign {
+        cc = -cc;
+    }
+    return INVSQRTPI*cc/sqrt(x);
+}
+
+/* R0/S0 on [0,2] */
+const R00: f64 = -6.25000000000000000000e-02; /* 0xBFB00000, 0x00000000 */
+const R01: f64 =  1.40705666955189706048e-03; /* 0x3F570D9F, 0x98472C61 */
+const R02: f64 = -1.59955631084035597520e-05; /* 0xBEF0C5C6, 0xBA169668 */
+const R03: f64 =  4.96727999609584448412e-08; /* 0x3E6AAAFA, 0x46CA0BD9 */
+const S01: f64 =  1.91537599538363460805e-02; /* 0x3F939D0B, 0x12637E53 */
+const S02: f64 =  1.85946785588630915560e-04; /* 0x3F285F56, 0xB9CDF664 */
+const S03: f64 =  1.17718464042623683263e-06; /* 0x3EB3BFF8, 0x333F8498 */
+const S04: f64 =  5.04636257076217042715e-09; /* 0x3E35AC88, 0xC97DFF2C */
+const S05: f64 =  1.23542274426137913908e-11; /* 0x3DAB2ACF, 0xCFB97ED8 */
+
+pub fn j1(x: f64) -> f64
+{
+    let mut z: f64;
+    let r: f64;
+    let s: f64;
+    let mut ix: u32;
+    let sign: bool;
+
+    ix = get_high_word(x);
+    sign = (ix>>31) != 0;
+    ix &= 0x7fffffff;
+    if ix >= 0x7ff00000 {
+        return 1.0/(x*x);
+    }
+    if ix >= 0x40000000 { /* |x| >= 2 */
+        return common(ix, fabs(x), false, sign);
+    }
+    if ix >= 0x38000000 {  /* |x| >= 2**-127 */
+        z = x*x;
+        r = z*(R00+z*(R01+z*(R02+z*R03)));
+        s = 1.0+z*(S01+z*(S02+z*(S03+z*(S04+z*S05))));
+        z = r/s;
+    } else {
+        /* avoid underflow, raise inexact if x!=0 */
+        z = x;
+    }
+    return (0.5 + z)*x;
+}
+
+const U0: [f64; 5] = [
+ -1.96057090646238940668e-01, /* 0xBFC91866, 0x143CBC8A */
+  5.04438716639811282616e-02, /* 0x3FA9D3C7, 0x76292CD1 */
+ -1.91256895875763547298e-03, /* 0xBF5F55E5, 0x4844F50F */
+  2.35252600561610495928e-05, /* 0x3EF8AB03, 0x8FA6B88E */
+ -9.19099158039878874504e-08, /* 0xBE78AC00, 0x569105B8 */
+];
+const V0: [f64; 5] = [
+  1.99167318236649903973e-02, /* 0x3F94650D, 0x3F4DA9F0 */
+  2.02552581025135171496e-04, /* 0x3F2A8C89, 0x6C257764 */
+  1.35608801097516229404e-06, /* 0x3EB6C05A, 0x894E8CA6 */
+  6.22741452364621501295e-09, /* 0x3E3ABF1D, 0x5BA69A86 */
+  1.66559246207992079114e-11, /* 0x3DB25039, 0xDACA772A */
+];
+
+pub fn y1(x: f64) -> f64
+{
+    let z: f64;
+    let u: f64;
+    let v: f64;
+    let ix: u32;
+    let lx: u32;
+
+    ix = get_high_word(x);
+    lx = get_low_word(x);
+
+    /* y1(nan)=nan, y1(<0)=nan, y1(0)=-inf, y1(inf)=0 */
+    if (ix<<1 | lx) == 0 {
+        return -1.0/0.0;
+    }
+    if (ix>>31) != 0 {
+        return 0.0/0.0;
+    }
+    if ix >= 0x7ff00000 {
+        return 1.0/x;
+    }
+
+    if ix >= 0x40000000 { /* x >= 2 */
+        return common(ix, x, true, false);
+    }
+    if ix < 0x3c900000 { /* x < 2**-54 */
+        return -TPI/x;
+    }
+    z = x*x;
+    u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
+    v = 1.0+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
+    return x*(u/v) + TPI*(j1(x)*log(x)-1.0/x);
+}
+
+/* For x >= 8, the asymptotic expansions of pone is
+ *      1 + 15/128 s^2 - 4725/2^15 s^4 - ...,   where s = 1/x.
+ * We approximate pone by
+ *      pone(x) = 1 + (R/S)
+ * where  R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10
+ *        S = 1 + ps0*s^2 + ... + ps4*s^10
+ * and
+ *      | pone(x)-1-R/S | <= 2  ** ( -60.06)
+ */
+
+const PR8: [f64; 6] = [ /* for x in [inf, 8]=1/[0,0.125] */
+  0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
+  1.17187499999988647970e-01, /* 0x3FBDFFFF, 0xFFFFFCCE */
+  1.32394806593073575129e+01, /* 0x402A7A9D, 0x357F7FCE */
+  4.12051854307378562225e+02, /* 0x4079C0D4, 0x652EA590 */
+  3.87474538913960532227e+03, /* 0x40AE457D, 0xA3A532CC */
+  7.91447954031891731574e+03, /* 0x40BEEA7A, 0xC32782DD */
+];
+const PS8: [f64; 5] = [
+  1.14207370375678408436e+02, /* 0x405C8D45, 0x8E656CAC */
+  3.65093083420853463394e+03, /* 0x40AC85DC, 0x964D274F */
+  3.69562060269033463555e+04, /* 0x40E20B86, 0x97C5BB7F */
+  9.76027935934950801311e+04, /* 0x40F7D42C, 0xB28F17BB */
+  3.08042720627888811578e+04, /* 0x40DE1511, 0x697A0B2D */
+];
+
+const PR5: [f64; 6] = [ /* for x in [8,4.5454]=1/[0.125,0.22001] */
+  1.31990519556243522749e-11, /* 0x3DAD0667, 0xDAE1CA7D */
+  1.17187493190614097638e-01, /* 0x3FBDFFFF, 0xE2C10043 */
+  6.80275127868432871736e+00, /* 0x401B3604, 0x6E6315E3 */
+  1.08308182990189109773e+02, /* 0x405B13B9, 0x452602ED */
+  5.17636139533199752805e+02, /* 0x40802D16, 0xD052D649 */
+  5.28715201363337541807e+02, /* 0x408085B8, 0xBB7E0CB7 */
+];
+const PS5: [f64; 5] = [
+  5.92805987221131331921e+01, /* 0x404DA3EA, 0xA8AF633D */
+  9.91401418733614377743e+02, /* 0x408EFB36, 0x1B066701 */
+  5.35326695291487976647e+03, /* 0x40B4E944, 0x5706B6FB */
+  7.84469031749551231769e+03, /* 0x40BEA4B0, 0xB8A5BB15 */
+  1.50404688810361062679e+03, /* 0x40978030, 0x036F5E51 */
+];
+
+const PR3: [f64; 6] = [
+  3.02503916137373618024e-09, /* 0x3E29FC21, 0xA7AD9EDD */
+  1.17186865567253592491e-01, /* 0x3FBDFFF5, 0x5B21D17B */
+  3.93297750033315640650e+00, /* 0x400F76BC, 0xE85EAD8A */
+  3.51194035591636932736e+01, /* 0x40418F48, 0x9DA6D129 */
+  9.10550110750781271918e+01, /* 0x4056C385, 0x4D2C1837 */
+  4.85590685197364919645e+01, /* 0x4048478F, 0x8EA83EE5 */
+];
+const PS3: [f64; 5] = [
+  3.47913095001251519989e+01, /* 0x40416549, 0xA134069C */
+  3.36762458747825746741e+02, /* 0x40750C33, 0x07F1A75F */
+  1.04687139975775130551e+03, /* 0x40905B7C, 0x5037D523 */
+  8.90811346398256432622e+02, /* 0x408BD67D, 0xA32E31E9 */
+  1.03787932439639277504e+02, /* 0x4059F26D, 0x7C2EED53 */
+];
+
+const PR2: [f64; 6] = [/* for x in [2.8570,2]=1/[0.3499,0.5] */
+  1.07710830106873743082e-07, /* 0x3E7CE9D4, 0xF65544F4 */
+  1.17176219462683348094e-01, /* 0x3FBDFF42, 0xBE760D83 */
+  2.36851496667608785174e+00, /* 0x4002F2B7, 0xF98FAEC0 */
+  1.22426109148261232917e+01, /* 0x40287C37, 0x7F71A964 */
+  1.76939711271687727390e+01, /* 0x4031B1A8, 0x177F8EE2 */
+  5.07352312588818499250e+00, /* 0x40144B49, 0xA574C1FE */
+];
+const PS2: [f64; 5] = [
+  2.14364859363821409488e+01, /* 0x40356FBD, 0x8AD5ECDC */
+  1.25290227168402751090e+02, /* 0x405F5293, 0x14F92CD5 */
+  2.32276469057162813669e+02, /* 0x406D08D8, 0xD5A2DBD9 */
+  1.17679373287147100768e+02, /* 0x405D6B7A, 0xDA1884A9 */
+  8.36463893371618283368e+00, /* 0x4020BAB1, 0xF44E5192 */
+];
+
+fn pone(x: f64) -> f64
+{
+    let p: &[f64; 6];
+    let q: &[f64; 5];
+    let z: f64;
+    let r: f64;
+    let s: f64;
+    let mut ix: u32;
+
+    ix = get_high_word(x);
+    ix &= 0x7fffffff;
+    if      ix >= 0x40200000 {p = &PR8; q = &PS8;}
+    else if ix >= 0x40122E8B {p = &PR5; q = &PS5;}
+    else if ix >= 0x4006DB6D {p = &PR3; q = &PS3;}
+    else /*ix >= 0x40000000*/{p = &PR2; q = &PS2;}
+    z = 1.0/(x*x);
+    r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
+    s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
+    return 1.0+ r/s;
+}
+
+/* For x >= 8, the asymptotic expansions of qone is
+ *      3/8 s - 105/1024 s^3 - ..., where s = 1/x.
+ * We approximate pone by
+ *      qone(x) = s*(0.375 + (R/S))
+ * where  R = qr1*s^2 + qr2*s^4 + ... + qr5*s^10
+ *        S = 1 + qs1*s^2 + ... + qs6*s^12
+ * and
+ *      | qone(x)/s -0.375-R/S | <= 2  ** ( -61.13)
+ */
+
+const QR8: [f64; 6] = [ /* for x in [inf, 8]=1/[0,0.125] */
+  0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
+ -1.02539062499992714161e-01, /* 0xBFBA3FFF, 0xFFFFFDF3 */
+ -1.62717534544589987888e+01, /* 0xC0304591, 0xA26779F7 */
+ -7.59601722513950107896e+02, /* 0xC087BCD0, 0x53E4B576 */
+ -1.18498066702429587167e+04, /* 0xC0C724E7, 0x40F87415 */
+ -4.84385124285750353010e+04, /* 0xC0E7A6D0, 0x65D09C6A */
+];
+const QS8: [f64; 6] = [
+  1.61395369700722909556e+02, /* 0x40642CA6, 0xDE5BCDE5 */
+  7.82538599923348465381e+03, /* 0x40BE9162, 0xD0D88419 */
+  1.33875336287249578163e+05, /* 0x4100579A, 0xB0B75E98 */
+  7.19657723683240939863e+05, /* 0x4125F653, 0x72869C19 */
+  6.66601232617776375264e+05, /* 0x412457D2, 0x7719AD5C */
+ -2.94490264303834643215e+05, /* 0xC111F969, 0x0EA5AA18 */
+];
+
+const QR5: [f64; 6] = [ /* for x in [8,4.5454]=1/[0.125,0.22001] */
+ -2.08979931141764104297e-11, /* 0xBDB6FA43, 0x1AA1A098 */
+ -1.02539050241375426231e-01, /* 0xBFBA3FFF, 0xCB597FEF */
+ -8.05644828123936029840e+00, /* 0xC0201CE6, 0xCA03AD4B */
+ -1.83669607474888380239e+02, /* 0xC066F56D, 0x6CA7B9B0 */
+ -1.37319376065508163265e+03, /* 0xC09574C6, 0x6931734F */
+ -2.61244440453215656817e+03, /* 0xC0A468E3, 0x88FDA79D */
+];
+const QS5: [f64; 6] = [
+  8.12765501384335777857e+01, /* 0x405451B2, 0xFF5A11B2 */
+  1.99179873460485964642e+03, /* 0x409F1F31, 0xE77BF839 */
+  1.74684851924908907677e+04, /* 0x40D10F1F, 0x0D64CE29 */
+  4.98514270910352279316e+04, /* 0x40E8576D, 0xAABAD197 */
+  2.79480751638918118260e+04, /* 0x40DB4B04, 0xCF7C364B */
+ -4.71918354795128470869e+03, /* 0xC0B26F2E, 0xFCFFA004 */
+];
+
+const QR3: [f64; 6] = [
+ -5.07831226461766561369e-09, /* 0xBE35CFA9, 0xD38FC84F */
+ -1.02537829820837089745e-01, /* 0xBFBA3FEB, 0x51AEED54 */
+ -4.61011581139473403113e+00, /* 0xC01270C2, 0x3302D9FF */
+ -5.78472216562783643212e+01, /* 0xC04CEC71, 0xC25D16DA */
+ -2.28244540737631695038e+02, /* 0xC06C87D3, 0x4718D55F */
+ -2.19210128478909325622e+02, /* 0xC06B66B9, 0x5F5C1BF6 */
+];
+const QS3: [f64; 6] = [
+  4.76651550323729509273e+01, /* 0x4047D523, 0xCCD367E4 */
+  6.73865112676699709482e+02, /* 0x40850EEB, 0xC031EE3E */
+  3.38015286679526343505e+03, /* 0x40AA684E, 0x448E7C9A */
+  5.54772909720722782367e+03, /* 0x40B5ABBA, 0xA61D54A6 */
+  1.90311919338810798763e+03, /* 0x409DBC7A, 0x0DD4DF4B */
+ -1.35201191444307340817e+02, /* 0xC060E670, 0x290A311F */
+];
+
+const QR2: [f64; 6] = [/* for x in [2.8570,2]=1/[0.3499,0.5] */
+ -1.78381727510958865572e-07, /* 0xBE87F126, 0x44C626D2 */
+ -1.02517042607985553460e-01, /* 0xBFBA3E8E, 0x9148B010 */
+ -2.75220568278187460720e+00, /* 0xC0060484, 0x69BB4EDA */
+ -1.96636162643703720221e+01, /* 0xC033A9E2, 0xC168907F */
+ -4.23253133372830490089e+01, /* 0xC04529A3, 0xDE104AAA */
+ -2.13719211703704061733e+01, /* 0xC0355F36, 0x39CF6E52 */
+];
+const QS2: [f64; 6] = [
+  2.95333629060523854548e+01, /* 0x403D888A, 0x78AE64FF */
+  2.52981549982190529136e+02, /* 0x406F9F68, 0xDB821CBA */
+  7.57502834868645436472e+02, /* 0x4087AC05, 0xCE49A0F7 */
+  7.39393205320467245656e+02, /* 0x40871B25, 0x48D4C029 */
+  1.55949003336666123687e+02, /* 0x40637E5E, 0x3C3ED8D4 */
+ -4.95949898822628210127e+00, /* 0xC013D686, 0xE71BE86B */
+];
+
+fn qone(x: f64) -> f64
+{
+    let p: &[f64; 6];
+    let q: &[f64; 6];
+    let s: f64;
+    let r: f64;
+    let z: f64;
+    let mut ix: u32;
+
+    ix = get_high_word(x);
+    ix &= 0x7fffffff;
+    if      ix >= 0x40200000 {p = &QR8; q = &QS8;}
+    else if ix >= 0x40122E8B {p = &QR5; q = &QS5;}
+    else if ix >= 0x4006DB6D {p = &QR3; q = &QS3;}
+    else /*ix >= 0x40000000*/{p = &QR2; q = &QS2;}
+    z = 1.0/(x*x);
+    r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
+    s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
+    return (0.375 + r/s)/x;
+}
diff --git a/src/math/j1f.rs b/src/math/j1f.rs
new file mode 100644
index 0000000..7cf9c45
--- /dev/null
+++ b/src/math/j1f.rs
@@ -0,0 +1,331 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/e_j1f.c */
+/*
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+/*
+ * ====================================================
+ * 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::{cosf, fabsf, logf, sinf, sqrtf};
+
+const INVSQRTPI: f32 = 5.6418961287e-01; /* 0x3f106ebb */
+const TPI: f32       = 6.3661974669e-01; /* 0x3f22f983 */
+
+fn common(ix: u32, x: f32, y1: bool, sign: bool) -> f32
+{
+    let z: f64;
+    let mut s: f64;
+    let c: f64;
+    let mut ss: f64;
+    let mut cc: f64;
+
+    s = sinf(x) as f64;
+    if y1 {
+        s = -s;
+    }
+    c = cosf(x) as f64;
+    cc = s-c;
+    if ix < 0x7f000000 {
+        ss = -s-c;
+        z = cosf(2.0*x) as f64;
+        if s*c > 0.0 {
+            cc = z/ss;
+        } else {
+            ss = z/cc;
+        }
+        if ix < 0x58800000 {
+            if y1 {
+                ss = -ss;
+            }
+            cc = (ponef(x) as f64)*cc-(qonef(x) as f64)*ss;
+        }
+    }
+    if sign {
+        cc = -cc;
+    }
+    return INVSQRTPI*(cc as f32)/sqrtf(x);
+}
+
+/* R0/S0 on [0,2] */
+const R00: f32 = -6.2500000000e-02; /* 0xbd800000 */
+const R01: f32 =  1.4070566976e-03; /* 0x3ab86cfd */
+const R02: f32 = -1.5995563444e-05; /* 0xb7862e36 */
+const R03: f32 =  4.9672799207e-08; /* 0x335557d2 */
+const S01: f32 =  1.9153760746e-02; /* 0x3c9ce859 */
+const S02: f32 =  1.8594678841e-04; /* 0x3942fab6 */
+const S03: f32 =  1.1771846857e-06; /* 0x359dffc2 */
+const S04: f32 =  5.0463624390e-09; /* 0x31ad6446 */
+const S05: f32 =  1.2354227016e-11; /* 0x2d59567e */
+
+pub fn j1f(x: f32) -> f32
+{
+    let mut z: f32;
+    let r: f32;
+    let s: f32;
+    let mut ix: u32;
+    let sign: bool;
+
+    ix = x.to_bits();
+    sign = (ix>>31) != 0;
+    ix &= 0x7fffffff;
+    if ix >= 0x7f800000 {
+        return 1.0/(x*x);
+    }
+    if ix >= 0x40000000 {  /* |x| >= 2 */
+        return common(ix, fabsf(x), false, sign);
+    }
+    if ix >= 0x39000000 {  /* |x| >= 2**-13 */
+        z = x*x;
+        r = z*(R00+z*(R01+z*(R02+z*R03)));
+        s = 1.0+z*(S01+z*(S02+z*(S03+z*(S04+z*S05))));
+        z = 0.5 + r/s;
+    } else {
+        z = 0.5;
+    }
+    return z*x;
+}
+
+const U0: [f32; 5] = [
+ -1.9605709612e-01, /* 0xbe48c331 */
+  5.0443872809e-02, /* 0x3d4e9e3c */
+ -1.9125689287e-03, /* 0xbafaaf2a */
+  2.3525259166e-05, /* 0x37c5581c */
+ -9.1909917899e-08, /* 0xb3c56003 */
+];
+const V0: [f32; 5] = [
+  1.9916731864e-02, /* 0x3ca3286a */
+  2.0255257550e-04, /* 0x3954644b */
+  1.3560879779e-06, /* 0x35b602d4 */
+  6.2274145840e-09, /* 0x31d5f8eb */
+  1.6655924903e-11, /* 0x2d9281cf */
+];
+
+pub fn y1f(x: f32) -> f32
+{
+    let z: f32;
+    let u: f32;
+    let v: f32;
+    let ix: u32;
+
+    ix = x.to_bits();
+    if (ix & 0x7fffffff) == 0 {
+        return -1.0/0.0;
+    }
+    if (ix>>31) != 0{
+        return 0.0/0.0;
+    }
+    if ix >= 0x7f800000 {
+        return 1.0/x;
+    }
+    if ix >= 0x40000000 {  /* |x| >= 2.0 */
+        return common(ix,x,true,false);
+    }
+    if ix < 0x33000000 {  /* x < 2**-25 */
+        return -TPI/x;
+    }
+    z = x*x;
+    u = U0[0]+z*(U0[1]+z*(U0[2]+z*(U0[3]+z*U0[4])));
+    v = 1.0+z*(V0[0]+z*(V0[1]+z*(V0[2]+z*(V0[3]+z*V0[4]))));
+    return x*(u/v) + TPI*(j1f(x)*logf(x)-1.0/x);
+}
+
+/* For x >= 8, the asymptotic expansions of pone is
+ *      1 + 15/128 s^2 - 4725/2^15 s^4 - ...,   where s = 1/x.
+ * We approximate pone by
+ *      pone(x) = 1 + (R/S)
+ * where  R = pr0 + pr1*s^2 + pr2*s^4 + ... + pr5*s^10
+ *        S = 1 + ps0*s^2 + ... + ps4*s^10
+ * and
+ *      | pone(x)-1-R/S | <= 2  ** ( -60.06)
+ */
+
+const PR8: [f32; 6] = [ /* for x in [inf, 8]=1/[0,0.125] */
+  0.0000000000e+00, /* 0x00000000 */
+  1.1718750000e-01, /* 0x3df00000 */
+  1.3239480972e+01, /* 0x4153d4ea */
+  4.1205184937e+02, /* 0x43ce06a3 */
+  3.8747453613e+03, /* 0x45722bed */
+  7.9144794922e+03, /* 0x45f753d6 */
+];
+const PS8: [f32; 5] = [
+  1.1420736694e+02, /* 0x42e46a2c */
+  3.6509309082e+03, /* 0x45642ee5 */
+  3.6956207031e+04, /* 0x47105c35 */
+  9.7602796875e+04, /* 0x47bea166 */
+  3.0804271484e+04, /* 0x46f0a88b */
+];
+
+const PR5: [f32; 6] = [ /* for x in [8,4.5454]=1/[0.125,0.22001] */
+  1.3199052094e-11, /* 0x2d68333f */
+  1.1718749255e-01, /* 0x3defffff */
+  6.8027510643e+00, /* 0x40d9b023 */
+  1.0830818176e+02, /* 0x42d89dca */
+  5.1763616943e+02, /* 0x440168b7 */
+  5.2871520996e+02, /* 0x44042dc6 */
+];
+const PS5: [f32; 5] = [
+  5.9280597687e+01, /* 0x426d1f55 */
+  9.9140142822e+02, /* 0x4477d9b1 */
+  5.3532670898e+03, /* 0x45a74a23 */
+  7.8446904297e+03, /* 0x45f52586 */
+  1.5040468750e+03, /* 0x44bc0180 */
+];
+
+const PR3: [f32; 6] = [
+  3.0250391081e-09, /* 0x314fe10d */
+  1.1718686670e-01, /* 0x3defffab */
+  3.9329774380e+00, /* 0x407bb5e7 */
+  3.5119403839e+01, /* 0x420c7a45 */
+  9.1055007935e+01, /* 0x42b61c2a */
+  4.8559066772e+01, /* 0x42423c7c */
+];
+const PS3: [f32; 5] = [
+  3.4791309357e+01, /* 0x420b2a4d */
+  3.3676245117e+02, /* 0x43a86198 */
+  1.0468714600e+03, /* 0x4482dbe3 */
+  8.9081134033e+02, /* 0x445eb3ed */
+  1.0378793335e+02, /* 0x42cf936c */
+];
+
+const PR2: [f32; 6] = [/* for x in [2.8570,2]=1/[0.3499,0.5] */
+  1.0771083225e-07, /* 0x33e74ea8 */
+  1.1717621982e-01, /* 0x3deffa16 */
+  2.3685150146e+00, /* 0x401795c0 */
+  1.2242610931e+01, /* 0x4143e1bc */
+  1.7693971634e+01, /* 0x418d8d41 */
+  5.0735230446e+00, /* 0x40a25a4d */
+];
+const PS2: [f32; 5] = [
+  2.1436485291e+01, /* 0x41ab7dec */
+  1.2529022980e+02, /* 0x42fa9499 */
+  2.3227647400e+02, /* 0x436846c7 */
+  1.1767937469e+02, /* 0x42eb5bd7 */
+  8.3646392822e+00, /* 0x4105d590 */
+];
+
+fn ponef(x: f32) -> f32
+{
+    let p: &[f32; 6];
+    let q: &[f32; 5];
+    let z: f32;
+    let r: f32;
+    let s: f32;
+    let mut ix: u32;
+
+    ix = x.to_bits();
+    ix &= 0x7fffffff;
+    if      ix >= 0x41000000 {p = &PR8; q = &PS8;}
+    else if ix >= 0x409173eb {p = &PR5; q = &PS5;}
+    else if ix >= 0x4036d917 {p = &PR3; q = &PS3;}
+    else /*ix >= 0x40000000*/{p = &PR2; q = &PS2;}
+    z = 1.0/(x*x);
+    r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
+    s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
+    return 1.0 + r/s;
+}
+
+/* For x >= 8, the asymptotic expansions of qone is
+ *      3/8 s - 105/1024 s^3 - ..., where s = 1/x.
+ * We approximate pone by
+ *      qone(x) = s*(0.375 + (R/S))
+ * where  R = qr1*s^2 + qr2*s^4 + ... + qr5*s^10
+ *        S = 1 + qs1*s^2 + ... + qs6*s^12
+ * and
+ *      | qone(x)/s -0.375-R/S | <= 2  ** ( -61.13)
+ */
+
+const QR8: [f32; 6] = [ /* for x in [inf, 8]=1/[0,0.125] */
+  0.0000000000e+00, /* 0x00000000 */
+ -1.0253906250e-01, /* 0xbdd20000 */
+ -1.6271753311e+01, /* 0xc1822c8d */
+ -7.5960174561e+02, /* 0xc43de683 */
+ -1.1849806641e+04, /* 0xc639273a */
+ -4.8438511719e+04, /* 0xc73d3683 */
+];
+const QS8: [f32; 6] = [
+  1.6139537048e+02, /* 0x43216537 */
+  7.8253862305e+03, /* 0x45f48b17 */
+  1.3387534375e+05, /* 0x4802bcd6 */
+  7.1965775000e+05, /* 0x492fb29c */
+  6.6660125000e+05, /* 0x4922be94 */
+ -2.9449025000e+05, /* 0xc88fcb48 */
+];
+
+const QR5: [f32; 6] = [ /* for x in [8,4.5454]=1/[0.125,0.22001] */
+ -2.0897993405e-11, /* 0xadb7d219 */
+ -1.0253904760e-01, /* 0xbdd1fffe */
+ -8.0564479828e+00, /* 0xc100e736 */
+ -1.8366960144e+02, /* 0xc337ab6b */
+ -1.3731937256e+03, /* 0xc4aba633 */
+ -2.6124443359e+03, /* 0xc523471c */
+];
+const QS5: [f32; 6] = [
+  8.1276550293e+01, /* 0x42a28d98 */
+  1.9917987061e+03, /* 0x44f8f98f */
+  1.7468484375e+04, /* 0x468878f8 */
+  4.9851425781e+04, /* 0x4742bb6d */
+  2.7948074219e+04, /* 0x46da5826 */
+ -4.7191835938e+03, /* 0xc5937978 */
+];
+
+const QR3: [f32; 6] = [
+ -5.0783124372e-09, /* 0xb1ae7d4f */
+ -1.0253783315e-01, /* 0xbdd1ff5b */
+ -4.6101160049e+00, /* 0xc0938612 */
+ -5.7847221375e+01, /* 0xc267638e */
+ -2.2824453735e+02, /* 0xc3643e9a */
+ -2.1921012878e+02, /* 0xc35b35cb */
+];
+const QS3: [f32; 6] = [
+  4.7665153503e+01, /* 0x423ea91e */
+  6.7386511230e+02, /* 0x4428775e */
+  3.3801528320e+03, /* 0x45534272 */
+  5.5477290039e+03, /* 0x45ad5dd5 */
+  1.9031191406e+03, /* 0x44ede3d0 */
+ -1.3520118713e+02, /* 0xc3073381 */
+];
+
+const QR2: [f32; 6] = [ /* for x in [2.8570,2]=1/[0.3499,0.5] */
+ -1.7838172539e-07, /* 0xb43f8932 */
+ -1.0251704603e-01, /* 0xbdd1f475 */
+ -2.7522056103e+00, /* 0xc0302423 */
+ -1.9663616180e+01, /* 0xc19d4f16 */
+ -4.2325313568e+01, /* 0xc2294d1f */
+ -2.1371921539e+01, /* 0xc1aaf9b2 */
+];
+const QS2: [f32; 6] = [
+  2.9533363342e+01, /* 0x41ec4454 */
+  2.5298155212e+02, /* 0x437cfb47 */
+  7.5750280762e+02, /* 0x443d602e */
+  7.3939318848e+02, /* 0x4438d92a */
+  1.5594900513e+02, /* 0x431bf2f2 */
+ -4.9594988823e+00, /* 0xc09eb437 */
+];
+
+fn qonef(x: f32) -> f32
+{
+    let p: &[f32; 6];
+    let q: &[f32; 6];
+    let s: f32;
+    let r: f32;
+    let z: f32;
+    let mut ix: u32;
+
+    ix = x.to_bits();
+    ix &= 0x7fffffff;
+    if      ix >= 0x41000000 {p = &QR8; q = &QS8;}
+    else if ix >= 0x409173eb {p = &QR5; q = &QS5;}
+    else if ix >= 0x4036d917 {p = &QR3; q = &QS3;}
+    else /*ix >= 0x40000000*/{p = &QR2; q = &QS2;}
+    z = 1.0/(x*x);
+    r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
+    s = 1.0+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));
+    return (0.375 + r/s)/x;
+} 
diff --git a/src/math/jn.rs b/src/math/jn.rs
new file mode 100644
index 0000000..7f7c06f
--- /dev/null
+++ b/src/math/jn.rs
@@ -0,0 +1,338 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/e_jn.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.
+ * ====================================================
+ */
+/*
+ * jn(n, x), yn(n, x)
+ * floating point Bessel's function of the 1st and 2nd kind
+ * of order n
+ *
+ * Special cases:
+ *      y0(0)=y1(0)=yn(n,0) = -inf with division by zero signal;
+ *      y0(-ve)=y1(-ve)=yn(n,-ve) are NaN with invalid signal.
+ * Note 2. About jn(n,x), yn(n,x)
+ *      For n=0, j0(x) is called,
+ *      for n=1, j1(x) is called,
+ *      for n<=x, forward recursion is used starting
+ *      from values of j0(x) and j1(x).
+ *      for n>x, a continued fraction approximation to
+ *      j(n,x)/j(n-1,x) is evaluated and then backward
+ *      recursion is used starting from a supposed value
+ *      for j(n,x). The resulting value of j(0,x) is
+ *      compared with the actual value to correct the
+ *      supposed value of j(n,x).
+ *
+ *      yn(n,x) is similar in all respects, except
+ *      that forward recursion is used for all
+ *      values of n>1.
+ */
+
+use super::{cos, fabs, get_high_word, get_low_word, log, j0, j1, sin, sqrt, y0, y1};
+
+const INVSQRTPI: f64 = 5.64189583547756279280e-01; /* 0x3FE20DD7, 0x50429B6D */
+
+pub fn jn(n: isize, mut x: f64) -> f64
+{
+    let mut ix: u32;
+    let lx: u32;
+    let nm1: isize;
+    let mut i: isize;
+    let mut sign: bool;
+    let mut a: f64;
+    let mut b: f64;
+    let mut temp: f64;
+
+    ix = get_high_word(x);
+    lx = get_low_word(x);
+    sign = (ix>>31) != 0;
+    ix &= 0x7fffffff;
+
+    // -lx == !lx + 1
+    if (ix | (lx|(!lx+1))>>31) > 0x7ff00000 { /* nan */
+        return x;
+    }
+
+    /* J(-n,x) = (-1)^n * J(n, x), J(n, -x) = (-1)^n * J(n, x)
+     * Thus, J(-n,x) = J(n,-x)
+     */
+    /* nm1 = |n|-1 is used instead of |n| to handle n==INT_MIN */
+    if n == 0 {
+        return j0(x);
+    }
+    if n < 0 {
+        nm1 = -(n+1);
+        x = -x;
+        sign = !sign;
+    } else {
+        nm1 = n-1;
+    }
+    if nm1 == 0 {
+        return j1(x);
+    }
+
+    sign &= (n & 1) != 0;  /* even n: 0, odd n: signbit(x) */
+    x = fabs(x);
+    if (ix|lx) == 0 || ix == 0x7ff00000 {  /* if x is 0 or inf */
+        b = 0.0;
+    } else if (nm1 as f64) < x {
+        /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
+        if ix >= 0x52d00000 { /* x > 2**302 */
+            /* (x >> n**2)
+             *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
+             *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
+             *      Let s=sin(x), c=cos(x),
+             *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
+             *
+             *             n    sin(xn)*sqt2    cos(xn)*sqt2
+             *          ----------------------------------
+             *             0     s-c             c+s
+             *             1    -s-c            -c+s
+             *             2    -s+c            -c-s
+             *             3     s+c             c-s
+             */
+            temp = match nm1&3 {
+                0 => -cos(x)+sin(x),
+                1 => -cos(x)-sin(x),
+                2 =>  cos(x)-sin(x),
+                3 | _ => cos(x)+sin(x),
+            };
+            b = INVSQRTPI*temp/sqrt(x);
+        } else {
+            a = j0(x);
+            b = j1(x);
+            i = 0;
+            while i < nm1 {
+                i += 1;
+                temp = b;
+                b = b*(2.0*(i as f64)/x) - a; /* avoid underflow */
+                a = temp;
+            }
+        }
+    } else {
+        if ix < 0x3e100000 { /* x < 2**-29 */
+            /* x is tiny, return the first Taylor expansion of J(n,x)
+             * J(n,x) = 1/n!*(x/2)^n  - ...
+             */
+            if nm1 > 32 {  /* underflow */
+                b = 0.0;
+            } else {
+                temp = x*0.5;
+                b = temp;
+                a = 1.0;
+                i = 2;
+                while i <= nm1 + 1 {
+                    a *= i as f64; /* a = n! */
+                    b *= temp;      /* b = (x/2)^n */
+                    i += 1;
+                }
+                b = b/a;
+            }
+        } else {
+            /* use backward recurrence */
+            /*                      x      x^2      x^2
+             *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
+             *                      2n  - 2(n+1) - 2(n+2)
+             *
+             *                      1      1        1
+             *  (for large x)   =  ----  ------   ------   .....
+             *                      2n   2(n+1)   2(n+2)
+             *                      -- - ------ - ------ -
+             *                       x     x         x
+             *
+             * Let w = 2n/x and h=2/x, then the above quotient
+             * is equal to the continued fraction:
+             *                  1
+             *      = -----------------------
+             *                     1
+             *         w - -----------------
+             *                        1
+             *              w+h - ---------
+             *                     w+2h - ...
+             *
+             * To determine how many terms needed, let
+             * Q(0) = w, Q(1) = w(w+h) - 1,
+             * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
+             * When Q(k) > 1e4      good for single
+             * When Q(k) > 1e9      good for double
+             * When Q(k) > 1e17     good for quadruple
+             */
+            /* determine k */
+            let mut t: f64;
+            let mut q0: f64;
+            let mut q1: f64;
+            let mut w: f64;
+            let h: f64;
+            let mut z: f64;
+            let mut tmp: f64;
+            let nf: f64;
+
+            let mut k: isize;
+
+            nf = (nm1 as f64) + 1.0;
+            w = 2.0*nf/x;
+            h = 2.0/x;
+            z = w+h;
+            q0 = w;
+            q1 = w*z - 1.0;
+            k = 1;
+            while q1 < 1.0e9 {
+                k += 1;
+                z += h;
+                tmp = z*q1 - q0;
+                q0 = q1;
+                q1 = tmp;
+            }
+            t = 0.0;
+            i = k;
+            while i >= 0 {
+                t = 1.0/(2.0*((i as f64)+nf)/x - t);
+                i -= 1;
+            }
+            a = t;
+            b = 1.0;
+            /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
+             *  Hence, if n*(log(2n/x)) > ...
+             *  single 8.8722839355e+01
+             *  double 7.09782712893383973096e+02
+             *  long double 1.1356523406294143949491931077970765006170e+04
+             *  then recurrent value may overflow and the result is
+             *  likely underflow to zero
+             */
+            tmp = nf*log(fabs(w));
+            if tmp < 7.09782712893383973096e+02 {
+                i = nm1;
+                while i > 0 {
+                    temp = b;
+                    b = b*(2.0*(i as f64))/x - a;
+                    a = temp;
+                    i -= 1;
+                }
+            } else {
+                i = nm1;
+                while i > 0 {
+                    temp = b;
+                    b = b*(2.0*(i as f64))/x - a;
+                    a = temp;
+                    /* scale b to avoid spurious overflow */
+                    let x1p500 = f64::from_bits(0x5f30000000000000); // 0x1p500 == 2^500
+                    if b > x1p500 {
+                        a /= b;
+                        t /= b;
+                        b  = 1.0;
+                    }
+                    i -= 1;
+                }
+            }
+            z = j0(x);
+            w = j1(x);
+            if fabs(z) >= fabs(w) {
+                b = t*z/b;
+            } else {
+                b = t*w/a;
+            }
+        }
+    }
+
+    if sign {
+        -b
+    } else {
+        b
+    }
+}
+
+
+pub fn yn(n: isize, x: f64) -> f64
+{
+    let mut ix: u32;
+    let lx: u32;
+    let mut ib: u32;
+    let nm1: isize;
+    let mut sign: bool;
+    let mut i: isize;
+    let mut a: f64;
+    let mut b: f64;
+    let mut temp: f64;
+
+    ix = get_high_word(x);
+    lx = get_low_word(x);
+    sign = (ix>>31) != 0;
+    ix &= 0x7fffffff;
+
+    // -lx == !lx + 1
+    if (ix | (lx|(!lx+1))>>31) > 0x7ff00000 { /* nan */
+        return x;
+    }
+    if sign && (ix|lx) != 0 { /* x < 0 */
+        return 0.0/0.0;
+    }
+    if ix == 0x7ff00000 {
+        return 0.0;
+    }
+
+    if n == 0 {
+        return y0(x);
+    }
+    if n < 0 {
+        nm1 = -(n+1);
+        sign = (n&1) != 0;
+    } else {
+        nm1 = n-1;
+        sign = false;
+    }
+    if nm1 == 0 {
+        if sign {
+            return -y1(x);
+        } else {
+            return y1(x);
+        }
+    }
+
+    if ix >= 0x52d00000 { /* x > 2**302 */
+        /* (x >> n**2)
+         *      Jn(x) = cos(x-(2n+1)*pi/4)*sqrt(2/x*pi)
+         *      Yn(x) = sin(x-(2n+1)*pi/4)*sqrt(2/x*pi)
+         *      Let s=sin(x), c=cos(x),
+         *          xn=x-(2n+1)*pi/4, sqt2 = sqrt(2),then
+         *
+         *             n    sin(xn)*sqt2    cos(xn)*sqt2
+         *          ----------------------------------
+         *             0     s-c             c+s
+         *             1    -s-c            -c+s
+         *             2    -s+c            -c-s
+         *             3     s+c             c-s
+         */
+        temp = match nm1&3 {
+            0 => -sin(x)-cos(x),
+            1 => -sin(x)+cos(x),
+            2 =>  sin(x)+cos(x),
+            3 | _ => sin(x)-cos(x),
+        };
+        b = INVSQRTPI*temp/sqrt(x);
+    } else {
+        a = y0(x);
+        b = y1(x);
+        /* quit if b is -inf */
+        ib = get_high_word(b);
+        i = 0;
+        while i < nm1 && ib != 0xfff00000 {
+            i += 1;
+            temp = b;
+            b = (2.0*(i as f64)/x)*b - a;
+            ib = get_high_word(b);
+            a = temp;
+        }
+    }
+
+    if sign {
+        -b
+    } else {
+        b
+    }
+} 
diff --git a/src/math/jnf.rs b/src/math/jnf.rs
new file mode 100644
index 0000000..4cd848a
--- /dev/null
+++ b/src/math/jnf.rs
@@ -0,0 +1,255 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/e_jnf.c */
+/*
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+/*
+ * ====================================================
+ * 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, j0f, j1f, logf, y0f, y1f};
+
+pub fn jnf(n: isize, mut x: f32) -> f32
+{
+    let mut ix: u32;
+    let mut nm1: isize;
+    let mut sign: bool;
+    let mut i: isize;
+    let mut a: f32;
+    let mut b: f32;
+    let mut temp: f32;
+
+    ix = x.to_bits();
+    sign = (ix>>31) != 0;
+    ix &= 0x7fffffff;
+    if ix > 0x7f800000 { /* nan */
+        return x;
+    }
+
+    /* J(-n,x) = J(n,-x), use |n|-1 to avoid overflow in -n */
+    if n == 0 {
+        return j0f(x);
+    }
+    if n < 0 {
+        nm1 = -(n+1);
+        x = -x;
+        sign = !sign;
+    } else {
+        nm1 = n-1;
+    }
+    if nm1 == 0 {
+        return j1f(x);
+    }
+
+    sign &= (n&1) != 0;  /* even n: 0, odd n: signbit(x) */
+    x = fabsf(x);
+    if ix == 0 || ix == 0x7f800000 {  /* if x is 0 or inf */
+        b = 0.0;
+    } else if (nm1 as f32) < x {
+        /* Safe to use J(n+1,x)=2n/x *J(n,x)-J(n-1,x) */
+        a = j0f(x);
+        b = j1f(x);
+        i = 0;
+        while i < nm1 {
+            i += 1;
+            temp = b;
+            b = b*(2.0*(i as f32)/x) - a;
+            a = temp;
+        }
+    } else {
+        if ix < 0x35800000 { /* x < 2**-20 */
+            /* x is tiny, return the first Taylor expansion of J(n,x)
+             * J(n,x) = 1/n!*(x/2)^n  - ...
+             */
+            if nm1 > 8 {  /* underflow */
+                nm1 = 8;
+            }
+            temp = 0.5 * x;
+            b = temp;
+            a = 1.0;
+            i = 2;
+            while i <= nm1 + 1 {
+                a *= i as f32;    /* a = n! */
+                b *= temp;        /* b = (x/2)^n */
+                i += 1;
+            }
+            b = b/a;
+        } else {
+            /* use backward recurrence */
+            /*                      x      x^2      x^2
+             *  J(n,x)/J(n-1,x) =  ----   ------   ------   .....
+             *                      2n  - 2(n+1) - 2(n+2)
+             *
+             *                      1      1        1
+             *  (for large x)   =  ----  ------   ------   .....
+             *                      2n   2(n+1)   2(n+2)
+             *                      -- - ------ - ------ -
+             *                       x     x         x
+             *
+             * Let w = 2n/x and h=2/x, then the above quotient
+             * is equal to the continued fraction:
+             *                  1
+             *      = -----------------------
+             *                     1
+             *         w - -----------------
+             *                        1
+             *              w+h - ---------
+             *                     w+2h - ...
+             *
+             * To determine how many terms needed, let
+             * Q(0) = w, Q(1) = w(w+h) - 1,
+             * Q(k) = (w+k*h)*Q(k-1) - Q(k-2),
+             * When Q(k) > 1e4      good for single
+             * When Q(k) > 1e9      good for double
+             * When Q(k) > 1e17     good for quadruple
+             */
+            /* determine k */
+            let mut t: f32;
+            let mut q0: f32;
+            let mut q1: f32;
+            let mut w: f32;
+            let h: f32;
+            let mut z: f32;
+            let mut tmp: f32;
+            let nf: f32;
+            let mut k: isize;
+
+            nf = (nm1 as f32)+1.0;
+            w = 2.0*(nf as f32)/x;
+            h = 2.0/x;
+            z = w+h;
+            q0 = w;
+            q1 = w*z - 1.0;
+            k = 1;
+            while q1 < 1.0e4 {
+                k += 1;
+                z += h;
+                tmp = z*q1 - q0;
+                q0 = q1;
+                q1 = tmp;
+            }
+            t = 0.0;
+            i = k;
+            while i >= 0 {
+                t = 1.0/(2.0*((i as f32)+nf)/x-t);
+                i -= 1;
+            }
+            a = t;
+            b = 1.0;
+            /*  estimate log((2/x)^n*n!) = n*log(2/x)+n*ln(n)
+             *  Hence, if n*(log(2n/x)) > ...
+             *  single 8.8722839355e+01
+             *  double 7.09782712893383973096e+02
+             *  long double 1.1356523406294143949491931077970765006170e+04
+             *  then recurrent value may overflow and the result is
+             *  likely underflow to zero
+             */
+            tmp = nf*logf(fabsf(w));
+            if tmp < 88.721679688 {
+                i = nm1;
+                while i > 0 {
+                    temp = b;
+                    b = 2.0*(i as f32)*b/x - a;
+                    a = temp;
+                    i -= 1;
+                }
+            } else {
+                i = nm1;
+                while i > 0 {
+                    temp = b;
+                    b = 2.0*(i as f32)*b/x - a;
+                    a = temp;
+                    /* scale b to avoid spurious overflow */
+                    let x1p60 = f32::from_bits(0x5d800000); // 0x1p60 == 2^60
+                    if b > x1p60 {
+                        a /= b;
+                        t /= b;
+                        b = 1.0;
+                    }
+                    i -= 1;
+                }
+            }
+            z = j0f(x);
+            w = j1f(x);
+            if fabsf(z) >= fabsf(w) {
+                b = t*z/b;
+            } else {
+                b = t*w/a;
+            }
+        }
+    }
+
+    if sign {
+        -b
+    } else {
+        b
+    }
+}
+
+pub fn ynf(n: isize, x: f32) -> f32
+{
+    let mut ix: u32;
+    let mut ib: u32;
+    let nm1: isize;
+    let mut sign: bool;
+    let mut i: isize;
+    let mut a: f32;
+    let mut b: f32;
+    let mut temp: f32;
+
+    ix = x.to_bits();
+    sign = (ix>>31) != 0;
+    ix &= 0x7fffffff;
+    if ix > 0x7f800000 { /* nan */
+        return x;
+    }
+    if sign && ix != 0 { /* x < 0 */
+        return 0.0/0.0;
+    }
+    if ix == 0x7f800000 {
+        return 0.0;
+    }
+
+    if n == 0 {
+        return y0f(x);
+    }
+    if n < 0 {
+        nm1 = -(n+1);
+        sign = (n&1) != 0;
+    } else {
+        nm1 = n-1;
+        sign = false;
+    }
+    if nm1 == 0 {
+        if sign {
+            return -y1f(x);
+        } else {
+            return y1f(x);
+        }
+    }
+
+    a = y0f(x);
+    b = y1f(x);
+    /* quit if b is -inf */
+    ib = b.to_bits();
+    i = 0;
+    while i < nm1 && ib != 0xff800000 {
+        i += 1;
+        temp = b;
+        b = (2.0*(i as f32)/x)*b - a;
+        ib = b.to_bits();
+        a = temp;
+    }
+
+    if sign {
+        -b
+    } else {
+        b
+    }
+}
diff --git a/src/math/lgamma.rs b/src/math/lgamma.rs
new file mode 100644
index 0000000..35b2526
--- /dev/null
+++ b/src/math/lgamma.rs
@@ -0,0 +1,309 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/e_lgamma_r.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.
+ * ====================================================
+ *
+ */
+/* lgamma_r(x, signgamp)
+ * Reentrant version of the logarithm of the Gamma function
+ * with user provide pointer for the sign of Gamma(x).
+ *
+ * Method:
+ *   1. Argument Reduction for 0 < x <= 8
+ *      Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
+ *      reduce x to a number in [1.5,2.5] by
+ *              lgamma(1+s) = log(s) + lgamma(s)
+ *      for example,
+ *              lgamma(7.3) = log(6.3) + lgamma(6.3)
+ *                          = log(6.3*5.3) + lgamma(5.3)
+ *                          = log(6.3*5.3*4.3*3.3*2.3) + lgamma(2.3)
+ *   2. Polynomial approximation of lgamma around its
+ *      minimun ymin=1.461632144968362245 to maintain monotonicity.
+ *      On [ymin-0.23, ymin+0.27] (i.e., [1.23164,1.73163]), use
+ *              Let z = x-ymin;
+ *              lgamma(x) = -1.214862905358496078218 + z^2*poly(z)
+ *      where
+ *              poly(z) is a 14 degree polynomial.
+ *   2. Rational approximation in the primary interval [2,3]
+ *      We use the following approximation:
+ *              s = x-2.0;
+ *              lgamma(x) = 0.5*s + s*P(s)/Q(s)
+ *      with accuracy
+ *              |P/Q - (lgamma(x)-0.5s)| < 2**-61.71
+ *      Our algorithms are based on the following observation
+ *
+ *                             zeta(2)-1    2    zeta(3)-1    3
+ * lgamma(2+s) = s*(1-Euler) + --------- * s  -  --------- * s  + ...
+ *                                 2                 3
+ *
+ *      where Euler = 0.5771... is the Euler constant, which is very
+ *      close to 0.5.
+ *
+ *   3. For x>=8, we have
+ *      lgamma(x)~(x-0.5)log(x)-x+0.5*log(2pi)+1/(12x)-1/(360x**3)+....
+ *      (better formula:
+ *         lgamma(x)~(x-0.5)*(log(x)-1)-.5*(log(2pi)-1) + ...)
+ *      Let z = 1/x, then we approximation
+ *              f(z) = lgamma(x) - (x-0.5)(log(x)-1)
+ *      by
+ *                                  3       5             11
+ *              w = w0 + w1*z + w2*z  + w3*z  + ... + w6*z
+ *      where
+ *              |w - f(z)| < 2**-58.74
+ *
+ *   4. For negative x, since (G is gamma function)
+ *              -x*G(-x)*G(x) = PI/sin(PI*x),
+ *      we have
+ *              G(x) = PI/(sin(PI*x)*(-x)*G(-x))
+ *      since G(-x) is positive, sign(G(x)) = sign(sin(PI*x)) for x<0
+ *      Hence, for x<0, signgam = sign(sin(PI*x)) and
+ *              lgamma(x) = log(|Gamma(x)|)
+ *                        = log(PI/(|x*sin(PI*x)|)) - lgamma(-x);
+ *      Note: one should avoid compute PI*(-x) directly in the
+ *            computation of sin(PI*(-x)).
+ *
+ *   5. Special Cases
+ *              lgamma(2+s) ~ s*(1-Euler) for tiny s
+ *              lgamma(1) = lgamma(2) = 0
+ *              lgamma(x) ~ -log(|x|) for tiny x
+ *              lgamma(0) = lgamma(neg.integer) = inf and raise divide-by-zero
+ *              lgamma(inf) = inf
+ *              lgamma(-inf) = inf (bug for bug compatible with C99!?)
+ *
+ */
+
+use super::{floor, k_cos, k_sin, log};
+
+const PI: f64  =  3.14159265358979311600e+00; /* 0x400921FB, 0x54442D18 */
+const A0: f64  =  7.72156649015328655494e-02; /* 0x3FB3C467, 0xE37DB0C8 */
+const A1: f64  =  3.22467033424113591611e-01; /* 0x3FD4A34C, 0xC4A60FAD */
+const A2: f64  =  6.73523010531292681824e-02; /* 0x3FB13E00, 0x1A5562A7 */
+const A3: f64  =  2.05808084325167332806e-02; /* 0x3F951322, 0xAC92547B */
+const A4: f64  =  7.38555086081402883957e-03; /* 0x3F7E404F, 0xB68FEFE8 */
+const A5: f64  =  2.89051383673415629091e-03; /* 0x3F67ADD8, 0xCCB7926B */
+const A6: f64  =  1.19270763183362067845e-03; /* 0x3F538A94, 0x116F3F5D */
+const A7: f64  =  5.10069792153511336608e-04; /* 0x3F40B6C6, 0x89B99C00 */
+const A8: f64  =  2.20862790713908385557e-04; /* 0x3F2CF2EC, 0xED10E54D */
+const A9: f64  =  1.08011567247583939954e-04; /* 0x3F1C5088, 0x987DFB07 */
+const A10: f64 =  2.52144565451257326939e-05; /* 0x3EFA7074, 0x428CFA52 */
+const A11: f64 =  4.48640949618915160150e-05; /* 0x3F07858E, 0x90A45837 */
+const TC: f64  =  1.46163214496836224576e+00; /* 0x3FF762D8, 0x6356BE3F */
+const TF: f64  = -1.21486290535849611461e-01; /* 0xBFBF19B9, 0xBCC38A42 */
+/* tt = -(tail of TF) */
+const TT: f64  = -3.63867699703950536541e-18; /* 0xBC50C7CA, 0xA48A971F */
+const T0: f64  =  4.83836122723810047042e-01; /* 0x3FDEF72B, 0xC8EE38A2 */
+const T1: f64  = -1.47587722994593911752e-01; /* 0xBFC2E427, 0x8DC6C509 */
+const T2: f64  =  6.46249402391333854778e-02; /* 0x3FB08B42, 0x94D5419B */
+const T3: f64  = -3.27885410759859649565e-02; /* 0xBFA0C9A8, 0xDF35B713 */
+const T4: f64  =  1.79706750811820387126e-02; /* 0x3F9266E7, 0x970AF9EC */
+const T5: f64  = -1.03142241298341437450e-02; /* 0xBF851F9F, 0xBA91EC6A */
+const T6: f64  =  6.10053870246291332635e-03; /* 0x3F78FCE0, 0xE370E344 */
+const T7: f64  = -3.68452016781138256760e-03; /* 0xBF6E2EFF, 0xB3E914D7 */
+const T8: f64  =  2.25964780900612472250e-03; /* 0x3F6282D3, 0x2E15C915 */
+const T9: f64  = -1.40346469989232843813e-03; /* 0xBF56FE8E, 0xBF2D1AF1 */
+const T10: f64 =  8.81081882437654011382e-04; /* 0x3F4CDF0C, 0xEF61A8E9 */
+const T11: f64 = -5.38595305356740546715e-04; /* 0xBF41A610, 0x9C73E0EC */
+const T12: f64 =  3.15632070903625950361e-04; /* 0x3F34AF6D, 0x6C0EBBF7 */
+const T13: f64 = -3.12754168375120860518e-04; /* 0xBF347F24, 0xECC38C38 */
+const T14: f64 =  3.35529192635519073543e-04; /* 0x3F35FD3E, 0xE8C2D3F4 */
+const U0: f64  = -7.72156649015328655494e-02; /* 0xBFB3C467, 0xE37DB0C8 */
+const U1: f64  =  6.32827064025093366517e-01; /* 0x3FE4401E, 0x8B005DFF */
+const U2: f64  =  1.45492250137234768737e+00; /* 0x3FF7475C, 0xD119BD6F */
+const U3: f64  =  9.77717527963372745603e-01; /* 0x3FEF4976, 0x44EA8450 */
+const U4: f64  =  2.28963728064692451092e-01; /* 0x3FCD4EAE, 0xF6010924 */
+const U5: f64  =  1.33810918536787660377e-02; /* 0x3F8B678B, 0xBF2BAB09 */
+const V1: f64  =  2.45597793713041134822e+00; /* 0x4003A5D7, 0xC2BD619C */
+const V2: f64  =  2.12848976379893395361e+00; /* 0x40010725, 0xA42B18F5 */
+const V3: f64  =  7.69285150456672783825e-01; /* 0x3FE89DFB, 0xE45050AF */
+const V4: f64  =  1.04222645593369134254e-01; /* 0x3FBAAE55, 0xD6537C88 */
+const V5: f64  =  3.21709242282423911810e-03; /* 0x3F6A5ABB, 0x57D0CF61 */
+const S0: f64  = -7.72156649015328655494e-02; /* 0xBFB3C467, 0xE37DB0C8 */
+const S1: f64  =  2.14982415960608852501e-01; /* 0x3FCB848B, 0x36E20878 */
+const S2: f64  =  3.25778796408930981787e-01; /* 0x3FD4D98F, 0x4F139F59 */
+const S3: f64  =  1.46350472652464452805e-01; /* 0x3FC2BB9C, 0xBEE5F2F7 */
+const S4: f64  =  2.66422703033638609560e-02; /* 0x3F9B481C, 0x7E939961 */
+const S5: f64  =  1.84028451407337715652e-03; /* 0x3F5E26B6, 0x7368F239 */
+const S6: f64  =  3.19475326584100867617e-05; /* 0x3F00BFEC, 0xDD17E945 */
+const R1: f64  =  1.39200533467621045958e+00; /* 0x3FF645A7, 0x62C4AB74 */
+const R2: f64  =  7.21935547567138069525e-01; /* 0x3FE71A18, 0x93D3DCDC */
+const R3: f64  =  1.71933865632803078993e-01; /* 0x3FC601ED, 0xCCFBDF27 */
+const R4: f64  =  1.86459191715652901344e-02; /* 0x3F9317EA, 0x742ED475 */
+const R5: f64  =  7.77942496381893596434e-04; /* 0x3F497DDA, 0xCA41A95B */
+const R6: f64  =  7.32668430744625636189e-06; /* 0x3EDEBAF7, 0xA5B38140 */
+const W0: f64  =  4.18938533204672725052e-01; /* 0x3FDACFE3, 0x90C97D69 */
+const W1: f64  =  8.33333333333329678849e-02; /* 0x3FB55555, 0x5555553B */
+const W2: f64  = -2.77777777728775536470e-03; /* 0xBF66C16C, 0x16B02E5C */
+const W3: f64  =  7.93650558643019558500e-04; /* 0x3F4A019F, 0x98CF38B6 */
+const W4: f64  = -5.95187557450339963135e-04; /* 0xBF4380CB, 0x8C0FE741 */
+const W5: f64  =  8.36339918996282139126e-04; /* 0x3F4B67BA, 0x4CDAD5D1 */
+const W6: f64  = -1.63092934096575273989e-03; /* 0xBF5AB89D, 0x0B9E43E4 */
+
+/* sin(PI*x) assuming x > 2^-100, if sin(PI*x)==0 the sign is arbitrary */
+fn sin_pi(mut x: f64) -> f64
+{
+    let mut n: isize;
+
+    /* spurious inexact if odd int */
+    x = 2.0*(x*0.5 - floor(x*0.5));  /* x mod 2.0 */
+
+    n = (x*4.0) as isize;
+    n = (n+1)/2;
+    x -= (n as f64)*0.5;
+    x *= PI;
+
+    match n {
+        1   => k_cos(x, 0.0),
+        2   => k_sin(-x, 0.0, 0),
+        3   => -k_cos(x, 0.0),
+        0|_ => k_sin(x, 0.0, 0),
+    }
+}
+
+pub fn lgamma(x: f64) -> f64 {
+    lgamma_r(x).0
+}
+
+pub fn lgamma_r(mut x: f64) -> (f64, isize)
+{
+    let u: u64 = x.to_bits();
+    let mut t: f64;
+    let y: f64;
+    let mut z: f64;
+    let nadj: f64;
+    let p: f64;
+    let p1: f64;
+    let p2: f64;
+    let p3: f64;
+    let q: f64;
+    let mut r: f64;
+    let w: f64;
+    let ix: u32;
+    let sign: bool;
+    let i: isize;
+    let mut signgam: isize;
+
+    /* purge off +-inf, NaN, +-0, tiny and negative arguments */
+    signgam = 1;
+    sign = (u>>63) != 0;
+    ix = ((u>>32) as u32) & 0x7fffffff;
+    if ix >= 0x7ff00000 {
+        return (x*x, signgam);
+    }
+    if ix < (0x3ff-70)<<20 {  /* |x|<2**-70, return -log(|x|) */
+        if sign {
+            x = -x;
+            signgam = -1;
+        }
+        return (-log(x), signgam);
+    }
+    if sign {
+        x = -x;
+        t = sin_pi(x);
+        if t == 0.0 { /* -integer */
+            return (1.0/(x-x), signgam);
+        }
+        if t > 0.0 {
+            signgam = -1;
+        } else {
+            t = -t;
+        }
+        nadj = log(PI/(t*x));
+    } else {
+        nadj = 0.0;
+    }
+
+    /* purge off 1 and 2 */
+    if (ix == 0x3ff00000 || ix == 0x40000000) && (u & 0xffffffff) == 0 {
+        r = 0.0;
+    }
+    /* for x < 2.0 */
+    else if ix < 0x40000000 {
+        if ix <= 0x3feccccc {   /* lgamma(x) = lgamma(x+1)-log(x) */
+            r = -log(x);
+            if ix >= 0x3FE76944 {
+                y = 1.0 - x;
+                i = 0;
+            } else if ix >= 0x3FCDA661 {
+                y = x - (TC-1.0);
+                i = 1;
+            } else {
+                y = x;
+                i = 2;
+            }
+        } else {
+            r = 0.0;
+            if ix >= 0x3FFBB4C3 {  /* [1.7316,2] */
+                y = 2.0 - x;
+                i = 0;
+            } else if ix >= 0x3FF3B4C4 {  /* [1.23,1.73] */
+                y = x - TC;
+                i = 1;
+            } else {
+                y = x - 1.0;
+                i = 2;
+            }
+        }
+        match i {
+            0 => {
+                z = y*y;
+                p1 = A0+z*(A2+z*(A4+z*(A6+z*(A8+z*A10))));
+                p2 = z*(A1+z*(A3+z*(A5+z*(A7+z*(A9+z*A11)))));
+                p = y*p1+p2;
+                r += p-0.5*y;
+            }
+            1 => {
+                z = y*y;
+                w = z*y;
+                p1 = T0+w*(T3+w*(T6+w*(T9 +w*T12)));    /* parallel comp */
+                p2 = T1+w*(T4+w*(T7+w*(T10+w*T13)));
+                p3 = T2+w*(T5+w*(T8+w*(T11+w*T14)));
+                p = z*p1-(TT-w*(p2+y*p3));
+                r += TF + p;
+            }
+            2 => {
+                p1 = y*(U0+y*(U1+y*(U2+y*(U3+y*(U4+y*U5)))));
+                p2 = 1.0+y*(V1+y*(V2+y*(V3+y*(V4+y*V5))));
+                r += -0.5*y + p1/p2;
+            }
+            #[cfg(feature = "checked")]
+            _ => unreachable!(),
+            #[cfg(not(feature = "checked"))]
+            _ => {}
+        }
+    } else if ix < 0x40200000 {  /* x < 8.0 */
+        i = x as isize;
+        y = x - (i as f64);
+        p = y*(S0+y*(S1+y*(S2+y*(S3+y*(S4+y*(S5+y*S6))))));
+        q = 1.0+y*(R1+y*(R2+y*(R3+y*(R4+y*(R5+y*R6)))));
+        r = 0.5*y+p/q;
+        z = 1.0;    /* lgamma(1+s) = log(s) + lgamma(s) */
+        // TODO: In C, this was implemented using switch jumps with fallthrough.
+        // Does this implementation have performance problems?
+        if i >= 7 { z *= y + 6.0; }
+        if i >= 6 { z *= y + 5.0; }
+        if i >= 5 { z *= y + 4.0; }
+        if i >= 4 { z *= y + 3.0; }
+        if i >= 3 {
+            z *= y + 2.0;
+            r += log(z);
+        }
+    } else if ix < 0x43900000 {  /* 8.0 <= x < 2**58 */
+        t = log(x);
+        z = 1.0/x;
+        y = z*z;
+        w = W0+z*(W1+y*(W2+y*(W3+y*(W4+y*(W5+y*W6)))));
+        r = (x-0.5)*(t-1.0)+w;
+    } else {                     /* 2**58 <= x <= inf */
+        r =  x*(log(x)-1.0);
+    }
+    if sign {
+        r = nadj - r;
+    }
+    return (r, signgam);
+}
diff --git a/src/math/lgammaf.rs b/src/math/lgammaf.rs
new file mode 100644
index 0000000..60effa3
--- /dev/null
+++ b/src/math/lgammaf.rs
@@ -0,0 +1,244 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/e_lgammaf_r.c */
+/*
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ */
+/*
+ * ====================================================
+ * 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::{floorf, k_cosf, k_sinf, logf};
+
+const PI: f32  =  3.1415927410e+00; /* 0x40490fdb */
+const A0: f32  =  7.7215664089e-02; /* 0x3d9e233f */
+const A1: f32  =  3.2246702909e-01; /* 0x3ea51a66 */
+const A2: f32  =  6.7352302372e-02; /* 0x3d89f001 */
+const A3: f32  =  2.0580807701e-02; /* 0x3ca89915 */
+const A4: f32  =  7.3855509982e-03; /* 0x3bf2027e */
+const A5: f32  =  2.8905137442e-03; /* 0x3b3d6ec6 */
+const A6: f32  =  1.1927076848e-03; /* 0x3a9c54a1 */
+const A7: f32  =  5.1006977446e-04; /* 0x3a05b634 */
+const A8: f32  =  2.2086278477e-04; /* 0x39679767 */
+const A9: f32  =  1.0801156895e-04; /* 0x38e28445 */
+const A10: f32 =  2.5214456400e-05; /* 0x37d383a2 */
+const A11: f32 =  4.4864096708e-05; /* 0x383c2c75 */
+const TC: f32  =  1.4616321325e+00; /* 0x3fbb16c3 */
+const TF: f32  = -1.2148628384e-01; /* 0xbdf8cdcd */
+/* TT = -(tail of TF) */
+const TT: f32  =  6.6971006518e-09; /* 0x31e61c52 */
+const T0: f32  =  4.8383611441e-01; /* 0x3ef7b95e */
+const T1: f32  = -1.4758771658e-01; /* 0xbe17213c */
+const T2: f32  =  6.4624942839e-02; /* 0x3d845a15 */
+const T3: f32  = -3.2788541168e-02; /* 0xbd064d47 */
+const T4: f32  =  1.7970675603e-02; /* 0x3c93373d */
+const T5: f32  = -1.0314224288e-02; /* 0xbc28fcfe */
+const T6: f32  =  6.1005386524e-03; /* 0x3bc7e707 */
+const T7: f32  = -3.6845202558e-03; /* 0xbb7177fe */
+const T8: f32  =  2.2596477065e-03; /* 0x3b141699 */
+const T9: f32  = -1.4034647029e-03; /* 0xbab7f476 */
+const T10: f32 =  8.8108185446e-04; /* 0x3a66f867 */
+const T11: f32 = -5.3859531181e-04; /* 0xba0d3085 */
+const T12: f32 =  3.1563205994e-04; /* 0x39a57b6b */
+const T13: f32 = -3.1275415677e-04; /* 0xb9a3f927 */
+const T14: f32 =  3.3552918467e-04; /* 0x39afe9f7 */
+const U0: f32  = -7.7215664089e-02; /* 0xbd9e233f */
+const U1: f32  =  6.3282704353e-01; /* 0x3f2200f4 */
+const U2: f32  =  1.4549225569e+00; /* 0x3fba3ae7 */
+const U3: f32  =  9.7771751881e-01; /* 0x3f7a4bb2 */
+const U4: f32  =  2.2896373272e-01; /* 0x3e6a7578 */
+const U5: f32  =  1.3381091878e-02; /* 0x3c5b3c5e */
+const V1: f32  =  2.4559779167e+00; /* 0x401d2ebe */
+const V2: f32  =  2.1284897327e+00; /* 0x4008392d */
+const V3: f32  =  7.6928514242e-01; /* 0x3f44efdf */
+const V4: f32  =  1.0422264785e-01; /* 0x3dd572af */
+const V5: f32  =  3.2170924824e-03; /* 0x3b52d5db */
+const S0: f32  = -7.7215664089e-02; /* 0xbd9e233f */
+const S1: f32  =  2.1498242021e-01; /* 0x3e5c245a */
+const S2: f32  =  3.2577878237e-01; /* 0x3ea6cc7a */
+const S3: f32  =  1.4635047317e-01; /* 0x3e15dce6 */
+const S4: f32  =  2.6642270386e-02; /* 0x3cda40e4 */
+const S5: f32  =  1.8402845599e-03; /* 0x3af135b4 */
+const S6: f32  =  3.1947532989e-05; /* 0x3805ff67 */
+const R1: f32  =  1.3920053244e+00; /* 0x3fb22d3b */
+const R2: f32  =  7.2193557024e-01; /* 0x3f38d0c5 */
+const R3: f32  =  1.7193385959e-01; /* 0x3e300f6e */
+const R4: f32  =  1.8645919859e-02; /* 0x3c98bf54 */
+const R5: f32  =  7.7794247773e-04; /* 0x3a4beed6 */
+const R6: f32  =  7.3266842264e-06; /* 0x36f5d7bd */
+const W0: f32  =  4.1893854737e-01; /* 0x3ed67f1d */
+const W1: f32  =  8.3333335817e-02; /* 0x3daaaaab */
+const W2: f32  = -2.7777778450e-03; /* 0xbb360b61 */
+const W3: f32  =  7.9365057172e-04; /* 0x3a500cfd */
+const W4: f32  = -5.9518753551e-04; /* 0xba1c065c */
+const W5: f32  =  8.3633989561e-04; /* 0x3a5b3dd2 */
+const W6: f32  = -1.6309292987e-03; /* 0xbad5c4e8 */
+
+/* sin(PI*x) assuming x > 2^-100, if sin(PI*x)==0 the sign is arbitrary */
+fn sin_pi(mut x: f32) -> f32
+{
+    let mut y: f64;
+    let mut n: isize;
+
+    /* spurious inexact if odd int */
+    x = 2.0*(x*0.5 - floorf(x*0.5));  /* x mod 2.0 */
+
+    n = (x*4.0) as isize;
+    n = (n+1)/2;
+    y = (x as f64) - (n as f64)*0.5;
+    y *= 3.14159265358979323846;
+    match n {
+        1   => k_cosf(y),
+        2   => k_sinf(-y),
+        3   => -k_cosf(y),
+        0|_ => k_sinf(y),
+    }
+}
+
+pub fn lgammaf(x: f32) -> f32 {
+    lgammaf_r(x).0
+}
+
+pub fn lgammaf_r(mut x: f32) -> (f32, isize)
+{
+    let u = x.to_bits();
+    let mut t: f32;
+    let y: f32;
+    let mut z: f32;
+    let nadj: f32;
+    let p: f32;
+    let p1: f32;
+    let p2: f32;
+    let p3: f32;
+    let q: f32;
+    let mut r: f32;
+    let w: f32;
+    let ix: u32;
+    let i: isize;
+    let sign: bool;
+    let mut signgam: isize;
+
+    /* purge off +-inf, NaN, +-0, tiny and negative arguments */
+    signgam = 1;
+    sign = (u>>31) != 0;
+    ix = u & 0x7fffffff;
+    if ix >= 0x7f800000 {
+        return (x*x, signgam);
+    }
+    if ix < 0x35000000 {  /* |x| < 2**-21, return -log(|x|) */
+        if sign {
+            signgam = -1;
+            x = -x;
+        }
+        return (-logf(x), signgam);
+    }
+    if sign {
+        x = -x;
+        t = sin_pi(x);
+        if t == 0.0 { /* -integer */
+            return (1.0/(x-x), signgam);
+        }
+        if t > 0.0 {
+            signgam = -1;
+        } else {
+            t = -t;
+        }
+        nadj = logf(PI/(t*x));
+    } else {
+        nadj = 0.0;
+    }
+
+    /* purge off 1 and 2 */
+    if ix == 0x3f800000 || ix == 0x40000000 {
+        r = 0.0;
+    }
+    /* for x < 2.0 */
+    else if ix < 0x40000000 {
+        if ix <= 0x3f666666 {  /* lgamma(x) = lgamma(x+1)-log(x) */
+            r = -logf(x);
+            if ix >= 0x3f3b4a20 {
+                y = 1.0 - x;
+                i = 0;
+            } else if ix >= 0x3e6d3308 {
+                y = x - (TC-1.0);
+                i = 1;
+            } else {
+                y = x;
+                i = 2;
+            }
+        } else {
+            r = 0.0;
+            if ix >= 0x3fdda618 {  /* [1.7316,2] */
+                y = 2.0 - x;
+                i = 0;
+            } else if ix >= 0x3F9da620 {  /* [1.23,1.73] */
+                y = x - TC;
+                i = 1;
+            } else {
+                y = x - 1.0;
+                i = 2;
+            }
+        }
+        match i {
+            0 => {
+                z = y*y;
+                p1 = A0+z*(A2+z*(A4+z*(A6+z*(A8+z*A10))));
+                p2 = z*(A1+z*(A3+z*(A5+z*(A7+z*(A9+z*A11)))));
+                p = y*p1+p2;
+                r += p - 0.5*y;
+            }
+            1 => {
+                z = y*y;
+                w = z*y;
+                p1 = T0+w*(T3+w*(T6+w*(T9 +w*T12)));    /* parallel comp */
+                p2 = T1+w*(T4+w*(T7+w*(T10+w*T13)));
+                p3 = T2+w*(T5+w*(T8+w*(T11+w*T14)));
+                p = z*p1-(TT-w*(p2+y*p3));
+                r += TF + p;
+            }
+            2 => {
+                p1 = y*(U0+y*(U1+y*(U2+y*(U3+y*(U4+y*U5)))));
+                p2 = 1.0+y*(V1+y*(V2+y*(V3+y*(V4+y*V5))));
+                r += -0.5*y + p1/p2;
+            }
+            #[cfg(feature = "checked")]
+            _ => unreachable!(),
+            #[cfg(not(feature = "checked"))]
+            _ => {}
+        }
+    } else if ix < 0x41000000 {  /* x < 8.0 */
+        i = x as isize;
+        y = x - (i as f32);
+        p = y*(S0+y*(S1+y*(S2+y*(S3+y*(S4+y*(S5+y*S6))))));
+        q = 1.0+y*(R1+y*(R2+y*(R3+y*(R4+y*(R5+y*R6)))));
+        r = 0.5*y+p/q;
+        z = 1.0;    /* lgamma(1+s) = log(s) + lgamma(s) */
+        // TODO: In C, this was implemented using switch jumps with fallthrough.
+        // Does this implementation have performance problems?
+        if i >= 7 { z *= y + 6.0; }
+        if i >= 6 { z *= y + 5.0; }
+        if i >= 5 { z *= y + 4.0; }
+        if i >= 4 { z *= y + 3.0; }
+        if i >= 3 {
+            z *= y + 2.0;
+            r += logf(z);
+        }
+    } else if ix < 0x5c800000 {  /* 8.0 <= x < 2**58 */
+        t = logf(x);
+        z = 1.0/x;
+        y = z*z;
+        w = W0+z*(W1+y*(W2+y*(W3+y*(W4+y*(W5+y*W6)))));
+        r = (x-0.5)*(t-1.0)+w;
+    } else {                     /* 2**58 <= x <= inf */
+        r =  x*(logf(x)-1.0);
+    }
+    if sign {
+        r = nadj - r;
+    }
+    return (r, signgam);
+}
diff --git a/src/math/mod.rs b/src/math/mod.rs
index be0918f..6e53b02 100644
--- a/src/math/mod.rs
+++ b/src/math/mod.rs
@@ -72,21 +72,33 @@ macro_rules! llvm_intrinsically_optimized {
 // Public modules
 mod acos;
 mod acosf;
+mod acosh;
+mod acoshf;
 mod asin;
 mod asinf;
+mod asinh;
+mod asinhf;
 mod atan;
 mod atan2;
 mod atan2f;
 mod atanf;
+mod atanh;
+mod atanhf;
 mod cbrt;
 mod cbrtf;
 mod ceil;
 mod ceilf;
+mod copysign;
+mod copysignf;
 mod cos;
 mod cosf;
 mod cosh;
 mod coshf;
+mod erf;
+mod erff;
 mod exp;
+mod exp10;
+mod exp10f;
 mod exp2;
 mod exp2f;
 mod expf;
@@ -102,10 +114,22 @@ mod fma;
 mod fmaf;
 mod fmod;
 mod fmodf;
+mod frexp;
+mod frexpf;
 mod hypot;
 mod hypotf;
 mod ldexp;
 mod ldexpf;
+mod ilogb;
+mod ilogbf;
+mod j0;
+mod j0f;
+mod j1;
+mod j1f;
+mod jn;
+mod jnf;
+mod lgamma;
+mod lgammaf;
 mod log;
 mod log10;
 mod log10f;
@@ -114,13 +138,19 @@ mod log1pf;
 mod log2;
 mod log2f;
 mod logf;
+mod modf;
+mod modff;
 mod pow;
 mod powf;
+mod remquo;
+mod remquof;
 mod round;
 mod roundf;
 mod scalbn;
 mod scalbnf;
 mod sin;
+mod sincos;
+mod sincosf;
 mod sinf;
 mod sinh;
 mod sinhf;
@@ -130,27 +160,43 @@ mod tan;
 mod tanf;
 mod tanh;
 mod tanhf;
+mod tgamma;
+mod tgammaf;
 mod trunc;
 mod truncf;
 
 // Use separated imports instead of {}-grouped imports for easier merging.
 pub use self::acos::acos;
 pub use self::acosf::acosf;
+pub use self::acosh::acosh;
+pub use self::acoshf::acoshf;
 pub use self::asin::asin;
 pub use self::asinf::asinf;
+pub use self::asinh::asinh;
+pub use self::asinhf::asinhf;
 pub use self::atan::atan;
 pub use self::atan2::atan2;
 pub use self::atan2f::atan2f;
 pub use self::atanf::atanf;
+pub use self::atanh::atanh;
+pub use self::atanhf::atanhf;
 pub use self::cbrt::cbrt;
 pub use self::cbrtf::cbrtf;
 pub use self::ceil::ceil;
 pub use self::ceilf::ceilf;
+pub use self::copysign::copysign;
+pub use self::copysignf::copysignf;
 pub use self::cos::cos;
 pub use self::cosf::cosf;
 pub use self::cosh::cosh;
 pub use self::coshf::coshf;
+pub use self::erf::erf;
+pub use self::erf::erfc;
+pub use self::erff::erff;
+pub use self::erff::erfcf;
 pub use self::exp::exp;
+pub use self::exp10::exp10;
+pub use self::exp10f::exp10f;
 pub use self::exp2::exp2;
 pub use self::exp2f::exp2f;
 pub use self::expf::expf;
@@ -166,10 +212,30 @@ pub use self::fma::fma;
 pub use self::fmaf::fmaf;
 pub use self::fmod::fmod;
 pub use self::fmodf::fmodf;
+pub use self::frexp::frexp;
+pub use self::frexpf::frexpf;
 pub use self::hypot::hypot;
 pub use self::hypotf::hypotf;
 pub use self::ldexp::ldexp;
 pub use self::ldexpf::ldexpf;
+pub use self::ilogb::ilogb;
+pub use self::ilogbf::ilogbf;
+pub use self::j0::j0;
+pub use self::j0::y0;
+pub use self::j0f::j0f;
+pub use self::j0f::y0f;
+pub use self::j1::j1;
+pub use self::j1::y1;
+pub use self::j1f::j1f;
+pub use self::j1f::y1f;
+pub use self::jn::jn;
+pub use self::jn::yn;
+pub use self::jnf::jnf;
+pub use self::jnf::ynf;
+pub use self::lgamma::lgamma;
+pub use self::lgamma::lgamma_r;
+pub use self::lgammaf::lgammaf;
+pub use self::lgammaf::lgammaf_r;
 pub use self::log::log;
 pub use self::log10::log10;
 pub use self::log10f::log10f;
@@ -178,13 +244,19 @@ pub use self::log1pf::log1pf;
 pub use self::log2::log2;
 pub use self::log2f::log2f;
 pub use self::logf::logf;
+pub use self::modf::modf;
+pub use self::modff::modff;
 pub use self::pow::pow;
 pub use self::powf::powf;
+pub use self::remquo::remquo;
+pub use self::remquof::remquof;
 pub use self::round::round;
 pub use self::roundf::roundf;
 pub use self::scalbn::scalbn;
 pub use self::scalbnf::scalbnf;
 pub use self::sin::sin;
+pub use self::sincos::sincos;
+pub use self::sincosf::sincosf;
 pub use self::sinf::sinf;
 pub use self::sinh::sinh;
 pub use self::sinhf::sinhf;
@@ -194,6 +266,8 @@ pub use self::tan::tan;
 pub use self::tanf::tanf;
 pub use self::tanh::tanh;
 pub use self::tanhf::tanhf;
+pub use self::tgamma::tgamma;
+pub use self::tgammaf::tgammaf;
 pub use self::trunc::trunc;
 pub use self::truncf::truncf;
 
diff --git a/src/math/modf.rs b/src/math/modf.rs
new file mode 100644
index 0000000..1ff8ee1
--- /dev/null
+++ b/src/math/modf.rs
@@ -0,0 +1,33 @@
+pub fn modf(x: f64) -> (f64, f64) {
+    let rv2: f64;
+    let mut u = x.to_bits();
+    let mask: u64;
+    let e = ((u>>52 & 0x7ff) as isize) - 0x3ff;
+
+    /* no fractional part */
+    if e >= 52 {
+        rv2 = x;
+        if e == 0x400 && (u<<12) != 0 { /* nan */
+            return (x, rv2);
+        }
+        u &= 1<<63;
+        return (f64::from_bits(u), rv2);
+    }
+
+    /* no integral part*/
+    if e < 0 {
+        u &= 1<<63;
+        rv2 = f64::from_bits(u);
+        return (x, rv2);
+    }
+
+    mask = ((!0)>>12)>>e;
+    if (u & mask) == 0 {
+        rv2 = x;
+        u &= 1<<63;
+        return (f64::from_bits(u), rv2);
+    }
+    u &= !mask;
+    rv2 = f64::from_bits(u);
+    return (x - rv2, rv2);
+} 
diff --git a/src/math/modff.rs b/src/math/modff.rs
new file mode 100644
index 0000000..5250e8d
--- /dev/null
+++ b/src/math/modff.rs
@@ -0,0 +1,32 @@
+pub fn modff(x: f32) -> (f32, f32) {
+    let rv2: f32;
+    let mut u: u32 = x.to_bits();
+    let mask: u32;
+    let e = ((u>>23 & 0xff) as isize) - 0x7f;
+
+    /* no fractional part */
+    if e >= 23 {
+        rv2 = x;
+        if e == 0x80 && (u<<9) != 0 { /* nan */
+            return (x, rv2);
+        }
+        u &= 0x80000000;
+        return (f32::from_bits(u), rv2);
+    }
+    /* no integral part */
+    if e < 0 {
+        u &= 0x80000000;
+        rv2 = f32::from_bits(u);
+        return (x, rv2);
+    }
+
+    mask = 0x007fffff>>e;
+    if (u & mask) == 0 {
+        rv2 = x;
+        u &= 0x80000000;
+        return (f32::from_bits(u), rv2);
+    }
+    u &= !mask;
+    rv2 = f32::from_bits(u);
+    return (x - rv2, rv2);
+}
diff --git a/src/math/remquo.rs b/src/math/remquo.rs
new file mode 100644
index 0000000..98f4b38
--- /dev/null
+++ b/src/math/remquo.rs
@@ -0,0 +1,98 @@
+pub fn remquo(mut x: f64, mut y: f64) -> (f64, isize)
+{
+    let ux: u64 = x.to_bits();
+    let mut uy: u64 = y.to_bits();
+    let mut ex = ((ux>>52) & 0x7ff) as isize;
+    let mut ey = ((uy>>52) & 0x7ff) as isize;
+    let sx = (ux>>63) != 0;
+    let sy = (uy>>63) != 0;
+    let mut q: u32;
+    let mut i: u64;
+    let mut uxi: u64 = ux;
+
+    if (uy<<1) == 0 || y.is_nan() || ex == 0x7ff {
+        return ((x*y)/(x*y), 0);
+    }
+    if (ux<<1) == 0 {
+        return (x, 0);
+    }
+
+    /* normalize x and y */
+    if ex == 0 {
+        i = uxi << 12;
+        while (i>>63) == 0 {
+            ex -= 1;
+            i <<= 1;
+        }
+        uxi <<= -ex + 1;
+    } else {
+        uxi &= (!0) >> 12;
+        uxi |= 1 << 52;
+    }
+    if ey == 0 {
+        i = uy<<12;
+        while (i>>63) == 0 {
+            ey -= 1;
+            i <<= 1;
+        }
+        uy <<= -ey + 1;
+    } else {
+        uy &= (!0) >> 12;
+        uy |= 1 << 52;
+    }
+
+    q = 0;
+
+    if ex+1 != ey {
+        if ex < ey {
+            return (x, 0);
+        }
+        /* x mod y */
+        while ex > ey {
+            i = uxi - uy;
+            if (i>>63) == 0 {
+                uxi = i;
+                q += 1;
+            }
+            uxi <<= 1;
+            q <<= 1;
+            ex -= 1;
+        }
+        i = uxi - uy;
+        if (i>>63) == 0 {
+            uxi = i;
+            q += 1;
+        }
+        if uxi == 0 {
+            ex = -60;
+        } else {
+            while (uxi>>52) == 0 {
+                uxi <<= 1;
+                ex -= 1;
+            }
+        }
+    }
+
+    /* scale result and decide between |x| and |x|-|y| */
+    if ex > 0 {
+        uxi -= 1 << 52;
+        uxi |= (ex as u64) << 52;
+    } else {
+        uxi >>= -ex + 1;
+    }
+    x = f64::from_bits(uxi);
+    if sy {
+        y = -y;
+    }
+    if ex == ey || (ex+1 == ey && (2.0*x > y || (2.0*x == y && (q%2) != 0))) {
+        x -= y;
+        q += 1;
+    }
+    q &= 0x7fffffff;
+    let quo = if sx ^ sy { -(q as isize) } else { q as isize };
+    if sx {
+        (-x, quo)
+    } else {
+        (x, quo)
+    }
+}
diff --git a/src/math/remquof.rs b/src/math/remquof.rs
new file mode 100644
index 0000000..4307e19
--- /dev/null
+++ b/src/math/remquof.rs
@@ -0,0 +1,97 @@
+pub fn remquof(mut x: f32, mut y: f32) -> (f32, isize)
+{
+    let ux: u32 = x.to_bits();
+    let mut uy: u32 = y.to_bits();
+    let mut ex = ((ux>>23) & 0xff) as isize;
+    let mut ey = ((uy>>23) & 0xff) as isize;
+    let sx = (ux>>31) != 0;
+    let sy = (uy>>31) != 0;
+    let mut q: u32;
+    let mut i: u32;
+    let mut uxi: u32 = ux;
+
+    if (uy<<1) == 0 || y.is_nan() || ex == 0xff {
+        return ((x*y)/(x*y), 0);
+    }
+    if (ux<<1) == 0 {
+        return (x, 0);
+    }
+
+    /* normalize x and y */
+    if ex == 0 {
+        i = uxi<<9;
+        while (i>>31) == 0 {
+            ex -= 1;
+            i <<= 1;
+        }
+        uxi <<= -ex + 1;
+    } else {
+        uxi &= (!0) >> 9;
+        uxi |= 1 << 23;
+    }
+    if ey == 0 {
+        i = uy<<9;
+        while (i>>31) == 0 {
+            ey -= 1;
+            i <<= 1;
+        }
+        uy <<= -ey + 1;
+    } else {
+        uy &= (!0) >> 9;
+        uy |= 1 << 23;
+    }
+
+    q = 0;
+    if ex+1 != ey {
+        if ex < ey {
+            return (x, 0);
+        }
+        /* x mod y */
+        while ex > ey {
+            i = uxi - uy;
+            if (i>>31) == 0 {
+                uxi = i;
+                q += 1;
+            }
+            uxi <<= 1;
+            q <<= 1;
+            ex -= 1;
+        }
+        i = uxi - uy;
+        if (i>>31) == 0 {
+            uxi = i;
+            q += 1;
+        }
+        if uxi == 0 {
+            ex = -30;
+        } else {
+            while (uxi>>23) == 0 {
+                uxi <<= 1;
+                ex -= 1;
+            }
+        }
+    }
+
+    /* scale result and decide between |x| and |x|-|y| */
+    if ex > 0 {
+        uxi -= 1 << 23;
+        uxi |= (ex as u32) << 23;
+    } else {
+        uxi >>= -ex + 1;
+    }
+    x = f32::from_bits(uxi);
+    if sy {
+        y = -y;
+    }
+    if ex == ey || (ex+1 == ey && (2.0*x > y || (2.0*x == y && (q%2) != 0))) {
+        x -= y;
+        q += 1;
+    }
+    q &= 0x7fffffff;
+    let quo = if sx^sy { -(q as isize) } else { q as isize };
+    if sx {
+        (-x, quo)
+    } else {
+        (x, quo)
+    }
+}
diff --git a/src/math/sincos.rs b/src/math/sincos.rs
new file mode 100644
index 0000000..c15ee46
--- /dev/null
+++ b/src/math/sincos.rs
@@ -0,0 +1,60 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_sin.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 super::{get_high_word, k_cos, k_sin, rem_pio2};
+
+pub fn sincos(x: f64) -> (f64, f64)
+{
+    let s: f64;
+    let c: f64;
+    let mut ix: u32;
+
+    ix = get_high_word(x);
+    ix &= 0x7fffffff;
+
+    /* |x| ~< pi/4 */
+    if ix <= 0x3fe921fb {
+        /* if |x| < 2**-27 * sqrt(2) */
+        if ix < 0x3e46a09e {
+            /* raise inexact if x!=0 and underflow if subnormal */
+            let x1p120 = f64::from_bits(0x4770000000000000); // 0x1p120 == 2^120
+            if ix < 0x00100000 {
+                force_eval!(x/x1p120);
+            } else {
+                force_eval!(x+x1p120);
+            }
+            return (x, 1.0);
+        }
+        return (k_sin(x, 0.0, 0), k_cos(x, 0.0));
+    }
+
+    /* sincos(Inf or NaN) is NaN */
+    if ix >= 0x7ff00000 {
+        let rv = x - x;
+        return (rv, rv);
+    }
+
+    /* argument reduction needed */
+    let (n, y0, y1) = rem_pio2(x);
+    s = k_sin(y0, y1, 1);
+    c = k_cos(y0, y1);
+    match n&3 {
+        0 => (s, c),
+        1 => (c, -s),
+        2 => (-s, -c),
+        3 => (-c, s),
+        #[cfg(feature = "checked")]
+        _ => unreachable!(),
+        #[cfg(not(feature = "checked"))]
+        _ => (0.0, 1.0),
+    }
+}
diff --git a/src/math/sincosf.rs b/src/math/sincosf.rs
new file mode 100644
index 0000000..911421d
--- /dev/null
+++ b/src/math/sincosf.rs
@@ -0,0 +1,122 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/s_sinf.c */
+/*
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
+ * 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.
+ * ====================================================
+ */
+
+use super::{k_cosf, k_sinf, rem_pio2f};
+
+/* Small multiples of pi/2 rounded to double precision. */
+const PI_2: f32 = 0.5 * 3.1415926535897931160E+00;
+const S1PIO2: f32 = 1.0*PI_2; /* 0x3FF921FB, 0x54442D18 */
+const S2PIO2: f32 = 2.0*PI_2; /* 0x400921FB, 0x54442D18 */
+const S3PIO2: f32 = 3.0*PI_2; /* 0x4012D97C, 0x7F3321D2 */
+const S4PIO2: f32 = 4.0*PI_2; /* 0x401921FB, 0x54442D18 */
+
+pub fn sincosf(x: f32) -> (f32, f32)
+{
+    let s: f32;
+    let c: f32;
+    let mut ix: u32;
+    let sign: bool;
+
+    ix = x.to_bits();
+    sign = (ix >> 31) != 0;
+    ix &= 0x7fffffff;
+
+    /* |x| ~<= pi/4 */
+    if ix <= 0x3f490fda {
+        /* |x| < 2**-12 */
+        if ix < 0x39800000 {
+            /* raise inexact if x!=0 and underflow if subnormal */
+
+            let x1p120 = f32::from_bits(0x7b800000); // 0x1p120 == 2^120
+            if ix < 0x00100000 {
+                force_eval!(x/x1p120);
+            } else {
+                force_eval!(x+x1p120);
+            }
+            return (x, 1.0);
+        }
+        return (k_sinf(x as f64), k_cosf(x as f64));
+    }
+
+    /* |x| ~<= 5*pi/4 */
+    if ix <= 0x407b53d1 {
+        if ix <= 0x4016cbe3 {  /* |x| ~<= 3pi/4 */
+            if sign {
+                s = -k_cosf((x + S1PIO2) as f64);
+                c = k_sinf((x + S1PIO2) as f64);
+            } else {
+                s = k_cosf((S1PIO2 - x) as f64);
+                c = k_sinf((S1PIO2 - x) as f64);
+            }
+        }
+        /* -sin(x+c) is not correct if x+c could be 0: -0 vs +0 */
+        else {
+            if sign {
+                s = k_sinf((x + S2PIO2) as f64);
+                c = k_cosf((x + S2PIO2) as f64);
+            } else {
+                s = k_sinf((x - S2PIO2) as f64);
+                c = k_cosf((x - S2PIO2) as f64);
+            }
+        }
+
+        return (s, c);
+    }
+
+    /* |x| ~<= 9*pi/4 */
+    if ix <= 0x40e231d5 {
+        if ix <= 0x40afeddf {  /* |x| ~<= 7*pi/4 */
+            if sign {
+                s = k_cosf((x + S3PIO2) as f64);
+                c = -k_sinf((x + S3PIO2) as f64);
+            } else {
+                s = -k_cosf((x - S3PIO2) as f64);
+                c = k_sinf((x - S3PIO2) as f64);
+            }
+        } else {
+            if sign {
+                s = k_cosf((x + S4PIO2) as f64);
+                c = k_sinf((x + S4PIO2) as f64);
+            } else {
+                s = k_cosf((x - S4PIO2) as f64);
+                c = k_sinf((x - S4PIO2) as f64);
+            }
+        }
+
+        return (s, c);
+    }
+
+    /* sin(Inf or NaN) is NaN */
+    if ix >= 0x7f800000 {
+        let rv = x - x;
+        return (rv, rv);
+    }
+
+    /* general argument reduction needed */
+    let (n, y) = rem_pio2f(x);
+    s = k_sinf(y);
+    c = k_cosf(y);
+    match n&3 {
+        0 => (s, c),
+        1 => (c, -s),
+        2 => (-s, -c),
+        3 => (-c, s),
+        #[cfg(feature = "checked")]
+        _ => unreachable!(),
+        #[cfg(not(feature = "checked"))]
+        _ => (0.0, 1.0),
+    }
+}
diff --git a/src/math/tgamma.rs b/src/math/tgamma.rs
new file mode 100644
index 0000000..598f46f
--- /dev/null
+++ b/src/math/tgamma.rs
@@ -0,0 +1,179 @@
+/*
+"A Precision Approximation of the Gamma Function" - Cornelius Lanczos (1964)
+"Lanczos Implementation of the Gamma Function" - Paul Godfrey (2001)
+"An Analysis of the Lanczos Gamma Approximation" - Glendon Ralph Pugh (2004)
+
+approximation method:
+
+                        (x - 0.5)         S(x)
+Gamma(x) = (x + g - 0.5)         *  ----------------
+                                    exp(x + g - 0.5)
+
+with
+                 a1      a2      a3            aN
+S(x) ~= [ a0 + ----- + ----- + ----- + ... + ----- ]
+               x + 1   x + 2   x + 3         x + N
+
+with a0, a1, a2, a3,.. aN constants which depend on g.
+
+for x < 0 the following reflection formula is used:
+
+Gamma(x)*Gamma(-x) = -pi/(x sin(pi x))
+
+most ideas and constants are from boost and python
+*/
+extern crate core;
+use super::{exp, floor, k_cos, k_sin, pow};
+
+const PI: f64 = 3.141592653589793238462643383279502884;
+
+/* sin(pi x) with x > 0x1p-100, if sin(pi*x)==0 the sign is arbitrary */
+fn sinpi(mut x: f64) -> f64
+{
+    let mut n: isize;
+
+    /* argument reduction: x = |x| mod 2 */
+    /* spurious inexact when x is odd int */
+    x = x * 0.5;
+    x = 2.0 * (x - floor(x));
+
+    /* reduce x into [-.25,.25] */
+    n = (4.0 * x) as isize;
+    n = (n+1)/2;
+    x -= (n as f64) * 0.5;
+
+    x *= PI;
+    match n {
+        1   => k_cos(x, 0.0),
+        2   => k_sin(-x, 0.0, 0),
+        3   => -k_cos(x, 0.0),
+        0|_ => k_sin(x, 0.0, 0),
+    }
+}
+
+const N: usize = 12;
+//static const double g = 6.024680040776729583740234375;
+const GMHALF: f64 = 5.524680040776729583740234375;
+const SNUM: [f64; N+1] = [
+    23531376880.410759688572007674451636754734846804940,
+    42919803642.649098768957899047001988850926355848959,
+    35711959237.355668049440185451547166705960488635843,
+    17921034426.037209699919755754458931112671403265390,
+    6039542586.3520280050642916443072979210699388420708,
+    1439720407.3117216736632230727949123939715485786772,
+    248874557.86205415651146038641322942321632125127801,
+    31426415.585400194380614231628318205362874684987640,
+    2876370.6289353724412254090516208496135991145378768,
+    186056.26539522349504029498971604569928220784236328,
+    8071.6720023658162106380029022722506138218516325024,
+    210.82427775157934587250973392071336271166969580291,
+    2.5066282746310002701649081771338373386264310793408,
+];
+const SDEN: [f64; N+1] = [
+    0.0, 39916800.0, 120543840.0, 150917976.0, 105258076.0,
+    45995730.0, 13339535.0, 2637558.0, 357423.0, 32670.0, 1925.0, 66.0, 1.0,
+];
+/* n! for small integer n */
+const FACT: [f64; 23] = [
+    1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0, 3628800.0,
+    39916800.0, 479001600.0, 6227020800.0, 87178291200.0, 1307674368000.0,
+    20922789888000.0, 355687428096000.0, 6402373705728000.0, 121645100408832000.0,
+    2432902008176640000.0, 51090942171709440000.0, 1124000727777607680000.0,
+];
+
+/* S(x) rational function for positive x */
+fn s(x: f64) -> f64
+{
+    let mut num: f64 = 0.0;
+    let mut den: f64 = 0.0;
+
+    /* to avoid overflow handle large x differently */
+    if x < 8.0 {
+        for i in (0..=N).rev() {
+            num = num * x + SNUM[i];
+            den = den * x + SDEN[i];
+        }
+    } else {
+        for i in 0..=N {
+            num = num / x + SNUM[i];
+            den = den / x + SDEN[i];
+        }
+    }
+    return num/den;
+}
+
+pub fn tgamma(mut x: f64) -> f64
+{
+    let u: u64 = x.to_bits();
+    let absx: f64;
+    let mut y: f64;
+    let mut dy: f64;
+    let mut z: f64;
+    let mut r: f64;
+    let ix: u32 = ((u >> 32) as u32) & 0x7fffffff;
+    let sign: bool = (u>>64) != 0;
+
+    /* special cases */
+    if ix >= 0x7ff00000 {
+        /* tgamma(nan)=nan, tgamma(inf)=inf, tgamma(-inf)=nan with invalid */
+        return x + core::f64::INFINITY;
+    }
+    if ix < ((0x3ff-54)<<20) {
+        /* |x| < 2^-54: tgamma(x) ~ 1/x, +-0 raises div-by-zero */
+        return 1.0/x;
+    }
+
+    /* integer arguments */
+    /* raise inexact when non-integer */
+    if x == floor(x) {
+        if sign {
+            return 0.0/0.0;
+        }
+        if x <= FACT.len() as f64 {
+            return FACT[(x as usize) - 1];
+        }
+    }
+
+    /* x >= 172: tgamma(x)=inf with overflow */
+    /* x =< -184: tgamma(x)=+-0 with underflow */
+    if ix >= 0x40670000 { /* |x| >= 184 */
+        if sign {
+            let x1p_126 = f64::from_bits(0x3810000000000000); // 0x1p-126 == 2^-126
+            force_eval!((x1p_126/x) as f32);
+            if floor(x) * 0.5 == floor(x * 0.5) {
+                return 0.0;
+            } else {
+                return -0.0;
+            }
+        }
+        let x1p1023 = f64::from_bits(0x7fe0000000000000); // 0x1p1023 == 2^1023
+        x *= x1p1023;
+        return x;
+    }
+
+    absx = if sign { -x } else { x };
+
+    /* handle the error of x + g - 0.5 */
+    y = absx + GMHALF;
+    if absx > GMHALF {
+        dy = y - absx;
+        dy -= GMHALF;
+    } else {
+        dy = y - GMHALF;
+        dy -= absx;
+    }
+
+    z = absx - 0.5;
+    r = s(absx) * exp(-y);
+    if x < 0.0 {
+        /* reflection formula for negative x */
+        /* sinpi(absx) is not 0, integers are already handled */
+        r = -PI / (sinpi(absx) * absx * r);
+        dy = -dy;
+        z = -z;
+    }
+    r += dy * (GMHALF+0.5) * r / y;
+    z = pow(y, 0.5*z);
+    y = r * z * z;
+    return y;
+}
diff --git a/src/math/tgammaf.rs b/src/math/tgammaf.rs
new file mode 100644
index 0000000..b9c799c
--- /dev/null
+++ b/src/math/tgammaf.rs
@@ -0,0 +1,5 @@
+use super::{tgamma};
+
+pub fn tgammaf(x: f32) -> f32 {
+    tgamma(x as f64) as f32
+}