fix bit shifting error
This commit is contained in:
Lukas Wirth
commit
a4602adbb5
10 changed files with 905 additions and 19 deletions
|
|
@ -70,7 +70,6 @@ pub trait F32Ext: private::Sealed {
|
|||
|
||||
fn exp(self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
fn exp2(self) -> Self;
|
||||
|
||||
fn ln(self) -> Self;
|
||||
|
|
@ -81,7 +80,6 @@ pub trait F32Ext: private::Sealed {
|
|||
|
||||
fn log10(self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
fn cbrt(self) -> Self;
|
||||
|
||||
fn hypot(self, other: Self) -> Self;
|
||||
|
|
@ -215,7 +213,6 @@ impl F32Ext for f32 {
|
|||
expf(self)
|
||||
}
|
||||
|
||||
#[cfg(todo)]
|
||||
#[inline]
|
||||
fn exp2(self) -> Self {
|
||||
exp2f(self)
|
||||
|
|
@ -241,7 +238,6 @@ impl F32Ext for f32 {
|
|||
log10f(self)
|
||||
}
|
||||
|
||||
#[cfg(todo)]
|
||||
#[inline]
|
||||
fn cbrt(self) -> Self {
|
||||
cbrtf(self)
|
||||
|
|
@ -389,7 +385,6 @@ pub trait F64Ext: private::Sealed {
|
|||
|
||||
fn exp(self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
fn exp2(self) -> Self;
|
||||
|
||||
fn ln(self) -> Self;
|
||||
|
|
@ -400,7 +395,6 @@ pub trait F64Ext: private::Sealed {
|
|||
|
||||
fn log10(self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
fn cbrt(self) -> Self;
|
||||
|
||||
fn hypot(self, other: Self) -> Self;
|
||||
|
|
@ -417,7 +411,6 @@ pub trait F64Ext: private::Sealed {
|
|||
#[cfg(todo)]
|
||||
fn asin(self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
fn acos(self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
|
|
@ -534,7 +527,6 @@ impl F64Ext for f64 {
|
|||
exp(self)
|
||||
}
|
||||
|
||||
#[cfg(todo)]
|
||||
#[inline]
|
||||
fn exp2(self) -> Self {
|
||||
exp2(self)
|
||||
|
|
@ -560,7 +552,6 @@ impl F64Ext for f64 {
|
|||
log10(self)
|
||||
}
|
||||
|
||||
#[cfg(todo)]
|
||||
#[inline]
|
||||
fn cbrt(self) -> Self {
|
||||
cbrt(self)
|
||||
|
|
@ -595,7 +586,6 @@ impl F64Ext for f64 {
|
|||
asin(self)
|
||||
}
|
||||
|
||||
#[cfg(todo)]
|
||||
#[inline]
|
||||
fn acos(self) -> Self {
|
||||
acos(self)
|
||||
|
|
|
|||
108
library/compiler-builtins/libm/src/math/acos.rs
Normal file
108
library/compiler-builtins/libm/src/math/acos.rs
Normal file
|
|
@ -0,0 +1,108 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/e_acos.c */
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* Developed at SunSoft, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
/* acos(x)
|
||||
* Method :
|
||||
* acos(x) = pi/2 - asin(x)
|
||||
* acos(-x) = pi/2 + asin(x)
|
||||
* For |x|<=0.5
|
||||
* acos(x) = pi/2 - (x + x*x^2*R(x^2)) (see asin.c)
|
||||
* For x>0.5
|
||||
* acos(x) = pi/2 - (pi/2 - 2asin(sqrt((1-x)/2)))
|
||||
* = 2asin(sqrt((1-x)/2))
|
||||
* = 2s + 2s*z*R(z) ...z=(1-x)/2, s=sqrt(z)
|
||||
* = 2f + (2c + 2s*z*R(z))
|
||||
* where f=hi part of s, and c = (z-f*f)/(s+f) is the correction term
|
||||
* for f so that f+c ~ sqrt(z).
|
||||
* For x<-0.5
|
||||
* acos(x) = pi - 2asin(sqrt((1-|x|)/2))
|
||||
* = pi - 0.5*(s+s*z*R(z)), where z=(1-|x|)/2,s=sqrt(z)
|
||||
*
|
||||
* Special cases:
|
||||
* if x is NaN, return x itself;
|
||||
* if |x|>1, return NaN with invalid signal.
|
||||
*
|
||||
* Function needed: sqrt
|
||||
*/
|
||||
|
||||
use super::sqrt;
|
||||
|
||||
const PIO2_HI: f64 = 1.57079632679489655800e+00; /* 0x3FF921FB, 0x54442D18 */
|
||||
const PIO2_LO: f64 = 6.12323399573676603587e-17; /* 0x3C91A626, 0x33145C07 */
|
||||
const PS0: f64 = 1.66666666666666657415e-01; /* 0x3FC55555, 0x55555555 */
|
||||
const PS1: f64 = -3.25565818622400915405e-01; /* 0xBFD4D612, 0x03EB6F7D */
|
||||
const PS2: f64 = 2.01212532134862925881e-01; /* 0x3FC9C155, 0x0E884455 */
|
||||
const PS3: f64 = -4.00555345006794114027e-02; /* 0xBFA48228, 0xB5688F3B */
|
||||
const PS4: f64 = 7.91534994289814532176e-04; /* 0x3F49EFE0, 0x7501B288 */
|
||||
const PS5: f64 = 3.47933107596021167570e-05; /* 0x3F023DE1, 0x0DFDF709 */
|
||||
const QS1: f64 = -2.40339491173441421878e+00; /* 0xC0033A27, 0x1C8A2D4B */
|
||||
const QS2: f64 = 2.02094576023350569471e+00; /* 0x40002AE5, 0x9C598AC8 */
|
||||
const QS3: f64 = -6.88283971605453293030e-01; /* 0xBFE6066C, 0x1B8D0159 */
|
||||
const QS4: f64 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
|
||||
|
||||
#[inline]
|
||||
fn r(z: f64) -> f64 {
|
||||
let p: f64 = z * (PS0 + z * (PS1 + z * (PS2 + z * (PS3 + z * (PS4 + z * PS5)))));
|
||||
let q: f64 = 1.0 + z * (QS1 + z * (QS2 + z * (QS3 + z * QS4)));
|
||||
return p / q;
|
||||
}
|
||||
|
||||
#[inline]
|
||||
pub fn acos(x: f64) -> f64 {
|
||||
let x1p_120f = f64::from_bits(0x3870000000000000); // 0x1p-120 === 2 ^ -120
|
||||
let z: f64;
|
||||
let w: f64;
|
||||
let s: f64;
|
||||
let c: f64;
|
||||
let df: f64;
|
||||
let hx: u32;
|
||||
let ix: u32;
|
||||
|
||||
hx = (x.to_bits() >> 32) as u32;
|
||||
ix = hx & 0x7fffffff;
|
||||
/* |x| >= 1 or nan */
|
||||
if ix >= 0x3ff00000 {
|
||||
let lx: u32 = x.to_bits() as u32;
|
||||
|
||||
if (ix - 0x3ff00000 | lx) == 0 {
|
||||
/* acos(1)=0, acos(-1)=pi */
|
||||
if (hx >> 31) != 0 {
|
||||
return 2. * PIO2_HI + x1p_120f;
|
||||
}
|
||||
return 0.;
|
||||
}
|
||||
return 0. / (x - x);
|
||||
}
|
||||
/* |x| < 0.5 */
|
||||
if ix < 0x3fe00000 {
|
||||
if ix <= 0x3c600000 {
|
||||
/* |x| < 2**-57 */
|
||||
return PIO2_HI + x1p_120f;
|
||||
}
|
||||
return PIO2_HI - (x - (PIO2_LO - x * r(x * x)));
|
||||
}
|
||||
/* x < -0.5 */
|
||||
if (hx >> 31) != 0 {
|
||||
z = (1.0 + x) * 0.5;
|
||||
s = sqrt(z);
|
||||
w = r(z) * s - PIO2_LO;
|
||||
return 2. * (PIO2_HI - (s + w));
|
||||
}
|
||||
/* x > 0.5 */
|
||||
z = (1.0 - x) * 0.5;
|
||||
s = sqrt(z);
|
||||
// Set the low 4 bytes to zero
|
||||
df = f64::from_bits(s.to_bits() & 0xff_ff_ff_ff_00_00_00_00);
|
||||
|
||||
c = (z - df * df) / (s + df);
|
||||
w = r(z) * s + c;
|
||||
return 2. * (df + w);
|
||||
}
|
||||
110
library/compiler-builtins/libm/src/math/cbrt.rs
Normal file
110
library/compiler-builtins/libm/src/math/cbrt.rs
Normal file
|
|
@ -0,0 +1,110 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/s_cbrt.c */
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*
|
||||
* Optimized by Bruce D. Evans.
|
||||
*/
|
||||
/* cbrt(x)
|
||||
* Return cube root of x
|
||||
*/
|
||||
|
||||
use core::f64;
|
||||
|
||||
const B1: u32 = 715094163; /* B1 = (1023-1023/3-0.03306235651)*2**20 */
|
||||
const B2: u32 = 696219795; /* B2 = (1023-1023/3-54/3-0.03306235651)*2**20 */
|
||||
|
||||
/* |1/cbrt(x) - p(x)| < 2**-23.5 (~[-7.93e-8, 7.929e-8]). */
|
||||
const P0: f64 = 1.87595182427177009643; /* 0x3ffe03e6, 0x0f61e692 */
|
||||
const P1: f64 = -1.88497979543377169875; /* 0xbffe28e0, 0x92f02420 */
|
||||
const P2: f64 = 1.621429720105354466140; /* 0x3ff9f160, 0x4a49d6c2 */
|
||||
const P3: f64 = -0.758397934778766047437; /* 0xbfe844cb, 0xbee751d9 */
|
||||
const P4: f64 = 0.145996192886612446982; /* 0x3fc2b000, 0xd4e4edd7 */
|
||||
|
||||
#[inline]
|
||||
pub fn cbrt(x: f64) -> f64 {
|
||||
let x1p54 = f64::from_bits(0x4350000000000000); // 0x1p54 === 2 ^ 54
|
||||
|
||||
let mut ui: u64 = x.to_bits();
|
||||
let mut r: f64;
|
||||
let s: f64;
|
||||
let mut t: f64;
|
||||
let w: f64;
|
||||
let mut hx: u32 = (ui >> 32) as u32 & 0x7fffffff;
|
||||
|
||||
if hx >= 0x7ff00000 {
|
||||
/* cbrt(NaN,INF) is itself */
|
||||
return x + x;
|
||||
}
|
||||
|
||||
/*
|
||||
* Rough cbrt to 5 bits:
|
||||
* cbrt(2**e*(1+m) ~= 2**(e/3)*(1+(e%3+m)/3)
|
||||
* where e is integral and >= 0, m is real and in [0, 1), and "/" and
|
||||
* "%" are integer division and modulus with rounding towards minus
|
||||
* infinity. The RHS is always >= the LHS and has a maximum relative
|
||||
* error of about 1 in 16. Adding a bias of -0.03306235651 to the
|
||||
* (e%3+m)/3 term reduces the error to about 1 in 32. With the IEEE
|
||||
* floating point representation, for finite positive normal values,
|
||||
* ordinary integer divison of the value in bits magically gives
|
||||
* almost exactly the RHS of the above provided we first subtract the
|
||||
* exponent bias (1023 for doubles) and later add it back. We do the
|
||||
* subtraction virtually to keep e >= 0 so that ordinary integer
|
||||
* division rounds towards minus infinity; this is also efficient.
|
||||
*/
|
||||
if hx < 0x00100000 {
|
||||
/* zero or subnormal? */
|
||||
ui = (x * x1p54).to_bits();
|
||||
hx = (ui >> 32) as u32 & 0x7fffffff;
|
||||
if hx == 0 {
|
||||
return x; /* cbrt(0) is itself */
|
||||
}
|
||||
hx = hx / 3 + B2;
|
||||
} else {
|
||||
hx = hx / 3 + B1;
|
||||
}
|
||||
ui &= 1 << 63;
|
||||
ui |= (hx as u64) << 32;
|
||||
t = f64::from_bits(ui);
|
||||
|
||||
/*
|
||||
* New cbrt to 23 bits:
|
||||
* cbrt(x) = t*cbrt(x/t**3) ~= t*P(t**3/x)
|
||||
* where P(r) is a polynomial of degree 4 that approximates 1/cbrt(r)
|
||||
* to within 2**-23.5 when |r - 1| < 1/10. The rough approximation
|
||||
* has produced t such than |t/cbrt(x) - 1| ~< 1/32, and cubing this
|
||||
* gives us bounds for r = t**3/x.
|
||||
*
|
||||
* Try to optimize for parallel evaluation as in __tanf.c.
|
||||
*/
|
||||
r = (t * t) * (t / x);
|
||||
t = t * ((P0 + r * (P1 + r * P2)) + ((r * r) * r) * (P3 + r * P4));
|
||||
|
||||
/*
|
||||
* Round t away from zero to 23 bits (sloppily except for ensuring that
|
||||
* the result is larger in magnitude than cbrt(x) but not much more than
|
||||
* 2 23-bit ulps larger). With rounding towards zero, the error bound
|
||||
* would be ~5/6 instead of ~4/6. With a maximum error of 2 23-bit ulps
|
||||
* in the rounded t, the infinite-precision error in the Newton
|
||||
* approximation barely affects third digit in the final error
|
||||
* 0.667; the error in the rounded t can be up to about 3 23-bit ulps
|
||||
* before the final error is larger than 0.667 ulps.
|
||||
*/
|
||||
ui = t.to_bits();
|
||||
ui = (ui + 0x80000000) & 0xffffffffc0000000;
|
||||
t = f64::from_bits(ui);
|
||||
|
||||
/* one step Newton iteration to 53 bits with error < 0.667 ulps */
|
||||
s = t * t; /* t*t is exact */
|
||||
r = x / s; /* error <= 0.5 ulps; |r| < |t| */
|
||||
w = t + t; /* t+t is exact */
|
||||
r = (r - t) / (w + r); /* r-t is exact; w+r ~= 3*t */
|
||||
t = t + t * r; /* error <= 0.5 + 0.5/3 + epsilon */
|
||||
t
|
||||
}
|
||||
72
library/compiler-builtins/libm/src/math/cbrtf.rs
Normal file
72
library/compiler-builtins/libm/src/math/cbrtf.rs
Normal file
|
|
@ -0,0 +1,72 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/s_cbrtf.c */
|
||||
/*
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
||||
* Debugged and optimized by Bruce D. Evans.
|
||||
*/
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
/* cbrtf(x)
|
||||
* Return cube root of x
|
||||
*/
|
||||
|
||||
use core::f32;
|
||||
|
||||
const B1: u32 = 709958130; /* B1 = (127-127.0/3-0.03306235651)*2**23 */
|
||||
const B2: u32 = 642849266; /* B2 = (127-127.0/3-24/3-0.03306235651)*2**23 */
|
||||
|
||||
#[inline]
|
||||
pub fn cbrtf(x: f32) -> f32 {
|
||||
let x1p24 = f32::from_bits(0x4b800000); // 0x1p24f === 2 ^ 24
|
||||
|
||||
let mut r: f64;
|
||||
let mut t: f64;
|
||||
let mut ui: u32 = x.to_bits();
|
||||
let mut hx: u32 = ui & 0x7fffffff;
|
||||
|
||||
if hx >= 0x7f800000 {
|
||||
/* cbrt(NaN,INF) is itself */
|
||||
return x + x;
|
||||
}
|
||||
|
||||
/* rough cbrt to 5 bits */
|
||||
if hx < 0x00800000 {
|
||||
/* zero or subnormal? */
|
||||
if hx == 0 {
|
||||
return x; /* cbrt(+-0) is itself */
|
||||
}
|
||||
ui = (x * x1p24).to_bits();
|
||||
hx = ui & 0x7fffffff;
|
||||
hx = hx / 3 + B2;
|
||||
} else {
|
||||
hx = hx / 3 + B1;
|
||||
}
|
||||
ui &= 0x80000000;
|
||||
ui |= hx;
|
||||
|
||||
/*
|
||||
* First step Newton iteration (solving t*t-x/t == 0) to 16 bits. In
|
||||
* double precision so that its terms can be arranged for efficiency
|
||||
* without causing overflow or underflow.
|
||||
*/
|
||||
t = f32::from_bits(ui) as f64;
|
||||
r = t * t * t;
|
||||
t = t * (x as f64 + x as f64 + r) / (x as f64 + r + r);
|
||||
|
||||
/*
|
||||
* Second step Newton iteration to 47 bits. In double precision for
|
||||
* efficiency and accuracy.
|
||||
*/
|
||||
r = t * t * t;
|
||||
t = t * (x as f64 + x as f64 + r) / (x as f64 + r + r);
|
||||
|
||||
/* rounding to 24 bits is perfect in round-to-nearest mode */
|
||||
t as f32
|
||||
}
|
||||
383
library/compiler-builtins/libm/src/math/exp2.rs
Normal file
383
library/compiler-builtins/libm/src/math/exp2.rs
Normal file
|
|
@ -0,0 +1,383 @@
|
|||
// origin: FreeBSD /usr/src/lib/msun/src/s_exp2.c */
|
||||
//-
|
||||
// Copyright (c) 2005 David Schultz <das@FreeBSD.ORG>
|
||||
// All rights reserved.
|
||||
//
|
||||
// Redistribution and use in source and binary forms, with or without
|
||||
// modification, are permitted provided that the following conditions
|
||||
// are met:
|
||||
// 1. Redistributions of source code must retain the above copyright
|
||||
// notice, this list of conditions and the following disclaimer.
|
||||
// 2. Redistributions in binary form must reproduce the above copyright
|
||||
// notice, this list of conditions and the following disclaimer in the
|
||||
// documentation and/or other materials provided with the distribution.
|
||||
//
|
||||
// THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
|
||||
// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
|
||||
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
|
||||
// ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
|
||||
// FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
|
||||
// DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
|
||||
// OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
|
||||
// HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
|
||||
// LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
|
||||
// OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
|
||||
// SUCH DAMAGE.
|
||||
|
||||
use super::scalbn::scalbn;
|
||||
|
||||
const TBLSIZE: usize = 256;
|
||||
|
||||
#[cfg_attr(rustfmt, rustfmt_skip)]
|
||||
static TBL: [u64; TBLSIZE * 2] = [
|
||||
// exp2(z + eps) eps
|
||||
0x3fe6a09e667f3d5d, 0x3d39880000000000,
|
||||
0x3fe6b052fa751744, 0x3cd8000000000000,
|
||||
0x3fe6c012750bd9fe, 0xbd28780000000000,
|
||||
0x3fe6cfdcddd476bf, 0x3d1ec00000000000,
|
||||
0x3fe6dfb23c651a29, 0xbcd8000000000000,
|
||||
0x3fe6ef9298593ae3, 0xbcbc000000000000,
|
||||
0x3fe6ff7df9519386, 0xbd2fd80000000000,
|
||||
0x3fe70f7466f42da3, 0xbd2c880000000000,
|
||||
0x3fe71f75e8ec5fc3, 0x3d13c00000000000,
|
||||
0x3fe72f8286eacf05, 0xbd38300000000000,
|
||||
0x3fe73f9a48a58152, 0xbd00c00000000000,
|
||||
0x3fe74fbd35d7ccfc, 0x3d2f880000000000,
|
||||
0x3fe75feb564267f1, 0x3d03e00000000000,
|
||||
0x3fe77024b1ab6d48, 0xbd27d00000000000,
|
||||
0x3fe780694fde5d38, 0xbcdd000000000000,
|
||||
0x3fe790b938ac1d00, 0x3ce3000000000000,
|
||||
0x3fe7a11473eb0178, 0xbced000000000000,
|
||||
0x3fe7b17b0976d060, 0x3d20400000000000,
|
||||
0x3fe7c1ed0130c133, 0x3ca0000000000000,
|
||||
0x3fe7d26a62ff8636, 0xbd26900000000000,
|
||||
0x3fe7e2f336cf4e3b, 0xbd02e00000000000,
|
||||
0x3fe7f3878491c3e8, 0xbd24580000000000,
|
||||
0x3fe80427543e1b4e, 0x3d33000000000000,
|
||||
0x3fe814d2add1071a, 0x3d0f000000000000,
|
||||
0x3fe82589994ccd7e, 0xbd21c00000000000,
|
||||
0x3fe8364c1eb942d0, 0x3d29d00000000000,
|
||||
0x3fe8471a4623cab5, 0x3d47100000000000,
|
||||
0x3fe857f4179f5bbc, 0x3d22600000000000,
|
||||
0x3fe868d99b4491af, 0xbd32c40000000000,
|
||||
0x3fe879cad931a395, 0xbd23000000000000,
|
||||
0x3fe88ac7d98a65b8, 0xbd2a800000000000,
|
||||
0x3fe89bd0a4785800, 0xbced000000000000,
|
||||
0x3fe8ace5422aa223, 0x3d33280000000000,
|
||||
0x3fe8be05bad619fa, 0x3d42b40000000000,
|
||||
0x3fe8cf3216b54383, 0xbd2ed00000000000,
|
||||
0x3fe8e06a5e08664c, 0xbd20500000000000,
|
||||
0x3fe8f1ae99157807, 0x3d28280000000000,
|
||||
0x3fe902fed0282c0e, 0xbd1cb00000000000,
|
||||
0x3fe9145b0b91ff96, 0xbd05e00000000000,
|
||||
0x3fe925c353aa2ff9, 0x3cf5400000000000,
|
||||
0x3fe93737b0cdc64a, 0x3d17200000000000,
|
||||
0x3fe948b82b5f98ae, 0xbd09000000000000,
|
||||
0x3fe95a44cbc852cb, 0x3d25680000000000,
|
||||
0x3fe96bdd9a766f21, 0xbd36d00000000000,
|
||||
0x3fe97d829fde4e2a, 0xbd01000000000000,
|
||||
0x3fe98f33e47a23a3, 0x3d2d000000000000,
|
||||
0x3fe9a0f170ca0604, 0xbd38a40000000000,
|
||||
0x3fe9b2bb4d53ff89, 0x3d355c0000000000,
|
||||
0x3fe9c49182a3f15b, 0x3d26b80000000000,
|
||||
0x3fe9d674194bb8c5, 0xbcec000000000000,
|
||||
0x3fe9e86319e3238e, 0x3d17d00000000000,
|
||||
0x3fe9fa5e8d07f302, 0x3d16400000000000,
|
||||
0x3fea0c667b5de54d, 0xbcf5000000000000,
|
||||
0x3fea1e7aed8eb8f6, 0x3d09e00000000000,
|
||||
0x3fea309bec4a2e27, 0x3d2ad80000000000,
|
||||
0x3fea42c980460a5d, 0xbd1af00000000000,
|
||||
0x3fea5503b23e259b, 0x3d0b600000000000,
|
||||
0x3fea674a8af46213, 0x3d38880000000000,
|
||||
0x3fea799e1330b3a7, 0x3d11200000000000,
|
||||
0x3fea8bfe53c12e8d, 0x3d06c00000000000,
|
||||
0x3fea9e6b5579fcd2, 0xbd29b80000000000,
|
||||
0x3feab0e521356fb8, 0x3d2b700000000000,
|
||||
0x3feac36bbfd3f381, 0x3cd9000000000000,
|
||||
0x3fead5ff3a3c2780, 0x3ce4000000000000,
|
||||
0x3feae89f995ad2a3, 0xbd2c900000000000,
|
||||
0x3feafb4ce622f367, 0x3d16500000000000,
|
||||
0x3feb0e07298db790, 0x3d2fd40000000000,
|
||||
0x3feb20ce6c9a89a9, 0x3d12700000000000,
|
||||
0x3feb33a2b84f1a4b, 0x3d4d470000000000,
|
||||
0x3feb468415b747e7, 0xbd38380000000000,
|
||||
0x3feb59728de5593a, 0x3c98000000000000,
|
||||
0x3feb6c6e29f1c56a, 0x3d0ad00000000000,
|
||||
0x3feb7f76f2fb5e50, 0x3cde800000000000,
|
||||
0x3feb928cf22749b2, 0xbd04c00000000000,
|
||||
0x3feba5b030a10603, 0xbd0d700000000000,
|
||||
0x3febb8e0b79a6f66, 0x3d0d900000000000,
|
||||
0x3febcc1e904bc1ff, 0x3d02a00000000000,
|
||||
0x3febdf69c3f3a16f, 0xbd1f780000000000,
|
||||
0x3febf2c25bd71db8, 0xbd10a00000000000,
|
||||
0x3fec06286141b2e9, 0xbd11400000000000,
|
||||
0x3fec199bdd8552e0, 0x3d0be00000000000,
|
||||
0x3fec2d1cd9fa64ee, 0xbd09400000000000,
|
||||
0x3fec40ab5fffd02f, 0xbd0ed00000000000,
|
||||
0x3fec544778fafd15, 0x3d39660000000000,
|
||||
0x3fec67f12e57d0cb, 0xbd1a100000000000,
|
||||
0x3fec7ba88988c1b6, 0xbd58458000000000,
|
||||
0x3fec8f6d9406e733, 0xbd1a480000000000,
|
||||
0x3feca3405751c4df, 0x3ccb000000000000,
|
||||
0x3fecb720dcef9094, 0x3d01400000000000,
|
||||
0x3feccb0f2e6d1689, 0x3cf0200000000000,
|
||||
0x3fecdf0b555dc412, 0x3cf3600000000000,
|
||||
0x3fecf3155b5bab3b, 0xbd06900000000000,
|
||||
0x3fed072d4a0789bc, 0x3d09a00000000000,
|
||||
0x3fed1b532b08c8fa, 0xbd15e00000000000,
|
||||
0x3fed2f87080d8a85, 0x3d1d280000000000,
|
||||
0x3fed43c8eacaa203, 0x3d01a00000000000,
|
||||
0x3fed5818dcfba491, 0x3cdf000000000000,
|
||||
0x3fed6c76e862e6a1, 0xbd03a00000000000,
|
||||
0x3fed80e316c9834e, 0xbd0cd80000000000,
|
||||
0x3fed955d71ff6090, 0x3cf4c00000000000,
|
||||
0x3feda9e603db32ae, 0x3cff900000000000,
|
||||
0x3fedbe7cd63a8325, 0x3ce9800000000000,
|
||||
0x3fedd321f301b445, 0xbcf5200000000000,
|
||||
0x3fede7d5641c05bf, 0xbd1d700000000000,
|
||||
0x3fedfc97337b9aec, 0xbd16140000000000,
|
||||
0x3fee11676b197d5e, 0x3d0b480000000000,
|
||||
0x3fee264614f5a3e7, 0x3d40ce0000000000,
|
||||
0x3fee3b333b16ee5c, 0x3d0c680000000000,
|
||||
0x3fee502ee78b3fb4, 0xbd09300000000000,
|
||||
0x3fee653924676d68, 0xbce5000000000000,
|
||||
0x3fee7a51fbc74c44, 0xbd07f80000000000,
|
||||
0x3fee8f7977cdb726, 0xbcf3700000000000,
|
||||
0x3feea4afa2a490e8, 0x3ce5d00000000000,
|
||||
0x3feeb9f4867ccae4, 0x3d161a0000000000,
|
||||
0x3feecf482d8e680d, 0x3cf5500000000000,
|
||||
0x3feee4aaa2188514, 0x3cc6400000000000,
|
||||
0x3feefa1bee615a13, 0xbcee800000000000,
|
||||
0x3fef0f9c1cb64106, 0xbcfa880000000000,
|
||||
0x3fef252b376bb963, 0xbd2c900000000000,
|
||||
0x3fef3ac948dd7275, 0x3caa000000000000,
|
||||
0x3fef50765b6e4524, 0xbcf4f00000000000,
|
||||
0x3fef6632798844fd, 0x3cca800000000000,
|
||||
0x3fef7bfdad9cbe38, 0x3cfabc0000000000,
|
||||
0x3fef91d802243c82, 0xbcd4600000000000,
|
||||
0x3fefa7c1819e908e, 0xbd0b0c0000000000,
|
||||
0x3fefbdba3692d511, 0xbcc0e00000000000,
|
||||
0x3fefd3c22b8f7194, 0xbd10de8000000000,
|
||||
0x3fefe9d96b2a23ee, 0x3cee430000000000,
|
||||
0x3ff0000000000000, 0x0,
|
||||
0x3ff00b1afa5abcbe, 0xbcb3400000000000,
|
||||
0x3ff0163da9fb3303, 0xbd12170000000000,
|
||||
0x3ff02168143b0282, 0x3cba400000000000,
|
||||
0x3ff02c9a3e77806c, 0x3cef980000000000,
|
||||
0x3ff037d42e11bbca, 0xbcc7400000000000,
|
||||
0x3ff04315e86e7f89, 0x3cd8300000000000,
|
||||
0x3ff04e5f72f65467, 0xbd1a3f0000000000,
|
||||
0x3ff059b0d315855a, 0xbd02840000000000,
|
||||
0x3ff0650a0e3c1f95, 0x3cf1600000000000,
|
||||
0x3ff0706b29ddf71a, 0x3d15240000000000,
|
||||
0x3ff07bd42b72a82d, 0xbce9a00000000000,
|
||||
0x3ff0874518759bd0, 0x3ce6400000000000,
|
||||
0x3ff092bdf66607c8, 0xbd00780000000000,
|
||||
0x3ff09e3ecac6f383, 0xbc98000000000000,
|
||||
0x3ff0a9c79b1f3930, 0x3cffa00000000000,
|
||||
0x3ff0b5586cf988fc, 0xbcfac80000000000,
|
||||
0x3ff0c0f145e46c8a, 0x3cd9c00000000000,
|
||||
0x3ff0cc922b724816, 0x3d05200000000000,
|
||||
0x3ff0d83b23395dd8, 0xbcfad00000000000,
|
||||
0x3ff0e3ec32d3d1f3, 0x3d1bac0000000000,
|
||||
0x3ff0efa55fdfa9a6, 0xbd04e80000000000,
|
||||
0x3ff0fb66affed2f0, 0xbd0d300000000000,
|
||||
0x3ff1073028d7234b, 0x3cf1500000000000,
|
||||
0x3ff11301d0125b5b, 0x3cec000000000000,
|
||||
0x3ff11edbab5e2af9, 0x3d16bc0000000000,
|
||||
0x3ff12abdc06c31d5, 0x3ce8400000000000,
|
||||
0x3ff136a814f2047d, 0xbd0ed00000000000,
|
||||
0x3ff1429aaea92de9, 0x3ce8e00000000000,
|
||||
0x3ff14e95934f3138, 0x3ceb400000000000,
|
||||
0x3ff15a98c8a58e71, 0x3d05300000000000,
|
||||
0x3ff166a45471c3df, 0x3d03380000000000,
|
||||
0x3ff172b83c7d5211, 0x3d28d40000000000,
|
||||
0x3ff17ed48695bb9f, 0xbd05d00000000000,
|
||||
0x3ff18af9388c8d93, 0xbd1c880000000000,
|
||||
0x3ff1972658375d66, 0x3d11f00000000000,
|
||||
0x3ff1a35beb6fcba7, 0x3d10480000000000,
|
||||
0x3ff1af99f81387e3, 0xbd47390000000000,
|
||||
0x3ff1bbe084045d54, 0x3d24e40000000000,
|
||||
0x3ff1c82f95281c43, 0xbd0a200000000000,
|
||||
0x3ff1d4873168b9b2, 0x3ce3800000000000,
|
||||
0x3ff1e0e75eb44031, 0x3ceac00000000000,
|
||||
0x3ff1ed5022fcd938, 0x3d01900000000000,
|
||||
0x3ff1f9c18438cdf7, 0xbd1b780000000000,
|
||||
0x3ff2063b88628d8f, 0x3d2d940000000000,
|
||||
0x3ff212be3578a81e, 0x3cd8000000000000,
|
||||
0x3ff21f49917ddd41, 0x3d2b340000000000,
|
||||
0x3ff22bdda2791323, 0x3d19f80000000000,
|
||||
0x3ff2387a6e7561e7, 0xbd19c80000000000,
|
||||
0x3ff2451ffb821427, 0x3d02300000000000,
|
||||
0x3ff251ce4fb2a602, 0xbd13480000000000,
|
||||
0x3ff25e85711eceb0, 0x3d12700000000000,
|
||||
0x3ff26b4565e27d16, 0x3d11d00000000000,
|
||||
0x3ff2780e341de00f, 0x3d31ee0000000000,
|
||||
0x3ff284dfe1f5633e, 0xbd14c00000000000,
|
||||
0x3ff291ba7591bb30, 0xbd13d80000000000,
|
||||
0x3ff29e9df51fdf09, 0x3d08b00000000000,
|
||||
0x3ff2ab8a66d10e9b, 0xbd227c0000000000,
|
||||
0x3ff2b87fd0dada3a, 0x3d2a340000000000,
|
||||
0x3ff2c57e39771af9, 0xbd10800000000000,
|
||||
0x3ff2d285a6e402d9, 0xbd0ed00000000000,
|
||||
0x3ff2df961f641579, 0xbcf4200000000000,
|
||||
0x3ff2ecafa93e2ecf, 0xbd24980000000000,
|
||||
0x3ff2f9d24abd8822, 0xbd16300000000000,
|
||||
0x3ff306fe0a31b625, 0xbd32360000000000,
|
||||
0x3ff31432edeea50b, 0xbd70df8000000000,
|
||||
0x3ff32170fc4cd7b8, 0xbd22480000000000,
|
||||
0x3ff32eb83ba8e9a2, 0xbd25980000000000,
|
||||
0x3ff33c08b2641766, 0x3d1ed00000000000,
|
||||
0x3ff3496266e3fa27, 0xbcdc000000000000,
|
||||
0x3ff356c55f929f0f, 0xbd30d80000000000,
|
||||
0x3ff36431a2de88b9, 0x3d22c80000000000,
|
||||
0x3ff371a7373aaa39, 0x3d20600000000000,
|
||||
0x3ff37f26231e74fe, 0xbd16600000000000,
|
||||
0x3ff38cae6d05d838, 0xbd0ae00000000000,
|
||||
0x3ff39a401b713ec3, 0xbd44720000000000,
|
||||
0x3ff3a7db34e5a020, 0x3d08200000000000,
|
||||
0x3ff3b57fbfec6e95, 0x3d3e800000000000,
|
||||
0x3ff3c32dc313a8f2, 0x3cef800000000000,
|
||||
0x3ff3d0e544ede122, 0xbd17a00000000000,
|
||||
0x3ff3dea64c1234bb, 0x3d26300000000000,
|
||||
0x3ff3ec70df1c4ecc, 0xbd48a60000000000,
|
||||
0x3ff3fa4504ac7e8c, 0xbd3cdc0000000000,
|
||||
0x3ff40822c367a0bb, 0x3d25b80000000000,
|
||||
0x3ff4160a21f72e95, 0x3d1ec00000000000,
|
||||
0x3ff423fb27094646, 0xbd13600000000000,
|
||||
0x3ff431f5d950a920, 0x3d23980000000000,
|
||||
0x3ff43ffa3f84b9eb, 0x3cfa000000000000,
|
||||
0x3ff44e0860618919, 0xbcf6c00000000000,
|
||||
0x3ff45c2042a7d201, 0xbd0bc00000000000,
|
||||
0x3ff46a41ed1d0016, 0xbd12800000000000,
|
||||
0x3ff4786d668b3326, 0x3d30e00000000000,
|
||||
0x3ff486a2b5c13c00, 0xbd2d400000000000,
|
||||
0x3ff494e1e192af04, 0x3d0c200000000000,
|
||||
0x3ff4a32af0d7d372, 0xbd1e500000000000,
|
||||
0x3ff4b17dea6db801, 0x3d07800000000000,
|
||||
0x3ff4bfdad53629e1, 0xbd13800000000000,
|
||||
0x3ff4ce41b817c132, 0x3d00800000000000,
|
||||
0x3ff4dcb299fddddb, 0x3d2c700000000000,
|
||||
0x3ff4eb2d81d8ab96, 0xbd1ce00000000000,
|
||||
0x3ff4f9b2769d2d02, 0x3d19200000000000,
|
||||
0x3ff508417f4531c1, 0xbd08c00000000000,
|
||||
0x3ff516daa2cf662a, 0xbcfa000000000000,
|
||||
0x3ff5257de83f51ea, 0x3d4a080000000000,
|
||||
0x3ff5342b569d4eda, 0xbd26d80000000000,
|
||||
0x3ff542e2f4f6ac1a, 0xbd32440000000000,
|
||||
0x3ff551a4ca5d94db, 0x3d483c0000000000,
|
||||
0x3ff56070dde9116b, 0x3d24b00000000000,
|
||||
0x3ff56f4736b529de, 0x3d415a0000000000,
|
||||
0x3ff57e27dbe2c40e, 0xbd29e00000000000,
|
||||
0x3ff58d12d497c76f, 0xbd23080000000000,
|
||||
0x3ff59c0827ff0b4c, 0x3d4dec0000000000,
|
||||
0x3ff5ab07dd485427, 0xbcc4000000000000,
|
||||
0x3ff5ba11fba87af4, 0x3d30080000000000,
|
||||
0x3ff5c9268a59460b, 0xbd26c80000000000,
|
||||
0x3ff5d84590998e3f, 0x3d469a0000000000,
|
||||
0x3ff5e76f15ad20e1, 0xbd1b400000000000,
|
||||
0x3ff5f6a320dcebca, 0x3d17700000000000,
|
||||
0x3ff605e1b976dcb8, 0x3d26f80000000000,
|
||||
0x3ff6152ae6cdf715, 0x3d01000000000000,
|
||||
0x3ff6247eb03a5531, 0xbd15d00000000000,
|
||||
0x3ff633dd1d1929b5, 0xbd12d00000000000,
|
||||
0x3ff6434634ccc313, 0xbcea800000000000,
|
||||
0x3ff652b9febc8efa, 0xbd28600000000000,
|
||||
0x3ff6623882553397, 0x3d71fe0000000000,
|
||||
0x3ff671c1c708328e, 0xbd37200000000000,
|
||||
0x3ff68155d44ca97e, 0x3ce6800000000000,
|
||||
0x3ff690f4b19e9471, 0xbd29780000000000,
|
||||
];
|
||||
|
||||
// exp2(x): compute the base 2 exponential of x
|
||||
//
|
||||
// Accuracy: Peak error < 0.503 ulp for normalized results.
|
||||
//
|
||||
// Method: (accurate tables)
|
||||
//
|
||||
// Reduce x:
|
||||
// x = k + y, for integer k and |y| <= 1/2.
|
||||
// Thus we have exp2(x) = 2**k * exp2(y).
|
||||
//
|
||||
// Reduce y:
|
||||
// y = i/TBLSIZE + z - eps[i] for integer i near y * TBLSIZE.
|
||||
// Thus we have exp2(y) = exp2(i/TBLSIZE) * exp2(z - eps[i]),
|
||||
// with |z - eps[i]| <= 2**-9 + 2**-39 for the table used.
|
||||
//
|
||||
// We compute exp2(i/TBLSIZE) via table lookup and exp2(z - eps[i]) via
|
||||
// a degree-5 minimax polynomial with maximum error under 1.3 * 2**-61.
|
||||
// The values in exp2t[] and eps[] are chosen such that
|
||||
// exp2t[i] = exp2(i/TBLSIZE + eps[i]), and eps[i] is a small offset such
|
||||
// that exp2t[i] is accurate to 2**-64.
|
||||
//
|
||||
// Note that the range of i is +-TBLSIZE/2, so we actually index the tables
|
||||
// by i0 = i + TBLSIZE/2. For cache efficiency, exp2t[] and eps[] are
|
||||
// virtual tables, interleaved in the real table tbl[].
|
||||
//
|
||||
// This method is due to Gal, with many details due to Gal and Bachelis:
|
||||
//
|
||||
// Gal, S. and Bachelis, B. An Accurate Elementary Mathematical Library
|
||||
// for the IEEE Floating Point Standard. TOMS 17(1), 26-46 (1991).
|
||||
pub fn exp2(mut x: f64) -> f64 {
|
||||
let redux = f64::from_bits(0x4338000000000000) / TBLSIZE as f64;
|
||||
let p1 = f64::from_bits(0x3fe62e42fefa39ef);
|
||||
let p2 = f64::from_bits(0x3fcebfbdff82c575);
|
||||
let p3 = f64::from_bits(0x3fac6b08d704a0a6);
|
||||
let p4 = f64::from_bits(0x3f83b2ab88f70400);
|
||||
let p5 = f64::from_bits(0x3f55d88003875c74);
|
||||
|
||||
// double_t r, t, z;
|
||||
// uint32_t ix, i0;
|
||||
// union {double f; uint64_t i;} u = {x};
|
||||
// union {uint32_t u; int32_t i;} k;
|
||||
let x1p1023 = f64::from_bits(0x7fe0000000000000);
|
||||
let x1p52 = f64::from_bits(0x4330000000000000);
|
||||
let _0x1p_149 = f64::from_bits(0xb6a0000000000000);
|
||||
|
||||
/* Filter out exceptional cases. */
|
||||
let ui = f64::to_bits(x);
|
||||
let ix = ui >> 32 & 0x7fffffff;
|
||||
if ix >= 0x408ff000 {
|
||||
/* |x| >= 1022 or nan */
|
||||
if ix >= 0x40900000 && ui >> 63 == 0 {
|
||||
/* x >= 1024 or nan */
|
||||
/* overflow */
|
||||
x *= x1p1023;
|
||||
return x;
|
||||
}
|
||||
if ix >= 0x7ff00000 {
|
||||
/* -inf or -nan */
|
||||
return -1.0 / x;
|
||||
}
|
||||
if ui >> 63 != 0 {
|
||||
/* x <= -1022 */
|
||||
/* underflow */
|
||||
if x <= -1075.0 || x - x1p52 + x1p52 != x {
|
||||
force_eval!((_0x1p_149 / x) as f32);
|
||||
}
|
||||
if x <= -1075.0 {
|
||||
return 0.0;
|
||||
}
|
||||
}
|
||||
} else if ix < 0x3c900000 {
|
||||
/* |x| < 0x1p-54 */
|
||||
return 1.0 + x;
|
||||
}
|
||||
|
||||
/* Reduce x, computing z, i0, and k. */
|
||||
let ui = f64::to_bits(x + redux);
|
||||
let mut i0 = ui as u32;
|
||||
i0 += TBLSIZE as u32 / 2;
|
||||
let ku = i0 / TBLSIZE as u32 * TBLSIZE as u32;
|
||||
let ki = ku as i32 / TBLSIZE as i32;
|
||||
i0 %= TBLSIZE as u32;
|
||||
let uf = f64::from_bits(ui) - redux;
|
||||
let mut z = x - uf;
|
||||
|
||||
/* Compute r = exp2(y) = exp2t[i0] * p(z - eps[i]). */
|
||||
let t = f64::from_bits(TBL[2 * i0 as usize]); /* exp2t[i0] */
|
||||
z -= f64::from_bits(TBL[2 * i0 as usize + 1]); /* eps[i0] */
|
||||
let r = t + t * z * (p1 + z * (p2 + z * (p3 + z * (p4 + z * p5))));
|
||||
|
||||
scalbn(r, ki)
|
||||
}
|
||||
130
library/compiler-builtins/libm/src/math/exp2f.rs
Normal file
130
library/compiler-builtins/libm/src/math/exp2f.rs
Normal file
|
|
@ -0,0 +1,130 @@
|
|||
// origin: FreeBSD /usr/src/lib/msun/src/s_exp2f.c
|
||||
//-
|
||||
// Copyright (c) 2005 David Schultz <das@FreeBSD.ORG>
|
||||
// All rights reserved.
|
||||
//
|
||||
// Redistribution and use in source and binary forms, with or without
|
||||
// modification, are permitted provided that the following conditions
|
||||
// are met:
|
||||
// 1. Redistributions of source code must retain the above copyright
|
||||
// notice, this list of conditions and the following disclaimer.
|
||||
// 2. Redistributions in binary form must reproduce the above copyright
|
||||
// notice, this list of conditions and the following disclaimer in the
|
||||
// documentation and/or other materials provided with the distribution.
|
||||
//
|
||||
// THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
|
||||
// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
|
||||
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
|
||||
// ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
|
||||
// FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
|
||||
// DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
|
||||
// OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
|
||||
// HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
|
||||
// LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
|
||||
// OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
|
||||
// SUCH DAMAGE.
|
||||
|
||||
const TBLSIZE: usize = 16;
|
||||
|
||||
static EXP2FT: [u64; TBLSIZE] = [
|
||||
0x3fe6a09e667f3bcd,
|
||||
0x3fe7a11473eb0187,
|
||||
0x3fe8ace5422aa0db,
|
||||
0x3fe9c49182a3f090,
|
||||
0x3feae89f995ad3ad,
|
||||
0x3fec199bdd85529c,
|
||||
0x3fed5818dcfba487,
|
||||
0x3feea4afa2a490da,
|
||||
0x3ff0000000000000,
|
||||
0x3ff0b5586cf9890f,
|
||||
0x3ff172b83c7d517b,
|
||||
0x3ff2387a6e756238,
|
||||
0x3ff306fe0a31b715,
|
||||
0x3ff3dea64c123422,
|
||||
0x3ff4bfdad5362a27,
|
||||
0x3ff5ab07dd485429,
|
||||
];
|
||||
|
||||
// exp2f(x): compute the base 2 exponential of x
|
||||
//
|
||||
// Accuracy: Peak error < 0.501 ulp; location of peak: -0.030110927.
|
||||
//
|
||||
// Method: (equally-spaced tables)
|
||||
//
|
||||
// Reduce x:
|
||||
// x = k + y, for integer k and |y| <= 1/2.
|
||||
// Thus we have exp2f(x) = 2**k * exp2(y).
|
||||
//
|
||||
// Reduce y:
|
||||
// y = i/TBLSIZE + z for integer i near y * TBLSIZE.
|
||||
// Thus we have exp2(y) = exp2(i/TBLSIZE) * exp2(z),
|
||||
// with |z| <= 2**-(TBLSIZE+1).
|
||||
//
|
||||
// We compute exp2(i/TBLSIZE) via table lookup and exp2(z) via a
|
||||
// degree-4 minimax polynomial with maximum error under 1.4 * 2**-33.
|
||||
// Using double precision for everything except the reduction makes
|
||||
// roundoff error insignificant and simplifies the scaling step.
|
||||
//
|
||||
// This method is due to Tang, but I do not use his suggested parameters:
|
||||
//
|
||||
// Tang, P. Table-driven Implementation of the Exponential Function
|
||||
// in IEEE Floating-Point Arithmetic. TOMS 15(2), 144-157 (1989).
|
||||
pub fn exp2f(mut x: f32) -> f32 {
|
||||
let redux = f32::from_bits(0x4b400000) / TBLSIZE as f32;
|
||||
let p1 = f32::from_bits(0x3f317218);
|
||||
let p2 = f32::from_bits(0x3e75fdf0);
|
||||
let p3 = f32::from_bits(0x3d6359a4);
|
||||
let p4 = f32::from_bits(0x3c1d964e);
|
||||
|
||||
// double_t t, r, z;
|
||||
// uint32_t ix, i0, k;
|
||||
|
||||
let x1p127 = f32::from_bits(0x7f000000);
|
||||
|
||||
/* Filter out exceptional cases. */
|
||||
let ui = f32::to_bits(x);
|
||||
let ix = ui & 0x7fffffff;
|
||||
if ix > 0x42fc0000 {
|
||||
/* |x| > 126 */
|
||||
if ix > 0x7f800000 {
|
||||
/* NaN */
|
||||
return x;
|
||||
}
|
||||
if ui >= 0x43000000 && ui < 0x80000000 {
|
||||
/* x >= 128 */
|
||||
x *= x1p127;
|
||||
return x;
|
||||
}
|
||||
if ui >= 0x80000000 {
|
||||
/* x < -126 */
|
||||
if ui >= 0xc3160000 || (ui & 0x0000ffff != 0) {
|
||||
force_eval!(f32::from_bits(0x80000001) / x);
|
||||
}
|
||||
if ui >= 0xc3160000 {
|
||||
/* x <= -150 */
|
||||
return 0.0;
|
||||
}
|
||||
}
|
||||
} else if ix <= 0x33000000 {
|
||||
/* |x| <= 0x1p-25 */
|
||||
return 1.0 + x;
|
||||
}
|
||||
|
||||
/* Reduce x, computing z, i0, and k. */
|
||||
let ui = f32::to_bits(x + redux);
|
||||
let mut i0 = ui;
|
||||
i0 += TBLSIZE as u32 / 2;
|
||||
let k = i0 / TBLSIZE as u32;
|
||||
let ukf = f64::from_bits(((0x3ff + k) as u64) << 52);
|
||||
i0 &= TBLSIZE as u32 - 1;
|
||||
let mut uf = f32::from_bits(ui);
|
||||
uf -= redux;
|
||||
let z: f64 = (x - uf) as f64;
|
||||
/* Compute r = exp2(y) = exp2ft[i0] * p(z). */
|
||||
let r: f64 = f64::from_bits(EXP2FT[i0 as usize]);
|
||||
let t: f64 = r as f64 * z;
|
||||
let r: f64 = r + t * (p1 as f64 + z * p2 as f64) + t * (z * z) * (p3 as f64 + z * p4 as f64);
|
||||
|
||||
/* Scale by 2**k */
|
||||
(r * ukf) as f32
|
||||
}
|
||||
|
|
@ -11,6 +11,7 @@ const Q3: f64 = -7.93650757867487942473e-05; /* BF14CE19 9EAADBB7 */
|
|||
const Q4: f64 = 4.00821782732936239552e-06; /* 3ED0CFCA 86E65239 */
|
||||
const Q5: f64 = -2.01099218183624371326e-07; /* BE8AFDB7 6E09C32D */
|
||||
|
||||
#[allow(warnings)]
|
||||
pub fn expm1(mut x: f64) -> f64 {
|
||||
let hi: f64;
|
||||
let lo: f64;
|
||||
|
|
@ -19,8 +20,8 @@ pub fn expm1(mut x: f64) -> f64 {
|
|||
let mut t: f64;
|
||||
let mut y: f64;
|
||||
|
||||
let mut ui = x.to_bits() >> 32;
|
||||
let hx = ui & 0x7fffffff;
|
||||
let mut ui = x.to_bits();
|
||||
let hx = ((ui >> 32) & 0x7fffffff) as u32;
|
||||
let sign = (ui >> 63) as i32;
|
||||
|
||||
/* filter out huge and non-finite argument */
|
||||
|
|
@ -63,7 +64,7 @@ pub fn expm1(mut x: f64) -> f64 {
|
|||
} else if hx < 0x3c900000 {
|
||||
/* |x| < 2**-54, return x */
|
||||
if hx < 0x00100000 {
|
||||
force_eval!(x as f32);
|
||||
force_eval!(x);
|
||||
}
|
||||
return x;
|
||||
} else {
|
||||
|
|
|
|||
80
library/compiler-builtins/libm/src/math/fmod.rs
Normal file
80
library/compiler-builtins/libm/src/math/fmod.rs
Normal file
|
|
@ -0,0 +1,80 @@
|
|||
use core::u64;
|
||||
|
||||
#[inline]
|
||||
pub fn fmod(x: f64, y: f64) -> f64 {
|
||||
let mut uxi = x.to_bits();
|
||||
let mut uyi = y.to_bits();
|
||||
let mut ex = (uxi >> 52 & 0x7ff) as i64;
|
||||
let mut ey = (uyi >> 52 & 0x7ff) as i64;
|
||||
let sx = uxi >> 63;
|
||||
let mut i;
|
||||
|
||||
if uyi << 1 == 0 || y.is_nan() || ex == 0x7ff {
|
||||
return (x * y) / (x * y);
|
||||
}
|
||||
if uxi << 1 <= uyi << 1 {
|
||||
if uxi << 1 == uyi << 1 {
|
||||
return 0.0 * x;
|
||||
}
|
||||
return x;
|
||||
}
|
||||
|
||||
/* normalize x and y */
|
||||
if ex == 0 {
|
||||
i = uxi << 12;
|
||||
while i >> 63 == 0 {
|
||||
ex -= 1;
|
||||
i <<= 1;
|
||||
}
|
||||
uxi <<= -ex + 1;
|
||||
} else {
|
||||
uxi &= u64::MAX >> 12;
|
||||
uxi |= 1 << 52;
|
||||
}
|
||||
if ey == 0 {
|
||||
i = uyi << 12;
|
||||
while i >> 63 == 0 {
|
||||
ey -= 1;
|
||||
i <<= 1;
|
||||
}
|
||||
uyi <<= -ey + 1;
|
||||
} else {
|
||||
uyi &= u64::MAX >> 12;
|
||||
uyi |= 1 << 52;
|
||||
}
|
||||
|
||||
/* x mod y */
|
||||
while ex > ey {
|
||||
i = uxi - uyi;
|
||||
if i >> 63 == 0 {
|
||||
if i == 0 {
|
||||
return 0.0 * x;
|
||||
}
|
||||
uxi = i;
|
||||
}
|
||||
uxi <<= 1;
|
||||
ex -= 1;
|
||||
}
|
||||
i = uxi - uyi;
|
||||
if i >> 63 == 0 {
|
||||
if i == 0 {
|
||||
return 0.0 * x;
|
||||
}
|
||||
uxi = i;
|
||||
}
|
||||
while uxi >> 52 == 0 {
|
||||
uxi <<= 1;
|
||||
ex -= 1;
|
||||
}
|
||||
|
||||
/* scale result */
|
||||
if ex > 0 {
|
||||
uxi -= 1 << 52;
|
||||
uxi |= (ex as u64) << 52;
|
||||
} else {
|
||||
uxi >>= -ex + 1;
|
||||
}
|
||||
uxi |= (sx as u64) << 63;
|
||||
|
||||
f64::from_bits(uxi)
|
||||
}
|
||||
|
|
@ -6,10 +6,15 @@ macro_rules! force_eval {
|
|||
};
|
||||
}
|
||||
|
||||
mod acos;
|
||||
mod cbrt;
|
||||
mod cbrtf;
|
||||
mod ceil;
|
||||
mod ceilf;
|
||||
mod cosf;
|
||||
mod exp;
|
||||
mod exp2;
|
||||
mod exp2f;
|
||||
mod expf;
|
||||
mod expm1;
|
||||
mod fabs;
|
||||
|
|
@ -18,6 +23,7 @@ mod fdim;
|
|||
mod fdimf;
|
||||
mod floor;
|
||||
mod floorf;
|
||||
mod fmod;
|
||||
mod fmodf;
|
||||
mod hypot;
|
||||
mod hypotf;
|
||||
|
|
@ -40,10 +46,15 @@ mod trunc;
|
|||
mod truncf;
|
||||
|
||||
// Use separated imports instead of {}-grouped imports for easier merging.
|
||||
pub use self::acos::acos;
|
||||
pub use self::cbrt::cbrt;
|
||||
pub use self::cbrtf::cbrtf;
|
||||
pub use self::ceil::ceil;
|
||||
pub use self::ceilf::ceilf;
|
||||
pub use self::cosf::cosf;
|
||||
pub use self::exp::exp;
|
||||
pub use self::exp2::exp2;
|
||||
pub use self::exp2f::exp2f;
|
||||
pub use self::expf::expf;
|
||||
pub use self::expm1::expm1;
|
||||
pub use self::fabs::fabs;
|
||||
|
|
@ -52,6 +63,7 @@ pub use self::fdim::fdim;
|
|||
pub use self::fdimf::fdimf;
|
||||
pub use self::floor::floor;
|
||||
pub use self::floorf::floorf;
|
||||
pub use self::fmod::fmod;
|
||||
pub use self::fmodf::fmodf;
|
||||
pub use self::hypot::hypot;
|
||||
pub use self::hypotf::hypotf;
|
||||
|
|
|
|||
|
|
@ -656,11 +656,11 @@ f32_f32! {
|
|||
truncf,
|
||||
// asinf,
|
||||
// atanf,
|
||||
// cbrtf,
|
||||
cbrtf,
|
||||
cosf,
|
||||
ceilf,
|
||||
// coshf,
|
||||
// exp2f,
|
||||
exp2f,
|
||||
expf,
|
||||
log10f,
|
||||
log1pf,
|
||||
|
|
@ -696,15 +696,15 @@ f32i32_f32! {
|
|||
|
||||
// With signature `fn(f64) -> f64`
|
||||
f64_f64! {
|
||||
// acos,
|
||||
acos,
|
||||
// asin,
|
||||
// atan,
|
||||
// cbrt,
|
||||
cbrt,
|
||||
ceil,
|
||||
// cos,
|
||||
// cosh,
|
||||
exp,
|
||||
// exp2,
|
||||
exp2,
|
||||
expm1,
|
||||
floor,
|
||||
log,
|
||||
|
|
@ -725,7 +725,7 @@ f64_f64! {
|
|||
f64f64_f64! {
|
||||
// atan2,
|
||||
fdim,
|
||||
// fmod,
|
||||
fmod,
|
||||
hypot,
|
||||
// pow,
|
||||
}
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue