rust-lang/rust
1
0
Fork 0
Code Issues Pull Requests Actions 1 Packages Projects Releases Wiki Activity

Merge branch 'master' into master

Browse Source
This commit is contained in:
Jorge Aparicio
2018-07-14 11:30:47 -05:00
committed by GitHub
parent 88faf19ca5 34ae08d2a2
commit ebb89e686d
20 changed files with 1199 additions and 20 deletions
Download Patch File Download Diff File Expand all files Collapse all files

10
library/compiler-builtins/libm/src/lib.rs
View File

@@ -87,7 +87,6 @@ pub trait F32Ext: private::Sealed {
#[cfg(todo)]
fn sin(self) -> Self;
#[cfg(todo)]
fn cos(self) -> Self;
#[cfg(todo)]
@@ -114,7 +113,6 @@ pub trait F32Ext: private::Sealed {
#[cfg(todo)]
fn exp_m1(self) -> Self;
#[cfg(todo)]
fn ln_1p(self) -> Self;
#[cfg(todo)]
@@ -253,7 +251,6 @@ impl F32Ext for f32 {
sinf(self)
}
#[cfg(todo)]
#[inline]
fn cos(self) -> Self {
cosf(self)
@@ -295,7 +292,6 @@ impl F32Ext for f32 {
expm1f(self)
}
#[cfg(todo)]
#[inline]
fn ln_1p(self) -> Self {
log1pf(self)
@@ -389,10 +385,8 @@ pub trait F64Ext: private::Sealed {
#[cfg(todo)]
fn exp2(self) -> Self;
#[cfg(todo)]
fn ln(self) -> Self;
#[cfg(todo)]
fn log(self, base: Self) -> Self;
fn log2(self) -> Self;
@@ -434,7 +428,6 @@ pub trait F64Ext: private::Sealed {
#[cfg(todo)]
fn exp_m1(self) -> Self;
#[cfg(todo)]
fn ln_1p(self) -> Self;
#[cfg(todo)]
@@ -539,13 +532,11 @@ impl F64Ext for f64 {
exp2(self)
}
#[cfg(todo)]
#[inline]
fn ln(self) -> Self {
log(self)
}
#[cfg(todo)]
#[inline]
fn log(self, base: Self) -> Self {
self.ln() / base.ln()
@@ -620,7 +611,6 @@ impl F64Ext for f64 {
expm1(self)
}
#[cfg(todo)]
#[inline]
fn ln_1p(self) -> Self {
log1p(self)

71
library/compiler-builtins/libm/src/math/cosf.rs Normal file
View File

@@ -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),
}
}

15
library/compiler-builtins/libm/src/math/expf.rs
View File

@@ -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];

13
library/compiler-builtins/libm/src/math/k_cosf.rs Normal file
View File

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

14
library/compiler-builtins/libm/src/math/k_sinf.rs Normal file
View File

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

117
library/compiler-builtins/libm/src/math/log.rs Normal file
View File

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

19
library/compiler-builtins/libm/src/math/log10.rs
View File

@@ -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 */

15
library/compiler-builtins/libm/src/math/log10f.rs
View File

@@ -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 */

142
library/compiler-builtins/libm/src/math/log1p.rs Normal file
View File

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

97
library/compiler-builtins/libm/src/math/log1pf.rs Normal file
View File

@@ -0,0 +1,97 @@
/* origin: FreeBSD /usr/src/lib/msun/src/s_log1pf.c */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
use core::f32;
const LN2_HI: f32 = 6.9313812256e-01; /* 0x3f317180 */
const LN2_LO: f32 = 9.0580006145e-06; /* 0x3717f7d1 */
/* |(log(1+s)-log(1-s))/s - Lg(s)| < 2**-34.24 (~[-4.95e-11, 4.97e-11]). */
const LG1: f32 = 0.66666662693; /* 0xaaaaaa.0p-24 */
const LG2: f32 = 0.40000972152; /* 0xccce13.0p-25 */
const LG3: f32 = 0.28498786688; /* 0x91e9ee.0p-25 */
const LG4: f32 = 0.24279078841; /* 0xf89e26.0p-26 */
pub fn log1pf(x: f32) -> f32 {
let mut ui: u32 = x.to_bits();
let hfsq: f32;
let mut f: f32 = 0.;
let mut c: f32 = 0.;
let s: f32;
let z: f32;
let r: f32;
let w: f32;
let t1: f32;
let t2: f32;
let dk: f32;
let ix: u32;
let mut iu: u32;
let mut k: i32;
ix = ui;
k = 1;
if ix < 0x3ed413d0 || (ix >> 31) > 0 {
/* 1+x < sqrt(2)+ */
if ix >= 0xbf800000 {
/* x <= -1.0 */
if x == -1. {
return x / 0.0; /* log1p(-1)=+inf */
}
return (x - x) / 0.0; /* log1p(x<-1)=NaN */
}
if ix << 1 < 0x33800000 << 1 {
/* |x| < 2**-24 */
/* underflow if subnormal */
if (ix & 0x7f800000) == 0 {
force_eval!(x * x);
}
return x;
}
if ix <= 0xbe95f619 {
/* sqrt(2)/2- <= 1+x < sqrt(2)+ */
k = 0;
c = 0.;
f = x;
}
} else if ix >= 0x7f800000 {
return x;
}
if k > 0 {
ui = (1. + x).to_bits();
iu = ui;
iu += 0x3f800000 - 0x3f3504f3;
k = (iu >> 23) as i32 - 0x7f;
/* correction term ~ log(1+x)-log(u), avoid underflow in c/u */
if k < 25 {
c = if k >= 2 {
1. - (f32::from_bits(ui) - x)
} else {
x - (f32::from_bits(ui) - 1.)
};
c /= f32::from_bits(ui);
} else {
c = 0.;
}
/* reduce u into [sqrt(2)/2, sqrt(2)] */
iu = (iu & 0x007fffff) + 0x3f3504f3;
ui = iu;
f = f32::from_bits(ui) - 1.;
}
s = f / (2.0 + f);
z = s * s;
w = z * z;
t1 = w * (LG2 + w * LG4);
t2 = z * (LG1 + w * LG3);
r = t2 + t1;
hfsq = 0.5 * f * f;
dk = k as f32;
return s * (hfsq + r) + (dk * LN2_LO + c) - hfsq + f + dk * LN2_HI;
}

19
library/compiler-builtins/libm/src/math/log2.rs
View File

@@ -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 */

15
library/compiler-builtins/libm/src/math/log2f.rs
View File

@@ -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 */

15
library/compiler-builtins/libm/src/math/logf.rs
View File

@@ -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]). */

46
library/compiler-builtins/libm/src/math/mod.rs
View File

@@ -7,6 +7,7 @@ macro_rules! force_eval {
}
mod ceilf;
mod cosf;
mod expf;
mod fabs;
mod fabsf;
@@ -16,8 +17,11 @@ mod fmod;
mod fmodf;
mod hypot;
mod hypotf;
mod log;
mod log10;
mod log10f;
mod log1p;
mod log1pf;
mod log2;
mod log2f;
mod logf;
@@ -31,11 +35,39 @@ 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;
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;
pub use self::log::log;
pub use self::log10::log10;
pub use self::log10f::log10f;
pub use self::log1p::log1p;
pub use self::log1pf::log1pf;
pub use self::log2::log2;
pub use self::log2f::log2f;
pub use self::logf::logf;
pub use self::powf::powf;
pub use self::round::round;
pub use self::roundf::roundf;
pub use self::scalbn::scalbn;
pub use self::scalbnf::scalbnf;
pub use self::sqrt::sqrt;
pub use self::sqrtf::sqrtf;
pub use self::trunc::trunc;
pub use self::truncf::truncf;
pub use self::{
ceilf::ceilf, expf::expf, fabs::fabs, fabsf::fabsf, floor::floor, floorf::floorf, fmod::fmod,
fmodf::fmodf, hypot::hypot, hypotf::hypotf, log10::log10, log10f::log10f, log2::log2,
log2f::log2f, logf::logf, powf::powf, round::round, roundf::roundf, scalbn::scalbn,
scalbnf::scalbnf, sqrt::sqrt, sqrtf::sqrtf, trunc::trunc, 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};

15
library/compiler-builtins/libm/src/math/powf.rs
View File

@@ -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];

450
library/compiler-builtins/libm/src/math/rem_pio2_large.rs Normal file
View File

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

46
library/compiler-builtins/libm/src/math/rem_pio2f.rs Normal file
View File

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

78
library/compiler-builtins/libm/src/math/sqrt.rs
View File

@@ -1,3 +1,81 @@
/* origin: FreeBSD /usr/src/lib/msun/src/e_sqrt.c */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* sqrt(x)
* Return correctly rounded sqrt.
* ------------------------------------------
* | Use the hardware sqrt if you have one |
* ------------------------------------------
* Method:
* Bit by bit method using integer arithmetic. (Slow, but portable)
* 1. Normalization
* Scale x to y in [1,4) with even powers of 2:
* find an integer k such that 1 <= (y=x*2^(2k)) < 4, then
* sqrt(x) = 2^k * sqrt(y)
* 2. Bit by bit computation
* Let q = sqrt(y) truncated to i bit after binary point (q = 1),
* i 0
* i+1 2
* s = 2*q , and y = 2 * ( y - q ). (1)
* i i i i
*
* To compute q from q , one checks whether
* i+1 i
*
* -(i+1) 2
* (q + 2 ) <= y. (2)
* i
* -(i+1)
* If (2) is false, then q = q ; otherwise q = q + 2 .
* i+1 i i+1 i
*
* With some algebric manipulation, it is not difficult to see
* that (2) is equivalent to
* -(i+1)
* s + 2 <= y (3)
* i i
*
* The advantage of (3) is that s and y can be computed by
* i i
* the following recurrence formula:
* if (3) is false
*
* s = s , y = y ; (4)
* i+1 i i+1 i
*
* otherwise,
* -i -(i+1)
* s = s + 2 , y = y - s - 2 (5)
* i+1 i i+1 i i
*
* One may easily use induction to prove (4) and (5).
* Note. Since the left hand side of (3) contain only i+2 bits,
* it does not necessary to do a full (53-bit) comparison
* in (3).
* 3. Final rounding
* After generating the 53 bits result, we compute one more bit.
* Together with the remainder, we can decide whether the
* result is exact, bigger than 1/2ulp, or less than 1/2ulp
* (it will never equal to 1/2ulp).
* The rounding mode can be detected by checking whether
* huge + tiny is equal to huge, and whether huge - tiny is
* equal to huge for some floating point number "huge" and "tiny".
*
* Special cases:
* sqrt(+-0) = +-0 ... exact
* sqrt(inf) = inf
* sqrt(-ve) = NaN ... with invalid signal
* sqrt(NaN) = NaN ... with invalid signal for signaling NaN
*/
use core::f64;
const TINY: f64 = 1.0e-300;

15
library/compiler-builtins/libm/src/math/sqrtf.rs
View File

@@ -1,3 +1,18 @@
/* origin: FreeBSD /usr/src/lib/msun/src/e_sqrtf.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 TINY: f32 = 1.0e-30;
#[inline]