Merge branch 'master' into master
This commit is contained in:
commit
d6edc4b6ab
34 changed files with 735 additions and 92 deletions
|
|
@ -5,6 +5,19 @@ This project adheres to [Semantic Versioning](http://semver.org/).
|
|||
|
||||
## [Unreleased]
|
||||
|
||||
### Added
|
||||
|
||||
- atan2f
|
||||
- cos
|
||||
- coshf
|
||||
- fmaf
|
||||
- sin
|
||||
- sinh
|
||||
- sinhf
|
||||
- tan
|
||||
- tanh
|
||||
- tanhf
|
||||
|
||||
## [v0.1.1] - 2018-07-14
|
||||
|
||||
### Added
|
||||
|
|
|
|||
|
|
@ -15,11 +15,6 @@ main() {
|
|||
# generate tests
|
||||
cargo run --package test-generator --target x86_64-unknown-linux-musl
|
||||
|
||||
if cargo fmt --version >/dev/null 2>&1; then
|
||||
# nicer syntax error messages (if any)
|
||||
cargo fmt
|
||||
fi
|
||||
|
||||
# run tests
|
||||
cross test --target $TARGET --release
|
||||
|
||||
|
|
|
|||
|
|
@ -14,18 +14,19 @@
|
|||
|
||||
mod math;
|
||||
|
||||
#[cfg(todo)]
|
||||
use core::{f32, f64};
|
||||
|
||||
pub use math::*;
|
||||
|
||||
/// Approximate equality with 1 ULP of tolerance
|
||||
#[doc(hidden)]
|
||||
#[inline]
|
||||
pub fn _eqf(a: u32, b: u32) -> bool {
|
||||
(a as i32).wrapping_sub(b as i32).abs() <= 1
|
||||
}
|
||||
|
||||
#[doc(hidden)]
|
||||
#[inline]
|
||||
pub fn _eq(a: u64, b: u64) -> bool {
|
||||
(a as i64).wrapping_sub(b as i64).abs() <= 1
|
||||
}
|
||||
|
|
@ -33,7 +34,7 @@ pub fn _eq(a: u64, b: u64) -> bool {
|
|||
/// Math support for `f32`
|
||||
///
|
||||
/// This trait is sealed and cannot be implemented outside of `libm`.
|
||||
pub trait F32Ext: private::Sealed {
|
||||
pub trait F32Ext: private::Sealed + Sized {
|
||||
fn floor(self) -> Self;
|
||||
|
||||
fn ceil(self) -> Self;
|
||||
|
|
@ -44,21 +45,17 @@ pub trait F32Ext: private::Sealed {
|
|||
|
||||
fn fdim(self, rhs: Self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
fn fract(self) -> Self;
|
||||
|
||||
fn abs(self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
fn signum(self) -> Self;
|
||||
// NOTE depends on unstable intrinsics::copysignf32
|
||||
// fn signum(self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
fn mul_add(self, a: Self, b: Self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
fn div_euc(self, rhs: Self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
fn mod_euc(self, rhs: Self) -> Self;
|
||||
|
||||
// NOTE depends on unstable intrinsics::powif32
|
||||
|
|
@ -98,9 +95,11 @@ pub trait F32Ext: private::Sealed {
|
|||
|
||||
fn atan2(self, other: Self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
#[inline]
|
||||
fn sin_cos(self) -> (Self, Self) {
|
||||
fn sin_cos(self) -> (Self, Self)
|
||||
where
|
||||
Self: Copy,
|
||||
{
|
||||
(self.sin(), self.cos())
|
||||
}
|
||||
|
||||
|
|
@ -114,13 +113,10 @@ pub trait F32Ext: private::Sealed {
|
|||
|
||||
fn tanh(self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
fn asinh(self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
fn acosh(self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
fn atanh(self) -> Self;
|
||||
}
|
||||
|
||||
|
|
@ -150,7 +146,6 @@ impl F32Ext for f32 {
|
|||
fdimf(self, rhs)
|
||||
}
|
||||
|
||||
#[cfg(todo)]
|
||||
#[inline]
|
||||
fn fract(self) -> Self {
|
||||
self - self.trunc()
|
||||
|
|
@ -161,13 +156,11 @@ impl F32Ext for f32 {
|
|||
fabsf(self)
|
||||
}
|
||||
|
||||
#[cfg(todo)]
|
||||
#[inline]
|
||||
fn mul_add(self, a: Self, b: Self) -> Self {
|
||||
fmaf(self, a, b)
|
||||
}
|
||||
|
||||
#[cfg(todo)]
|
||||
#[inline]
|
||||
fn div_euc(self, rhs: Self) -> Self {
|
||||
let q = (self / rhs).trunc();
|
||||
|
|
@ -177,7 +170,6 @@ impl F32Ext for f32 {
|
|||
q
|
||||
}
|
||||
|
||||
#[cfg(todo)]
|
||||
#[inline]
|
||||
fn mod_euc(self, rhs: f32) -> f32 {
|
||||
let r = self % rhs;
|
||||
|
|
@ -298,7 +290,6 @@ impl F32Ext for f32 {
|
|||
tanhf(self)
|
||||
}
|
||||
|
||||
#[cfg(todo)]
|
||||
#[inline]
|
||||
fn asinh(self) -> Self {
|
||||
if self == f32::NEG_INFINITY {
|
||||
|
|
@ -308,7 +299,6 @@ impl F32Ext for f32 {
|
|||
}
|
||||
}
|
||||
|
||||
#[cfg(todo)]
|
||||
#[inline]
|
||||
fn acosh(self) -> Self {
|
||||
match self {
|
||||
|
|
@ -317,7 +307,6 @@ impl F32Ext for f32 {
|
|||
}
|
||||
}
|
||||
|
||||
#[cfg(todo)]
|
||||
#[inline]
|
||||
fn atanh(self) -> Self {
|
||||
0.5 * ((2.0 * self) / (1.0 - self)).ln_1p()
|
||||
|
|
@ -327,7 +316,7 @@ impl F32Ext for f32 {
|
|||
/// Math support for `f64`
|
||||
///
|
||||
/// This trait is sealed and cannot be implemented outside of `libm`.
|
||||
pub trait F64Ext: private::Sealed {
|
||||
pub trait F64Ext: private::Sealed + Sized {
|
||||
fn floor(self) -> Self;
|
||||
|
||||
fn ceil(self) -> Self;
|
||||
|
|
@ -338,20 +327,17 @@ pub trait F64Ext: private::Sealed {
|
|||
|
||||
fn fdim(self, rhs: Self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
fn fract(self) -> Self;
|
||||
|
||||
fn abs(self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
fn signum(self) -> Self;
|
||||
// NOTE depends on unstable intrinsics::copysignf64
|
||||
// fn signum(self) -> Self;
|
||||
|
||||
fn mul_add(self, a: Self, b: Self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
fn div_euc(self, rhs: Self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
fn mod_euc(self, rhs: Self) -> Self;
|
||||
|
||||
// NOTE depends on unstable intrinsics::powif64
|
||||
|
|
@ -384,19 +370,19 @@ pub trait F64Ext: private::Sealed {
|
|||
|
||||
fn tan(self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
fn asin(self) -> Self;
|
||||
|
||||
fn acos(self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
fn atan(self) -> Self;
|
||||
|
||||
fn atan2(self, other: Self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
#[inline]
|
||||
fn sin_cos(self) -> (Self, Self) {
|
||||
fn sin_cos(self) -> (Self, Self)
|
||||
where
|
||||
Self: Copy,
|
||||
{
|
||||
(self.sin(), self.cos())
|
||||
}
|
||||
|
||||
|
|
@ -406,19 +392,14 @@ pub trait F64Ext: private::Sealed {
|
|||
|
||||
fn sinh(self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
fn cosh(self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
fn tanh(self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
fn asinh(self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
fn acosh(self) -> Self;
|
||||
|
||||
#[cfg(todo)]
|
||||
fn atanh(self) -> Self;
|
||||
}
|
||||
|
||||
|
|
@ -447,7 +428,7 @@ impl F64Ext for f64 {
|
|||
fn fdim(self, rhs: Self) -> Self {
|
||||
fdim(self, rhs)
|
||||
}
|
||||
#[cfg(todo)]
|
||||
|
||||
#[inline]
|
||||
fn fract(self) -> Self {
|
||||
self - self.trunc()
|
||||
|
|
@ -463,7 +444,6 @@ impl F64Ext for f64 {
|
|||
fma(self, a, b)
|
||||
}
|
||||
|
||||
#[cfg(todo)]
|
||||
#[inline]
|
||||
fn div_euc(self, rhs: Self) -> Self {
|
||||
let q = (self / rhs).trunc();
|
||||
|
|
@ -473,9 +453,8 @@ impl F64Ext for f64 {
|
|||
q
|
||||
}
|
||||
|
||||
#[cfg(todo)]
|
||||
#[inline]
|
||||
fn mod_euc(self, rhs: f32) -> f32 {
|
||||
fn mod_euc(self, rhs: f64) -> f64 {
|
||||
let r = self % rhs;
|
||||
if r < 0.0 {
|
||||
r + rhs.abs()
|
||||
|
|
@ -550,7 +529,6 @@ impl F64Ext for f64 {
|
|||
tan(self)
|
||||
}
|
||||
|
||||
#[cfg(todo)]
|
||||
#[inline]
|
||||
fn asin(self) -> Self {
|
||||
asin(self)
|
||||
|
|
@ -561,7 +539,6 @@ impl F64Ext for f64 {
|
|||
acos(self)
|
||||
}
|
||||
|
||||
#[cfg(todo)]
|
||||
#[inline]
|
||||
fn atan(self) -> Self {
|
||||
atan(self)
|
||||
|
|
@ -587,19 +564,16 @@ impl F64Ext for f64 {
|
|||
sinh(self)
|
||||
}
|
||||
|
||||
#[cfg(todo)]
|
||||
#[inline]
|
||||
fn cosh(self) -> Self {
|
||||
cosh(self)
|
||||
}
|
||||
|
||||
#[cfg(todo)]
|
||||
#[inline]
|
||||
fn tanh(self) -> Self {
|
||||
tanh(self)
|
||||
}
|
||||
|
||||
#[cfg(todo)]
|
||||
#[inline]
|
||||
fn asinh(self) -> Self {
|
||||
if self == f64::NEG_INFINITY {
|
||||
|
|
@ -609,7 +583,6 @@ impl F64Ext for f64 {
|
|||
}
|
||||
}
|
||||
|
||||
#[cfg(todo)]
|
||||
#[inline]
|
||||
fn acosh(self) -> Self {
|
||||
match self {
|
||||
|
|
@ -618,7 +591,6 @@ impl F64Ext for f64 {
|
|||
}
|
||||
}
|
||||
|
||||
#[cfg(todo)]
|
||||
#[inline]
|
||||
fn atanh(self) -> Self {
|
||||
0.5 * ((2.0 * self) / (1.0 - self)).ln_1p()
|
||||
|
|
|
|||
|
|
@ -1,3 +1,18 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/e_acosf.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::sqrtf::sqrtf;
|
||||
|
||||
const PIO2_HI: f32 = 1.5707962513e+00; /* 0x3fc90fda */
|
||||
|
|
@ -7,6 +22,7 @@ const P_S1: f32 = -4.2743422091e-02;
|
|||
const P_S2: f32 = -8.6563630030e-03;
|
||||
const Q_S1: f32 = -7.0662963390e-01;
|
||||
|
||||
#[inline]
|
||||
fn r(z: f32) -> f32 {
|
||||
let p = z * (P_S0 + z * (P_S1 + z * P_S2));
|
||||
let q = 1. + z * Q_S1;
|
||||
|
|
|
|||
|
|
@ -55,12 +55,14 @@ const Q_S2: f64 = 2.02094576023350569471e+00; /* 0x40002AE5, 0x9C598AC8 */
|
|||
const Q_S3: f64 = -6.88283971605453293030e-01; /* 0xBFE6066C, 0x1B8D0159 */
|
||||
const Q_S4: f64 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
|
||||
|
||||
#[inline]
|
||||
fn comp_r(z: f64) -> f64 {
|
||||
let p = z * (P_S0 + z * (P_S1 + z * (P_S2 + z * (P_S3 + z * (P_S4 + z * P_S5)))));
|
||||
let q = 1.0 + z * (Q_S1 + z * (Q_S2 + z * (Q_S3 + z * Q_S4)));
|
||||
return p / q;
|
||||
}
|
||||
|
||||
#[inline]
|
||||
pub fn asin(mut x: f64) -> f64 {
|
||||
let z: f64;
|
||||
let r: f64;
|
||||
|
|
|
|||
|
|
@ -1,3 +1,18 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/e_asinf.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::fabsf;
|
||||
use super::sqrt::sqrt;
|
||||
|
||||
|
|
@ -9,6 +24,7 @@ const P_S1: f32 = -4.2743422091e-02;
|
|||
const P_S2: f32 = -8.6563630030e-03;
|
||||
const Q_S1: f32 = -7.0662963390e-01;
|
||||
|
||||
#[inline]
|
||||
fn r(z: f32) -> f32 {
|
||||
let p = z * (P_S0 + z * (P_S1 + z * P_S2));
|
||||
let q = 1. + z * Q_S1;
|
||||
|
|
|
|||
170
library/compiler-builtins/libm/src/math/atan.rs
Normal file
170
library/compiler-builtins/libm/src/math/atan.rs
Normal file
|
|
@ -0,0 +1,170 @@
|
|||
/* atan(x)
|
||||
* Method
|
||||
* 1. Reduce x to positive by atan(x) = -atan(-x).
|
||||
* 2. According to the integer k=4t+0.25 chopped, t=x, the argument
|
||||
* is further reduced to one of the following intervals and the
|
||||
* arctangent of t is evaluated by the corresponding formula:
|
||||
*
|
||||
* [0,7/16] atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
|
||||
* [7/16,11/16] atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) )
|
||||
* [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) )
|
||||
* [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) )
|
||||
* [39/16,INF] atan(x) = atan(INF) + atan( -1/t )
|
||||
*
|
||||
* 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.
|
||||
*/
|
||||
|
||||
use super::fabs;
|
||||
use core::f64;
|
||||
|
||||
const ATANHI: [f64; 4] = [
|
||||
4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */
|
||||
7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */
|
||||
9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */
|
||||
1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */
|
||||
];
|
||||
|
||||
const ATANLO: [f64; 4] = [
|
||||
2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */
|
||||
3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */
|
||||
1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */
|
||||
6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */
|
||||
];
|
||||
|
||||
const AT: [f64; 11] = [
|
||||
3.33333333333329318027e-01, /* 0x3FD55555, 0x5555550D */
|
||||
-1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */
|
||||
1.42857142725034663711e-01, /* 0x3FC24924, 0x920083FF */
|
||||
-1.11111104054623557880e-01, /* 0xBFBC71C6, 0xFE231671 */
|
||||
9.09088713343650656196e-02, /* 0x3FB745CD, 0xC54C206E */
|
||||
-7.69187620504482999495e-02, /* 0xBFB3B0F2, 0xAF749A6D */
|
||||
6.66107313738753120669e-02, /* 0x3FB10D66, 0xA0D03D51 */
|
||||
-5.83357013379057348645e-02, /* 0xBFADDE2D, 0x52DEFD9A */
|
||||
4.97687799461593236017e-02, /* 0x3FA97B4B, 0x24760DEB */
|
||||
-3.65315727442169155270e-02, /* 0xBFA2B444, 0x2C6A6C2F */
|
||||
1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */
|
||||
];
|
||||
|
||||
#[inline]
|
||||
pub fn atan(x: f64) -> f64 {
|
||||
let mut x = x;
|
||||
let mut ix = (x.to_bits() >> 32) as u32;
|
||||
let sign = ix >> 31;
|
||||
ix &= 0x7fff_ffff;
|
||||
if ix >= 0x4410_0000 {
|
||||
if x.is_nan() {
|
||||
return x;
|
||||
}
|
||||
|
||||
let z = ATANHI[3] + f64::from_bits(0x0380_0000); // 0x1p-120f
|
||||
return if sign != 0 { -z } else { z };
|
||||
}
|
||||
|
||||
let id = if ix < 0x3fdc_0000 {
|
||||
/* |x| < 0.4375 */
|
||||
if ix < 0x3e40_0000 {
|
||||
/* |x| < 2^-27 */
|
||||
if ix < 0x0010_0000 {
|
||||
/* raise underflow for subnormal x */
|
||||
force_eval!(x as f32);
|
||||
}
|
||||
|
||||
return x;
|
||||
}
|
||||
|
||||
-1
|
||||
} else {
|
||||
x = fabs(x);
|
||||
if ix < 0x3ff30000 {
|
||||
/* |x| < 1.1875 */
|
||||
if ix < 0x3fe60000 {
|
||||
/* 7/16 <= |x| < 11/16 */
|
||||
x = (2. * x - 1.) / (2. + x);
|
||||
0
|
||||
} else {
|
||||
/* 11/16 <= |x| < 19/16 */
|
||||
x = (x - 1.) / (x + 1.);
|
||||
1
|
||||
}
|
||||
} else {
|
||||
if ix < 0x40038000 {
|
||||
/* |x| < 2.4375 */
|
||||
x = (x - 1.5) / (1. + 1.5 * x);
|
||||
2
|
||||
} else {
|
||||
/* 2.4375 <= |x| < 2^66 */
|
||||
x = -1. / x;
|
||||
3
|
||||
}
|
||||
}
|
||||
};
|
||||
|
||||
let z = x * x;
|
||||
let w = z * z;
|
||||
/* break sum from i=0 to 10 AT[i]z**(i+1) into odd and even poly */
|
||||
let s1 = z * (AT[0] + w * (AT[2] + w * (AT[4] + w * (AT[6] + w * (AT[8] + w * AT[10])))));
|
||||
let s2 = w * (AT[1] + w * (AT[3] + w * (AT[5] + w * (AT[7] + w * AT[9]))));
|
||||
|
||||
if id < 0 {
|
||||
return x - x * (s1 + s2);
|
||||
}
|
||||
|
||||
let z = ATANHI[id as usize] - (x * (s1 + s2) - ATANLO[id as usize] - x);
|
||||
|
||||
if sign != 0 {
|
||||
-z
|
||||
} else {
|
||||
z
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use super::atan;
|
||||
use core::f64;
|
||||
|
||||
#[test]
|
||||
fn sanity_check() {
|
||||
for (input, answer) in [
|
||||
(3.0_f64.sqrt() / 3.0, f64::consts::FRAC_PI_6),
|
||||
(1.0, f64::consts::FRAC_PI_4),
|
||||
(3.0_f64.sqrt(), f64::consts::FRAC_PI_3),
|
||||
(-3.0_f64.sqrt() / 3.0, -f64::consts::FRAC_PI_6),
|
||||
(-1.0, -f64::consts::FRAC_PI_4),
|
||||
(-3.0_f64.sqrt(), -f64::consts::FRAC_PI_3),
|
||||
].iter()
|
||||
{
|
||||
assert!(
|
||||
(atan(*input) - answer) / answer < 1e-5,
|
||||
"\natan({:.4}/16) = {:.4}, actual: {}",
|
||||
input * 16.0,
|
||||
answer,
|
||||
atan(*input)
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn zero() {
|
||||
assert_eq!(atan(0.0), 0.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn infinity() {
|
||||
assert_eq!(atan(f64::INFINITY), f64::consts::FRAC_PI_2);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn minus_infinity() {
|
||||
assert_eq!(atan(f64::NEG_INFINITY), -f64::consts::FRAC_PI_2);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn nan() {
|
||||
assert!(atan(f64::NAN).is_nan());
|
||||
}
|
||||
}
|
||||
|
|
@ -1,3 +1,18 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/s_atanf.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;
|
||||
|
||||
const ATAN_HI: [f32; 4] = [
|
||||
|
|
|
|||
|
|
@ -41,6 +41,7 @@ use super::{k_cos, k_sin, rem_pio2};
|
|||
// Accuracy:
|
||||
// TRIG(x) returns trig(x) nearly rounded
|
||||
//
|
||||
#[inline]
|
||||
pub fn cos(x: f64) -> f64 {
|
||||
let ix = (f64::to_bits(x) >> 32) as u32 & 0x7fffffff;
|
||||
|
||||
|
|
|
|||
|
|
@ -1,3 +1,19 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/s_cosf.c */
|
||||
/*
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
||||
* Optimized by Bruce D. Evans.
|
||||
*/
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
use super::{k_cosf, k_sinf, rem_pio2f};
|
||||
|
||||
use core::f64::consts::FRAC_PI_2;
|
||||
|
|
|
|||
33
library/compiler-builtins/libm/src/math/cosh.rs
Normal file
33
library/compiler-builtins/libm/src/math/cosh.rs
Normal file
|
|
@ -0,0 +1,33 @@
|
|||
use super::exp;
|
||||
use super::expm1;
|
||||
use super::k_expo2;
|
||||
|
||||
#[inline]
|
||||
pub fn cosh(mut x: f64) -> f64 {
|
||||
/* |x| */
|
||||
let mut ix = x.to_bits();
|
||||
ix &= 0x7fffffffffffffff;
|
||||
x = f64::from_bits(ix);
|
||||
let w = ix >> 32;
|
||||
|
||||
/* |x| < log(2) */
|
||||
if w < 0x3fe62e42 {
|
||||
if w < 0x3ff00000 - (26 << 20) {
|
||||
let x1p120 = f64::from_bits(0x4770000000000000);
|
||||
force_eval!(x + x1p120);
|
||||
return 1.;
|
||||
}
|
||||
let t = expm1(x); // exponential minus 1
|
||||
return 1. + t * t / (2. * (1. + t));
|
||||
}
|
||||
|
||||
/* |x| < log(DBL_MAX) */
|
||||
if w < 0x40862e42 {
|
||||
let t = exp(x);
|
||||
/* note: if x>log(0x1p26) then the 1/t is not needed */
|
||||
return 0.5 * (t + 1. / t);
|
||||
}
|
||||
|
||||
/* |x| > log(DBL_MAX) or nan */
|
||||
k_expo2(x)
|
||||
}
|
||||
|
|
@ -1,3 +1,15 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/s_expm1.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::f64;
|
||||
|
||||
const O_THRESHOLD: f64 = 7.09782712893383973096e+02; /* 0x40862E42, 0xFEFA39EF */
|
||||
|
|
@ -11,7 +23,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)]
|
||||
#[inline]
|
||||
pub fn expm1(mut x: f64) -> f64 {
|
||||
let hi: f64;
|
||||
let lo: f64;
|
||||
|
|
|
|||
|
|
@ -1,3 +1,18 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/s_expm1f.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 O_THRESHOLD: f32 = 8.8721679688e+01; /* 0x42b17180 */
|
||||
const LN2_HI: f32 = 6.9313812256e-01; /* 0x3f317180 */
|
||||
const LN2_LO: f32 = 9.0580006145e-06; /* 0x3717f7d1 */
|
||||
|
|
|
|||
|
|
@ -1,7 +1,8 @@
|
|||
use super::{combine_words, exp};
|
||||
|
||||
/* exp(x)/2 for x >= log(DBL_MAX), slightly better than 0.5*exp(x/2)*exp(x/2) */
|
||||
pub(crate) fn expo2(x: f64) -> f64 {
|
||||
#[inline]
|
||||
pub fn expo2(x: f64) -> f64 {
|
||||
/* k is such that k*ln2 has minimal relative error and x - kln2 > log(DBL_MIN) */
|
||||
const K: i32 = 2043;
|
||||
let kln2 = f64::from_bits(0x40962066151add8b);
|
||||
|
|
|
|||
33
library/compiler-builtins/libm/src/math/fenv.rs
Normal file
33
library/compiler-builtins/libm/src/math/fenv.rs
Normal file
|
|
@ -0,0 +1,33 @@
|
|||
// src: musl/src/fenv/fenv.c
|
||||
/* Dummy functions for archs lacking fenv implementation */
|
||||
|
||||
pub const FE_UNDERFLOW: i32 = 0;
|
||||
pub const FE_INEXACT: i32 = 0;
|
||||
|
||||
pub const FE_TONEAREST: i32 = 0;
|
||||
pub const FE_TOWARDZERO: i32 = 0;
|
||||
|
||||
#[inline]
|
||||
pub fn feclearexcept(_mask: i32) -> i32 {
|
||||
0
|
||||
}
|
||||
|
||||
#[inline]
|
||||
pub fn feraiseexcept(_mask: i32) -> i32 {
|
||||
0
|
||||
}
|
||||
|
||||
#[inline]
|
||||
pub fn fetestexcept(_mask: i32) -> i32 {
|
||||
0
|
||||
}
|
||||
|
||||
#[inline]
|
||||
pub fn fegetround() -> i32 {
|
||||
FE_TONEAREST
|
||||
}
|
||||
|
||||
#[inline]
|
||||
pub fn fesetround(_r: i32) -> i32 {
|
||||
0
|
||||
}
|
||||
100
library/compiler-builtins/libm/src/math/fmaf.rs
Normal file
100
library/compiler-builtins/libm/src/math/fmaf.rs
Normal file
|
|
@ -0,0 +1,100 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/s_fmaf.c */
|
||||
/*-
|
||||
* Copyright (c) 2005-2011 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 core::f32;
|
||||
use core::ptr::read_volatile;
|
||||
|
||||
use super::fenv::{
|
||||
feclearexcept, fegetround, feraiseexcept, fesetround, fetestexcept, FE_INEXACT, FE_TONEAREST,
|
||||
FE_TOWARDZERO, FE_UNDERFLOW,
|
||||
};
|
||||
|
||||
/*
|
||||
* Fused multiply-add: Compute x * y + z with a single rounding error.
|
||||
*
|
||||
* A double has more than twice as much precision than a float, so
|
||||
* direct double-precision arithmetic suffices, except where double
|
||||
* rounding occurs.
|
||||
*/
|
||||
#[inline]
|
||||
pub fn fmaf(x: f32, y: f32, mut z: f32) -> f32 {
|
||||
let xy: f64;
|
||||
let mut result: f64;
|
||||
let mut ui: u64;
|
||||
let e: i32;
|
||||
|
||||
xy = x as f64 * y as f64;
|
||||
result = xy + z as f64;
|
||||
ui = result.to_bits();
|
||||
e = (ui >> 52) as i32 & 0x7ff;
|
||||
/* Common case: The double precision result is fine. */
|
||||
if (
|
||||
/* not a halfway case */
|
||||
ui & 0x1fffffff) != 0x10000000 ||
|
||||
/* NaN */
|
||||
e == 0x7ff ||
|
||||
/* exact */
|
||||
(result - xy == z as f64 && result - z as f64 == xy) ||
|
||||
/* not round-to-nearest */
|
||||
fegetround() != FE_TONEAREST
|
||||
{
|
||||
/*
|
||||
underflow may not be raised correctly, example:
|
||||
fmaf(0x1p-120f, 0x1p-120f, 0x1p-149f)
|
||||
*/
|
||||
if e < 0x3ff - 126 && e >= 0x3ff - 149 && fetestexcept(FE_INEXACT) != 0 {
|
||||
feclearexcept(FE_INEXACT);
|
||||
// prevent `xy + vz` from being CSE'd with `xy + z` above
|
||||
let vz: f32 = unsafe { read_volatile(&z) };
|
||||
result = xy + vz as f64;
|
||||
if fetestexcept(FE_INEXACT) != 0 {
|
||||
feraiseexcept(FE_UNDERFLOW);
|
||||
} else {
|
||||
feraiseexcept(FE_INEXACT);
|
||||
}
|
||||
}
|
||||
z = result as f32;
|
||||
return z;
|
||||
}
|
||||
|
||||
/*
|
||||
* If result is inexact, and exactly halfway between two float values,
|
||||
* we need to adjust the low-order bit in the direction of the error.
|
||||
*/
|
||||
fesetround(FE_TOWARDZERO);
|
||||
// prevent `vxy + z` from being CSE'd with `xy + z` above
|
||||
let vxy: f64 = unsafe { read_volatile(&xy) };
|
||||
let mut adjusted_result: f64 = vxy + z as f64;
|
||||
fesetround(FE_TONEAREST);
|
||||
if result == adjusted_result {
|
||||
ui = adjusted_result.to_bits();
|
||||
ui += 1;
|
||||
adjusted_result = f64::from_bits(ui);
|
||||
}
|
||||
z = adjusted_result as f32;
|
||||
z
|
||||
}
|
||||
|
|
@ -1,3 +1,19 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/k_cosf.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.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
/* |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 */
|
||||
|
|
@ -5,7 +21,7 @@ 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 {
|
||||
pub fn k_cosf(x: f64) -> f32 {
|
||||
let z = x * x;
|
||||
let w = z * z;
|
||||
let r = C2 + z * C3;
|
||||
|
|
|
|||
14
library/compiler-builtins/libm/src/math/k_expo2.rs
Normal file
14
library/compiler-builtins/libm/src/math/k_expo2.rs
Normal file
|
|
@ -0,0 +1,14 @@
|
|||
use super::exp;
|
||||
|
||||
/* k is such that k*ln2 has minimal relative error and x - kln2 > log(FLT_MIN) */
|
||||
const K: i32 = 2043;
|
||||
|
||||
/* expf(x)/2 for x >= log(FLT_MAX), slightly better than 0.5f*expf(x/2)*expf(x/2) */
|
||||
#[inline]
|
||||
pub(crate) fn k_expo2(x: f64) -> f64 {
|
||||
let k_ln2 = f64::from_bits(0x40962066151add8b);
|
||||
/* note that k is odd and scale*scale overflows */
|
||||
let scale = f64::from_bits(((((0x3ff + K / 2) as u32) << 20) as u64) << 32);
|
||||
/* exp(x - k ln2) * 2**(k-1) */
|
||||
exp(x - k_ln2) * scale * scale
|
||||
}
|
||||
|
|
@ -5,7 +5,7 @@ const K: i32 = 235;
|
|||
|
||||
/* expf(x)/2 for x >= log(FLT_MAX), slightly better than 0.5f*expf(x/2)*expf(x/2) */
|
||||
#[inline]
|
||||
pub(crate) fn k_expo2f(x: f32) -> f32 {
|
||||
pub fn k_expo2f(x: f32) -> f32 {
|
||||
let k_ln2 = f32::from_bits(0x4322e3bc);
|
||||
/* note that k is odd and scale*scale overflows */
|
||||
let scale = f32::from_bits(((0x7f + K / 2) as u32) << 23);
|
||||
|
|
|
|||
|
|
@ -1,3 +1,19 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/k_sinf.c */
|
||||
/*
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
||||
* Optimized by Bruce D. Evans.
|
||||
*/
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
/* |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 */
|
||||
|
|
@ -5,7 +21,7 @@ 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 {
|
||||
pub fn k_sinf(x: f64) -> f32 {
|
||||
let z = x * x;
|
||||
let w = z * z;
|
||||
let r = S3 + z * S4;
|
||||
|
|
|
|||
|
|
@ -58,7 +58,8 @@ static T: [f64; 13] = [
|
|||
const PIO4: f64 = 7.85398163397448278999e-01; /* 3FE921FB, 54442D18 */
|
||||
const PIO4_LO: f64 = 3.06161699786838301793e-17; /* 3C81A626, 33145C07 */
|
||||
|
||||
pub(crate) fn k_tan(mut x: f64, mut y: f64, odd: i32) -> f64 {
|
||||
#[inline]
|
||||
pub fn k_tan(mut x: f64, mut y: f64, odd: i32) -> f64 {
|
||||
let hx = (f64::to_bits(x) >> 32) as u32;
|
||||
let big = (hx & 0x7fffffff) >= 0x3FE59428; /* |x| >= 0.6744 */
|
||||
if big {
|
||||
|
|
|
|||
|
|
@ -1,3 +1,14 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/k_tan.c */
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright 2004 Sun Microsystems, Inc. All Rights Reserved.
|
||||
*
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
/* |tan(x)/x - t(x)| < 2**-25.5 (~[-2e-08, 2e-08]). */
|
||||
const T: [f64; 6] = [
|
||||
0.333331395030791399758, /* 0x15554d3418c99f.0p-54 */
|
||||
|
|
@ -9,7 +20,7 @@ const T: [f64; 6] = [
|
|||
];
|
||||
|
||||
#[inline]
|
||||
pub(crate) fn k_tanf(x: f64, odd: bool) -> f32 {
|
||||
pub fn k_tanf(x: f64, odd: bool) -> f32 {
|
||||
let z = x * x;
|
||||
/*
|
||||
* Split up the polynomial into small independent terms to give
|
||||
|
|
|
|||
|
|
@ -10,7 +10,6 @@ macro_rules! force_eval {
|
|||
mod acos;
|
||||
mod acosf;
|
||||
mod asin;
|
||||
mod asinf;
|
||||
mod atan2;
|
||||
mod atan2f;
|
||||
mod atanf;
|
||||
|
|
@ -20,6 +19,7 @@ mod ceil;
|
|||
mod ceilf;
|
||||
mod cos;
|
||||
mod cosf;
|
||||
mod cosh;
|
||||
mod coshf;
|
||||
mod exp;
|
||||
mod exp2;
|
||||
|
|
@ -34,6 +34,7 @@ mod fdimf;
|
|||
mod floor;
|
||||
mod floorf;
|
||||
mod fma;
|
||||
mod fmaf;
|
||||
mod fmod;
|
||||
mod fmodf;
|
||||
mod hypot;
|
||||
|
|
@ -59,6 +60,7 @@ mod sqrt;
|
|||
mod sqrtf;
|
||||
mod tan;
|
||||
mod tanf;
|
||||
mod tanh;
|
||||
mod tanhf;
|
||||
mod trunc;
|
||||
mod truncf;
|
||||
|
|
@ -69,6 +71,7 @@ pub use self::acosf::acosf;
|
|||
pub use self::asin::asin;
|
||||
pub use self::asinf::asinf;
|
||||
pub use self::atan2::atan2;
|
||||
pub use self::atan::atan;
|
||||
pub use self::atan2f::atan2f;
|
||||
pub use self::atanf::atanf;
|
||||
pub use self::cbrt::cbrt;
|
||||
|
|
@ -77,6 +80,7 @@ pub use self::ceil::ceil;
|
|||
pub use self::ceilf::ceilf;
|
||||
pub use self::cos::cos;
|
||||
pub use self::cosf::cosf;
|
||||
pub use self::cosh::cosh;
|
||||
pub use self::coshf::coshf;
|
||||
pub use self::exp::exp;
|
||||
pub use self::exp2::exp2;
|
||||
|
|
@ -91,6 +95,7 @@ pub use self::fdimf::fdimf;
|
|||
pub use self::floor::floor;
|
||||
pub use self::floorf::floorf;
|
||||
pub use self::fma::fma;
|
||||
pub use self::fmaf::fmaf;
|
||||
pub use self::fmod::fmod;
|
||||
pub use self::fmodf::fmodf;
|
||||
pub use self::hypot::hypot;
|
||||
|
|
@ -116,14 +121,17 @@ pub use self::sqrt::sqrt;
|
|||
pub use self::sqrtf::sqrtf;
|
||||
pub use self::tan::tan;
|
||||
pub use self::tanf::tanf;
|
||||
pub use self::tanh::tanh;
|
||||
pub use self::tanhf::tanhf;
|
||||
pub use self::trunc::trunc;
|
||||
pub use self::truncf::truncf;
|
||||
|
||||
// Private modules
|
||||
mod expo2;
|
||||
mod fenv;
|
||||
mod k_cos;
|
||||
mod k_cosf;
|
||||
mod k_expo2;
|
||||
mod k_expo2f;
|
||||
mod k_sin;
|
||||
mod k_sinf;
|
||||
|
|
@ -137,6 +145,7 @@ mod rem_pio2f;
|
|||
use self::expo2::expo2;
|
||||
use self::k_cos::k_cos;
|
||||
use self::k_cosf::k_cosf;
|
||||
use self::k_expo2::k_expo2;
|
||||
use self::k_expo2f::k_expo2f;
|
||||
use self::k_sin::k_sin;
|
||||
use self::k_sinf::k_sinf;
|
||||
|
|
@ -147,17 +156,18 @@ use self::rem_pio2_large::rem_pio2_large;
|
|||
use self::rem_pio2f::rem_pio2f;
|
||||
|
||||
#[inline]
|
||||
pub fn get_high_word(x: f64) -> u32 {
|
||||
fn get_high_word(x: f64) -> u32 {
|
||||
(x.to_bits() >> 32) as u32
|
||||
}
|
||||
|
||||
#[inline]
|
||||
pub fn get_low_word(x: f64) -> u32 {
|
||||
fn get_low_word(x: f64) -> u32 {
|
||||
x.to_bits() as u32
|
||||
}
|
||||
|
||||
#[allow(dead_code)]
|
||||
#[inline]
|
||||
pub fn with_set_high_word(f: f64, hi: u32) -> f64 {
|
||||
fn with_set_high_word(f: f64, hi: u32) -> f64 {
|
||||
let mut tmp = f.to_bits();
|
||||
tmp &= 0x00000000_ffffffff;
|
||||
tmp |= (hi as u64) << 32;
|
||||
|
|
@ -165,7 +175,7 @@ pub fn with_set_high_word(f: f64, hi: u32) -> f64 {
|
|||
}
|
||||
|
||||
#[inline]
|
||||
pub fn with_set_low_word(f: f64, lo: u32) -> f64 {
|
||||
fn with_set_low_word(f: f64, lo: u32) -> f64 {
|
||||
let mut tmp = f.to_bits();
|
||||
tmp &= 0xffffffff_00000000;
|
||||
tmp |= lo as u64;
|
||||
|
|
|
|||
|
|
@ -42,12 +42,14 @@ const PIO2_3T: f64 = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
|
|||
// use rem_pio2_large() for large x
|
||||
//
|
||||
// caller must handle the case when reduction is not needed: |x| ~<= pi/4 */
|
||||
#[inline]
|
||||
pub fn rem_pio2(x: f64) -> (i32, f64, f64) {
|
||||
let x1p24 = f64::from_bits(0x4170000000000000);
|
||||
|
||||
let sign = (f64::to_bits(x) >> 63) as i32;
|
||||
let ix = (f64::to_bits(x) >> 32) as u32 & 0x7fffffff;
|
||||
|
||||
#[inline]
|
||||
fn medium(x: f64, ix: u32) -> (i32, f64, f64) {
|
||||
/* rint(x/(pi/2)), Assume round-to-nearest. */
|
||||
let f_n = x as f64 * INV_PIO2 + TO_INT - TO_INT;
|
||||
|
|
|
|||
|
|
@ -1,3 +1,15 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/k_rem_pio2.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.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
use super::floor;
|
||||
use super::scalbn;
|
||||
|
||||
|
|
@ -210,7 +222,7 @@ const PIO2: [f64; 8] = [
|
|||
/// 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 {
|
||||
pub 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)
|
||||
|
||||
|
|
|
|||
|
|
@ -1,3 +1,19 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/e_rem_pio2f.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.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
use super::rem_pio2_large;
|
||||
|
||||
use core::f64;
|
||||
|
|
@ -16,7 +32,7 @@ const PIO2_1T: f64 = 1.58932547735281966916e-08; /* 0x3E5110b4, 0x611A6263 */
|
|||
/// 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) {
|
||||
pub fn rem_pio2f(x: f32) -> (i32, f64) {
|
||||
let x64 = x as f64;
|
||||
|
||||
let mut tx: [f64; 1] = [0.];
|
||||
|
|
|
|||
|
|
@ -40,6 +40,7 @@ use super::{k_cos, k_sin, rem_pio2};
|
|||
//
|
||||
// Accuracy:
|
||||
// TRIG(x) returns trig(x) nearly rounded
|
||||
#[inline]
|
||||
pub fn sin(x: f64) -> f64 {
|
||||
let x1p120 = f64::from_bits(0x4770000000000000); // 0x1p120f === 2 ^ 120
|
||||
|
||||
|
|
|
|||
|
|
@ -1,3 +1,19 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/s_sinf.c */
|
||||
/*
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
||||
* Optimized by Bruce D. Evans.
|
||||
*/
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
use super::{k_cosf, k_sinf, rem_pio2f};
|
||||
|
||||
use core::f64::consts::FRAC_PI_2;
|
||||
|
|
|
|||
|
|
@ -4,6 +4,7 @@ use super::{expm1, expo2};
|
|||
// = (exp(x)-1 + (exp(x)-1)/exp(x))/2
|
||||
// = x + x^3/6 + o(x^5)
|
||||
//
|
||||
#[inline]
|
||||
pub fn sinh(x: f64) -> f64 {
|
||||
// union {double f; uint64_t i;} u = {.f = x};
|
||||
// uint32_t w;
|
||||
|
|
|
|||
|
|
@ -39,6 +39,7 @@ use super::{k_tan, rem_pio2};
|
|||
//
|
||||
// Accuracy:
|
||||
// TRIG(x) returns trig(x) nearly rounded
|
||||
#[inline]
|
||||
pub fn tan(x: f64) -> f64 {
|
||||
let x1p120 = f32::from_bits(0x7b800000); // 0x1p120f === 2 ^ 120
|
||||
|
||||
|
|
|
|||
|
|
@ -1,3 +1,19 @@
|
|||
/* origin: FreeBSD /usr/src/lib/msun/src/s_tanf.c */
|
||||
/*
|
||||
* Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
|
||||
* Optimized by Bruce D. Evans.
|
||||
*/
|
||||
/*
|
||||
* ====================================================
|
||||
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
||||
*
|
||||
* Developed at SunPro, a Sun Microsystems, Inc. business.
|
||||
* Permission to use, copy, modify, and distribute this
|
||||
* software is freely granted, provided that this notice
|
||||
* is preserved.
|
||||
* ====================================================
|
||||
*/
|
||||
|
||||
use super::{k_tanf, rem_pio2f};
|
||||
|
||||
use core::f64::consts::FRAC_PI_2;
|
||||
|
|
|
|||
53
library/compiler-builtins/libm/src/math/tanh.rs
Normal file
53
library/compiler-builtins/libm/src/math/tanh.rs
Normal file
|
|
@ -0,0 +1,53 @@
|
|||
use super::expm1;
|
||||
|
||||
/* tanh(x) = (exp(x) - exp(-x))/(exp(x) + exp(-x))
|
||||
* = (exp(2*x) - 1)/(exp(2*x) - 1 + 2)
|
||||
* = (1 - exp(-2*x))/(exp(-2*x) - 1 + 2)
|
||||
*/
|
||||
#[inline]
|
||||
pub fn tanh(mut x: f64) -> f64 {
|
||||
let mut uf: f64 = x;
|
||||
let mut ui: u64 = f64::to_bits(uf);
|
||||
|
||||
let w: u32;
|
||||
let sign: bool;
|
||||
let mut t: f64;
|
||||
|
||||
/* x = |x| */
|
||||
sign = ui >> 63 != 0;
|
||||
ui &= !1 / 2;
|
||||
uf = f64::from_bits(ui);
|
||||
x = uf;
|
||||
w = (ui >> 32) as u32;
|
||||
|
||||
if w > 0x3fe193ea {
|
||||
/* |x| > log(3)/2 ~= 0.5493 or nan */
|
||||
if w > 0x40340000 {
|
||||
/* |x| > 20 or nan */
|
||||
/* note: this branch avoids raising overflow */
|
||||
t = 1.0 - 0.0 / x;
|
||||
} else {
|
||||
t = expm1(2.0 * x);
|
||||
t = 1.0 - 2.0 / (t + 2.0);
|
||||
}
|
||||
} else if w > 0x3fd058ae {
|
||||
/* |x| > log(5/3)/2 ~= 0.2554 */
|
||||
t = expm1(2.0 * x);
|
||||
t = t / (t + 2.0);
|
||||
} else if w >= 0x00100000 {
|
||||
/* |x| >= 0x1p-1022, up to 2ulp error in [0.1,0.2554] */
|
||||
t = expm1(-2.0 * x);
|
||||
t = -t / (t + 2.0);
|
||||
} else {
|
||||
/* |x| is subnormal */
|
||||
/* note: the branch above would not raise underflow in [0x1p-1023,0x1p-1022) */
|
||||
force_eval!(x as f32);
|
||||
t = x;
|
||||
}
|
||||
|
||||
if sign {
|
||||
-t
|
||||
} else {
|
||||
t
|
||||
}
|
||||
}
|
||||
|
|
@ -6,3 +6,4 @@ publish = false
|
|||
|
||||
[dependencies]
|
||||
rand = "0.5.3"
|
||||
itertools = "0.7.8"
|
||||
|
|
|
|||
|
|
@ -4,13 +4,15 @@
|
|||
// NOTE usually the only thing you need to do to test a new math function is to add it to one of the
|
||||
// macro invocations found in the bottom of this file.
|
||||
|
||||
#[macro_use]
|
||||
extern crate itertools;
|
||||
extern crate rand;
|
||||
|
||||
use std::error::Error;
|
||||
use std::fmt::Write as _0;
|
||||
use std::fs::{self, File};
|
||||
use std::io::Write as _1;
|
||||
use std::{i16, u16, u32, u64, u8};
|
||||
use std::{f32, f64, i16, u16, u32, u64, u8};
|
||||
|
||||
use rand::{Rng, SeedableRng, XorShiftRng};
|
||||
|
||||
|
|
@ -34,6 +36,30 @@ fn f64(rng: &mut XorShiftRng) -> f64 {
|
|||
f64::from_bits(sign + exponent + mantissa)
|
||||
}
|
||||
|
||||
const EDGE_CASES32: &[f32] = &[
|
||||
-0.,
|
||||
0.,
|
||||
f32::EPSILON,
|
||||
f32::INFINITY,
|
||||
f32::MAX,
|
||||
f32::MIN,
|
||||
f32::MIN_POSITIVE,
|
||||
f32::NAN,
|
||||
f32::NEG_INFINITY,
|
||||
];
|
||||
|
||||
const EDGE_CASES64: &[f64] = &[
|
||||
-0.,
|
||||
0.,
|
||||
f64::EPSILON,
|
||||
f64::INFINITY,
|
||||
f64::MAX,
|
||||
f64::MIN,
|
||||
f64::MIN_POSITIVE,
|
||||
f64::NAN,
|
||||
f64::NEG_INFINITY,
|
||||
];
|
||||
|
||||
// fn(f32) -> f32
|
||||
macro_rules! f32_f32 {
|
||||
($($intr:ident,)*) => {
|
||||
|
|
@ -45,8 +71,9 @@ macro_rules! f32_f32 {
|
|||
|
||||
$(
|
||||
let mut cases = String::new();
|
||||
for _ in 0..NTESTS {
|
||||
let inp = f32(rng);
|
||||
|
||||
// random inputs
|
||||
for inp in EDGE_CASES32.iter().cloned().chain((0..NTESTS).map(|_| f32(rng))) {
|
||||
let out = unsafe { $intr(inp) };
|
||||
|
||||
let inp = inp.to_bits();
|
||||
|
|
@ -112,11 +139,17 @@ macro_rules! f32f32_f32 {
|
|||
$(fn $intr(_: f32, _: f32) -> f32;)*
|
||||
}
|
||||
|
||||
let mut rng2 = rng.clone();
|
||||
let mut rng3 = rng.clone();
|
||||
$(
|
||||
let mut cases = String::new();
|
||||
for _ in 0..NTESTS {
|
||||
let i1 = f32(rng);
|
||||
let i2 = f32(rng);
|
||||
for (i1, i2) in iproduct!(
|
||||
EDGE_CASES32.iter().cloned(),
|
||||
EDGE_CASES32.iter().cloned()
|
||||
).chain(EDGE_CASES32.iter().map(|i1| (*i1, f32(rng))))
|
||||
.chain(EDGE_CASES32.iter().map(|i2| (f32(&mut rng2), *i2)))
|
||||
.chain((0..NTESTS).map(|_| (f32(&mut rng3), f32(&mut rng3))))
|
||||
{
|
||||
let out = unsafe { $intr(i1, i2) };
|
||||
|
||||
let i1 = i1.to_bits();
|
||||
|
|
@ -186,12 +219,16 @@ macro_rules! f32f32f32_f32 {
|
|||
$(fn $intr(_: f32, _: f32, _: f32) -> f32;)*
|
||||
}
|
||||
|
||||
let mut rng2 = rng.clone();
|
||||
$(
|
||||
let mut cases = String::new();
|
||||
for _ in 0..NTESTS {
|
||||
let i1 = f32(rng);
|
||||
let i2 = f32(rng);
|
||||
let i3 = f32(rng);
|
||||
for (i1, i2, i3) in iproduct!(
|
||||
EDGE_CASES32.iter().cloned(),
|
||||
EDGE_CASES32.iter().cloned(),
|
||||
EDGE_CASES32.iter().cloned()
|
||||
).chain(EDGE_CASES32.iter().map(|i1| (*i1, f32(rng), f32(rng))))
|
||||
.chain((0..NTESTS).map(|_| (f32(&mut rng2), f32(&mut rng2), f32(&mut rng2))))
|
||||
{
|
||||
let out = unsafe { $intr(i1, i2, i3) };
|
||||
|
||||
let i1 = i1.to_bits();
|
||||
|
|
@ -266,10 +303,10 @@ macro_rules! f32i32_f32 {
|
|||
$(fn $intr(_: f32, _: i32) -> f32;)*
|
||||
}
|
||||
|
||||
let mut rng2 = rng.clone();
|
||||
$(
|
||||
let mut cases = String::new();
|
||||
for _ in 0..NTESTS {
|
||||
let i1 = f32(rng);
|
||||
for i1 in EDGE_CASES32.iter().cloned().chain((0..NTESTS).map(|_| f32(&mut rng2))) {
|
||||
let i2 = rng.gen_range(i16::MIN, i16::MAX);
|
||||
let out = unsafe { $intr(i1, i2 as i32) };
|
||||
|
||||
|
|
@ -342,8 +379,7 @@ macro_rules! f64_f64 {
|
|||
|
||||
$(
|
||||
let mut cases = String::new();
|
||||
for _ in 0..NTESTS {
|
||||
let inp = f64(rng);
|
||||
for inp in EDGE_CASES64.iter().cloned().chain((0..NTESTS).map(|_| f64(rng))) {
|
||||
let out = unsafe { $intr(inp) };
|
||||
|
||||
let inp = inp.to_bits();
|
||||
|
|
@ -412,11 +448,17 @@ macro_rules! f64f64_f64 {
|
|||
$(fn $intr(_: f64, _: f64) -> f64;)*
|
||||
}
|
||||
|
||||
let mut rng2 = rng.clone();
|
||||
let mut rng3 = rng.clone();
|
||||
$(
|
||||
let mut cases = String::new();
|
||||
for _ in 0..NTESTS {
|
||||
let i1 = f64(rng);
|
||||
let i2 = f64(rng);
|
||||
for (i1, i2) in iproduct!(
|
||||
EDGE_CASES64.iter().cloned(),
|
||||
EDGE_CASES64.iter().cloned()
|
||||
).chain(EDGE_CASES64.iter().map(|i1| (*i1, f64(rng))))
|
||||
.chain(EDGE_CASES64.iter().map(|i2| (f64(&mut rng2), *i2)))
|
||||
.chain((0..NTESTS).map(|_| (f64(&mut rng3), f64(&mut rng3))))
|
||||
{
|
||||
let out = unsafe { $intr(i1, i2) };
|
||||
|
||||
let i1 = i1.to_bits();
|
||||
|
|
@ -485,12 +527,16 @@ macro_rules! f64f64f64_f64 {
|
|||
$(fn $intr(_: f64, _: f64, _: f64) -> f64;)*
|
||||
}
|
||||
|
||||
let mut rng2 = rng.clone();
|
||||
$(
|
||||
let mut cases = String::new();
|
||||
for _ in 0..NTESTS {
|
||||
let i1 = f64(rng);
|
||||
let i2 = f64(rng);
|
||||
let i3 = f64(rng);
|
||||
for (i1, i2, i3) in iproduct!(
|
||||
EDGE_CASES64.iter().cloned(),
|
||||
EDGE_CASES64.iter().cloned(),
|
||||
EDGE_CASES64.iter().cloned()
|
||||
).chain(EDGE_CASES64.iter().map(|i1| (*i1, f64(rng), f64(rng))))
|
||||
.chain((0..NTESTS).map(|_| (f64(&mut rng2), f64(&mut rng2), f64(&mut rng2))))
|
||||
{
|
||||
let out = unsafe { $intr(i1, i2, i3) };
|
||||
|
||||
let i1 = i1.to_bits();
|
||||
|
|
@ -565,10 +611,10 @@ macro_rules! f64i32_f64 {
|
|||
$(fn $intr(_: f64, _: i32) -> f64;)*
|
||||
}
|
||||
|
||||
let mut rng2 = rng.clone();
|
||||
$(
|
||||
let mut cases = String::new();
|
||||
for _ in 0..NTESTS {
|
||||
let i1 = f64(rng);
|
||||
for i1 in EDGE_CASES64.iter().cloned().chain((0..NTESTS).map(|_| f64(&mut rng2))) {
|
||||
let i2 = rng.gen_range(i16::MIN, i16::MAX);
|
||||
let out = unsafe { $intr(i1, i2 as i32) };
|
||||
|
||||
|
|
@ -687,7 +733,7 @@ f32f32_f32! {
|
|||
|
||||
// With signature `fn(f32, f32, f32) -> f32`
|
||||
f32f32f32_f32! {
|
||||
// fmaf,
|
||||
fmaf,
|
||||
}
|
||||
|
||||
// With signature `fn(f32, i32) -> f32`
|
||||
|
|
@ -699,11 +745,11 @@ f32i32_f32! {
|
|||
f64_f64! {
|
||||
acos,
|
||||
asin,
|
||||
// atan,
|
||||
atan,
|
||||
cbrt,
|
||||
ceil,
|
||||
cos,
|
||||
// cosh,
|
||||
cosh,
|
||||
exp,
|
||||
exp2,
|
||||
expm1,
|
||||
|
|
@ -717,7 +763,7 @@ f64_f64! {
|
|||
sinh,
|
||||
sqrt,
|
||||
tan,
|
||||
// tanh,
|
||||
tanh,
|
||||
trunc,
|
||||
fabs,
|
||||
}
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue