Ported several remaining math functions from musl

Please note that these aren't tested yet.

Author: m1el

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;
+}

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;
+}

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;
+}

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
+    }
+}

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
+    }
+}

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
+    }
+}

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
+    }
+}

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)
+}

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)
+}