Merge branch 'master' into separate-import

Author: japaric

libm/src/lib.rs

@@ -87,7 +87,6 @@ pub trait F32Ext: private::Sealed {
     #[cfg(todo)]
     fn sin(self) -> Self;
 
-    #[cfg(todo)]
     fn cos(self) -> Self;
 
     #[cfg(todo)]
@@ -252,7 +251,6 @@ impl F32Ext for f32 {
         sinf(self)
     }
 
-    #[cfg(todo)]
     #[inline]
     fn cos(self) -> Self {
         cosf(self)

libm/src/math/cosf.rs

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

libm/src/math/expf.rs

@@ -1,3 +1,18 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/e_expf.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::scalbnf;
 
 const HALF: [f32; 2] = [0.5, -0.5];

libm/src/math/k_cosf.rs

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

libm/src/math/k_sinf.rs

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

libm/src/math/log10.rs

@@ -1,3 +1,22 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/e_log10.c */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+/*
+ * Return the base 10 logarithm of x.  See log.c for most comments.
+ *
+ * Reduce x to 2^k (1+f) and calculate r = log(1+f) - f + f*f/2
+ * as in log.c, then combine and scale in extra precision:
+ *    log10(x) = (f - f*f/2 + r)/log(10) + k*log10(2)
+ */
+
 use core::f64;
 
 const IVLN10HI: f64 = 4.34294481878168880939e-01; /* 0x3fdbcb7b, 0x15200000 */

libm/src/math/log10f.rs

@@ -1,3 +1,18 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/e_log10f.c */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+/*
+ * See comments in log10.c.
+ */
+
 use core::f32;
 
 const IVLN10HI: f32 = 4.3432617188e-01; /* 0x3ede6000 */

libm/src/math/log2.rs

@@ -1,3 +1,22 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/e_log2.c */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+/*
+ * Return the base 2 logarithm of x.  See log.c for most comments.
+ *
+ * Reduce x to 2^k (1+f) and calculate r = log(1+f) - f + f*f/2
+ * as in log.c, then combine and scale in extra precision:
+ *    log2(x) = (f - f*f/2 + r)/log(2) + k
+ */
+
 use core::f64;
 
 const IVLN2HI: f64 = 1.44269504072144627571e+00; /* 0x3ff71547, 0x65200000 */

libm/src/math/log2f.rs

@@ -1,3 +1,18 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/e_log2f.c */
+/*
+ * ====================================================
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+ *
+ * Developed at SunPro, a Sun Microsystems, Inc. business.
+ * Permission to use, copy, modify, and distribute this
+ * software is freely granted, provided that this notice
+ * is preserved.
+ * ====================================================
+ */
+/*
+ * See comments in log2.c.
+ */
+
 use core::f32;
 
 const IVLN2HI: f32 = 1.4428710938e+00; /* 0x3fb8b000 */

libm/src/math/logf.rs

@@ -1,3 +1,18 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/e_logf.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.
+ * ====================================================
+ */
+
 const LN2_HI: f32 = 6.9313812256e-01; /* 0x3f317180 */
 const LN2_LO: f32 = 9.0580006145e-06; /* 0x3717f7d1 */
 /* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */

libm/src/math/mod.rs

@@ -7,6 +7,7 @@ macro_rules! force_eval {
 }
 
 mod ceilf;
+mod cosf;
 mod expf;
 mod fabs;
 mod fabsf;
@@ -33,10 +34,9 @@ mod sqrtf;
 mod trunc;
 mod truncf;
 
-//mod service;
-
 // Use separated imports instead of {}-grouped imports for easier merging.
 pub use self::ceilf::ceilf;
+pub use self::cosf::cosf;
 pub use self::expf::expf;
 pub use self::fabs::fabs;
 pub use self::fabsf::fabsf;
@@ -63,6 +63,13 @@ pub use self::sqrtf::sqrtf;
 pub use self::trunc::trunc;
 pub use self::truncf::truncf;
 
+mod k_cosf;
+mod k_sinf;
+mod rem_pio2_large;
+mod rem_pio2f;
+
+use self::{k_cosf::k_cosf, k_sinf::k_sinf, rem_pio2_large::rem_pio2_large, rem_pio2f::rem_pio2f};
+
 fn isnanf(x: f32) -> bool {
     x.to_bits() & 0x7fffffff > 0x7f800000
 }

libm/src/math/powf.rs

@@ -1,3 +1,18 @@
+/* origin: FreeBSD /usr/src/lib/msun/src/e_powf.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, scalbnf, sqrtf};
 
 const BP: [f32; 2] = [1.0, 1.5];