diff --git a/library/compiler-builtins/libm/CHANGELOG.md b/library/compiler-builtins/libm/CHANGELOG.md new file mode 100644 index 000000000000..c3e74a814338 --- /dev/null +++ b/library/compiler-builtins/libm/CHANGELOG.md @@ -0,0 +1,12 @@ +# Change Log + +All notable changes to this project will be documented in this file. +This project adheres to [Semantic Versioning](http://semver.org/). + +## [Unreleased] + +## v0.1.0 - 2018-07-13 + +- Initial release + +[Unreleased]: https://github.com/japaric/libm/compare/v0.1.0...HEAD diff --git a/library/compiler-builtins/libm/CONTRIBUTING.md b/library/compiler-builtins/libm/CONTRIBUTING.md new file mode 100644 index 000000000000..1f0a05a3f8d9 --- /dev/null +++ b/library/compiler-builtins/libm/CONTRIBUTING.md @@ -0,0 +1,87 @@ +# How to contribute + +- Pick your favorite math function from the [issue tracker]. +- Look for the C implementation of the function in the [MUSL source code][src]. +- Copy paste the C code into a Rust file in the `src/math` directory and adjust `src/math/mod.rs` + accordingly. Also, uncomment the corresponding trait method in `src/lib.rs`. +- Run `cargo watch check` and fix the compiler errors. +- Tweak the bottom of `test-generator/src/main.rs` to add your function to the test suite. +- If you can, run the full test suite locally (see the [testing](#testing) section below). If you + can't, no problem! Your PR will be fully tested automatically. Though you may still want to add + and run some unit tests. See the bottom of [`src/math/truncf.rs`] for an example of such tests; + you can run unit tests with the `cargo test --lib` command. +- Send us a pull request! +- :tada: + +[issue tracker]: https://github.com/japaric/libm/issues +[src]: https://git.musl-libc.org/cgit/musl/tree/src/math +[`src/math/truncf.rs`]: https://github.com/japaric/libm/blob/master/src/math/truncf.rs + +Check [PR #65] for an example. + +[PR #65]: https://github.com/japaric/libm/pull/65 + +## Tips and tricks + +- *IMPORTANT* The code in this crate will end up being used in the `core` crate so it can **not** + have any external dependencies (other than `core` itself). + +- Only use relative imports within the `math` directory / module, e.g. `use self::fabs::fabs` or +`use super::isnanf`. Absolute imports from core are OK, e.g. `use core::u64`. + +- To reinterpret a float as an integer use the `to_bits` method. The MUSL code uses the + `GET_FLOAT_WORD` macro, or a union, to do this operation. + +- To reinterpret an integer as a float use the `f32::from_bits` constructor. The MUSL code uses the + `SET_FLOAT_WORD` macro, or a union, to do this operation. + +- You may encounter weird literals like `0x1p127f` in the MUSL code. These are hexadecimal floating + point literals. Rust (the language) doesn't support these kind of literals. The best way I have + found to deal with these literals is to turn them into their integer representation using the + [`hexf!`] macro and then turn them back into floats. See below: + +[`hexf!`]: https://crates.io/crates/hexf + +``` rust +// Step 1: write a program to convert the float into its integer representation +#[macro_use] +extern crate hexf; + +fn main() { + println!("{:#x}", hexf32!("0x1.0p127").to_bits()); +} +``` + +``` console +$ # Step 2: run the program +$ cargo run +0x7f000000 +``` + +``` rust +// Step 3: copy paste the output into libm +let x1p127 = f32::from_bits(0x7f000000); // 0x1p127f === 2 ^ 12 +``` + +- Rust code panics on arithmetic overflows when not optimized. You may need to use the [`Wrapping`] + newtype to avoid this problem. + +[`Wrapping`]: https://doc.rust-lang.org/std/num/struct.Wrapping.html + +## Testing + +The test suite of this crate can only be run on x86_64 Linux systems using the following commands: + +``` console +$ # The test suite depends on the `cross` tool so install it if you don't have it +$ cargo install cross + +$ # and the `cross` tool requires docker to be running +$ systemctl start docker + +$ # execute the test suite for the x86_64 target +$ TARGET=x86_64-unknown-linux-gnu bash ci/script.sh + +$ # execute the test suite for the ARMv7 target +$ TARGET=armv7-unknown-linux-gnueabihf bash ci/script.sh +``` diff --git a/library/compiler-builtins/libm/Cargo.toml b/library/compiler-builtins/libm/Cargo.toml index e3498eed0381..5a1ae4a6c92b 100644 --- a/library/compiler-builtins/libm/Cargo.toml +++ b/library/compiler-builtins/libm/Cargo.toml @@ -1,7 +1,13 @@ [package] -name = "libm" -version = "0.1.0" authors = ["Jorge Aparicio "] +categories = ["no-std"] +description = "libm in pure Rust" +documentation = "https://docs.rs/libm" +keywords = ["libm", "math"] +license = "MIT OR Apache-2.0" +name = "libm" +repository = "https://github.com/japaric/libm" +version = "0.1.0" [workspace] members = ["cb", "test-generator"] \ No newline at end of file diff --git a/library/compiler-builtins/libm/README.md b/library/compiler-builtins/libm/README.md index cb29baf63a7b..02de9765a50b 100644 --- a/library/compiler-builtins/libm/README.md +++ b/library/compiler-builtins/libm/README.md @@ -6,62 +6,41 @@ A port of [MUSL]'s libm to Rust. ## Goals -The short term goal of this library is to enable math support (e.g. `sin`, `atan2`) for the -`wasm32-unknown-unknown` target. The longer term goal is to enable math support in the `core` crate. +The short term goal of this library is to [enable math support (e.g. `sin`, `atan2`) for the +`wasm32-unknown-unknown` target][wasm] (cf. [rust-lang-nursery/compiler-builtins][pr]). The longer +term goal is to enable [math support in the `core` crate][core]. -## Testing +[wasm]: https://github.com/japaric/libm/milestone/1 +[pr]: https://github.com/rust-lang-nursery/compiler-builtins/pull/248 +[core]: https://github.com/japaric/libm/milestone/2 -The test suite of this crate can only be run on x86_64 Linux systems. +## Already usable +This crate is [on crates.io] and can be used today in stable `#![no_std]` programs like this: + +[on crates.io]: https://crates.io/crates/libm + +``` rust +#![no_std] + +extern crate libm; + +use libm::F32Ext; // adds methods to `f32` + +fn foo(x: f32) { + let y = x.sqrt(); + let z = libm::truncf(x); +} ``` -$ # The test suite depends on the `cross` tool so install it if you don't have it -$ cargo install cross -$ # and the `cross` tool requires docker to be running -$ systemctl start docker +Not all the math functions are available at the moment. Check the [API docs] to learn what's +currently supported. -$ # execute the test suite for the x86_64 target -$ TARGET=x86_64-unknown-linux-gnu bash ci/script.sh - -$ # execute the test suite for the ARMv7 target -$ TARGET=armv7-unknown-linux-gnueabihf bash ci/script.sh -``` +[API docs]: https://docs.rs/libm ## Contributing -- Pick your favorite math function from the [issue tracker]. -- Look for the C implementation of the function in the [MUSL source code][src]. -- Copy paste the C code into a Rust file in the `src/math` directory and adjust `src/math/mod.rs` - accordingly. Also, uncomment the corresponding trait method in `src/lib.rs`. -- Run `cargo watch check` and fix the compiler errors. -- Tweak the bottom of `test-generator/src/main.rs` to add your function to the test suite. -- If you can, run the test suite locally. If you can't, no problem! Your PR will be tested - automatically. -- Send us a pull request! -- :tada: - -[issue tracker]: https://github.com/japaric/libm/issues -[src]: https://git.musl-libc.org/cgit/musl/tree/src/math - -Check [PR #2] for an example. - -[PR #2]: https://github.com/japaric/libm/pull/2 - -### Notes - -- Only use relative imports within the `math` directory / module, e.g. `use self::fabs::fabs` or -`use super::isnanf`. Absolute imports from core are OK, e.g. `use core::u64`. - -- To reinterpret a float as an integer use the `to_bits` method. The MUSL code uses the - `GET_FLOAT_WORD` macro, or a union, to do this operation. - -- To reinterpret an integer as a float use the `f32::from_bits` constructor. The MUSL code uses the - `SET_FLOAT_WORD` macro, or a union, to do this operation. - -- Rust code panics on arithmetic overflows when not optimized. You may need to use the [`Wrapping`] - newtype to avoid this problem. - -[`Wrapping`]: https://doc.rust-lang.org/std/num/struct.Wrapping.html +Please check [CONTRIBUTING.md](CONTRIBUTING.md) ## License diff --git a/library/compiler-builtins/libm/src/lib.rs b/library/compiler-builtins/libm/src/lib.rs index b87e82b2ac6a..b18687c1e333 100644 --- a/library/compiler-builtins/libm/src/lib.rs +++ b/library/compiler-builtins/libm/src/lib.rs @@ -32,9 +32,8 @@ pub fn _eq(a: u64, b: u64) -> bool { /// Math support for `f32` /// -/// NOTE this meant to be a closed extension trait. The only stable way to use this trait is to -/// import it to access its methods. -pub trait F32Ext { +/// This trait is sealed and cannot be implemented outside of `libm`. +pub trait F32Ext: private::Sealed { #[cfg(todo)] fn floor(self) -> Self; @@ -43,7 +42,6 @@ pub trait F32Ext { fn round(self) -> Self; - #[cfg(todo)] fn trunc(self) -> Self; #[cfg(todo)] @@ -70,16 +68,13 @@ pub trait F32Ext { fn sqrt(self) -> Self; - #[cfg(todo)] fn exp(self) -> Self; #[cfg(todo)] fn exp2(self) -> Self; - #[cfg(todo)] fn ln(self) -> Self; - #[cfg(todo)] fn log(self, base: Self) -> Self; #[cfg(todo)] @@ -91,7 +86,6 @@ pub trait F32Ext { #[cfg(todo)] fn cbrt(self) -> Self; - #[cfg(todo)] fn hypot(self, other: Self) -> Self; #[cfg(todo)] @@ -164,7 +158,6 @@ impl F32Ext for f32 { roundf(self) } - #[cfg(todo)] #[inline] fn trunc(self) -> Self { truncf(self) @@ -218,7 +211,6 @@ impl F32Ext for f32 { sqrtf(self) } - #[cfg(todo)] #[inline] fn exp(self) -> Self { expf(self) @@ -230,13 +222,11 @@ impl F32Ext for f32 { exp2f(self) } - #[cfg(todo)] #[inline] fn ln(self) -> Self { logf(self) } - #[cfg(todo)] #[inline] fn log(self, base: Self) -> Self { self.ln() / base.ln() @@ -260,7 +250,6 @@ impl F32Ext for f32 { cbrtf(self) } - #[cfg(todo)] #[inline] fn hypot(self, other: Self) -> Self { hypotf(self, other) @@ -364,12 +353,10 @@ impl F32Ext for f32 { } } -/// Math support for `f32` +/// Math support for `f64` /// -/// NOTE this meant to be a closed extension trait. The only stable way to use this trait is to -/// import it to access its methods. -pub trait F64Ext { - #[cfg(todo)] +/// This trait is sealed and cannot be implemented outside of `libm`. +pub trait F64Ext: private::Sealed { fn floor(self) -> Self; #[cfg(todo)] @@ -377,7 +364,6 @@ pub trait F64Ext { fn round(self) -> Self; - #[cfg(todo)] fn trunc(self) -> Self; #[cfg(todo)] @@ -403,7 +389,6 @@ pub trait F64Ext { #[cfg(todo)] fn powf(self, n: Self) -> Self; - #[cfg(todo)] fn sqrt(self) -> Self; #[cfg(todo)] @@ -427,7 +412,6 @@ pub trait F64Ext { #[cfg(todo)] fn cbrt(self) -> Self; - #[cfg(todo)] fn hypot(self, other: Self) -> Self; #[cfg(todo)] @@ -483,7 +467,6 @@ pub trait F64Ext { } impl F64Ext for f64 { - #[cfg(todo)] #[inline] fn floor(self) -> Self { floor(self) @@ -500,7 +483,6 @@ impl F64Ext for f64 { round(self) } - #[cfg(todo)] #[inline] fn trunc(self) -> Self { trunc(self) @@ -550,7 +532,6 @@ impl F64Ext for f64 { pow(self, n) } - #[cfg(todo)] #[inline] fn sqrt(self) -> Self { sqrt(self) @@ -598,7 +579,6 @@ impl F64Ext for f64 { cbrt(self) } - #[cfg(todo)] #[inline] fn hypot(self, other: Self) -> Self { hypot(self, other) @@ -701,3 +681,10 @@ impl F64Ext for f64 { 0.5 * ((2.0 * self) / (1.0 - self)).ln_1p() } } + +mod private { + pub trait Sealed {} + + impl Sealed for f32 {} + impl Sealed for f64 {} +} diff --git a/library/compiler-builtins/libm/src/math/expf.rs b/library/compiler-builtins/libm/src/math/expf.rs new file mode 100644 index 000000000000..1b645654e369 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/expf.rs @@ -0,0 +1,76 @@ +use super::scalbnf; + +const HALF : [f32; 2] = [0.5,-0.5]; +const LN2_HI : f32 = 6.9314575195e-01; /* 0x3f317200 */ +const LN2_LO : f32 = 1.4286067653e-06; /* 0x35bfbe8e */ +const INV_LN2 : f32 = 1.4426950216e+00; /* 0x3fb8aa3b */ +/* + * Domain [-0.34568, 0.34568], range ~[-4.278e-9, 4.447e-9]: + * |x*(exp(x)+1)/(exp(x)-1) - p(x)| < 2**-27.74 + */ +const P1 : f32 = 1.6666625440e-1; /* 0xaaaa8f.0p-26 */ +const P2 : f32 = -2.7667332906e-3; /* -0xb55215.0p-32 */ + +#[inline] +pub fn expf(mut x: f32) -> f32 { + let x1p127 = f32::from_bits(0x7f000000); // 0x1p127f === 2 ^ 127 + let x1p_126 = f32::from_bits(0x800000); // 0x1p-126f === 2 ^ -126 /*original 0x1p-149f ??????????? */ + + let mut hx = x.to_bits(); + let sign = (hx >> 31) as i32; /* sign bit of x */ + let signb : bool = sign != 0; + hx &= 0x7fffffff; /* high word of |x| */ + + /* special cases */ + if hx >= 0x42aeac50 { /* if |x| >= -87.33655f or NaN */ + if hx > 0x7f800000 {/* NaN */ + return x; + } + if (hx >= 0x42b17218) && (!signb) { /* x >= 88.722839f */ + /* overflow */ + x *= x1p127; + return x; + } + if signb { + /* underflow */ + force_eval!(-x1p_126/x); + if hx >= 0x42cff1b5 { /* x <= -103.972084f */ + return 0. + } + } + } + + /* argument reduction */ + let k : i32; + let hi : f32; + let lo : f32; + if hx > 0x3eb17218 { /* if |x| > 0.5 ln2 */ + if hx > 0x3f851592 { /* if |x| > 1.5 ln2 */ + k = (INV_LN2*x + HALF[sign as usize]) as i32; + } else { + k = 1 - sign - sign; + } + let kf = k as f32; + hi = x - kf*LN2_HI; /* k*ln2hi is exact here */ + lo = kf*LN2_LO; + x = hi - lo; + } else if hx > 0x39000000 { /* |x| > 2**-14 */ + k = 0; + hi = x; + lo = 0.; + } else { + /* raise inexact */ + force_eval!(x1p127 + x); + return 1. + x; + } + + /* x is now in primary range */ + let xx = x*x; + let c = x - xx*(P1+xx*P2); + let y = 1. + (x*c/(2.-c) - lo + hi); + if k == 0 { + y + } else { + scalbnf(y, k) + } +} diff --git a/library/compiler-builtins/libm/src/math/floor.rs b/library/compiler-builtins/libm/src/math/floor.rs new file mode 100644 index 000000000000..a5fb1757413b --- /dev/null +++ b/library/compiler-builtins/libm/src/math/floor.rs @@ -0,0 +1,29 @@ +use core::f64; + +const TOINT : f64 = 1. / f64::EPSILON; + +#[inline] +pub fn floor(x : f64) -> f64 { + let ui = x.to_bits(); + let e = ((ui >> 52) & 0x7ff) as i32; + + if (e >= 0x3ff+52) || (x == 0.) { + return x; + } + /* y = int(x) - x, where int(x) is an integer neighbor of x */ + let y = if (ui >> 63) != 0 { + x - TOINT + TOINT - x + } else { + x + TOINT - TOINT - x + }; + /* special case because of non-nearest rounding modes */ + if e <= 0x3ff-1 { + force_eval!(y); + return if (ui >> 63) != 0 { -1. } else { 0. }; + } + if y > 0. { + x + y - 1. + } else { + x + y + } +} diff --git a/library/compiler-builtins/libm/src/math/hypot.rs b/library/compiler-builtins/libm/src/math/hypot.rs new file mode 100644 index 000000000000..dcc17d914fb9 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/hypot.rs @@ -0,0 +1,74 @@ +use core::f64; + +use super::sqrt; + +const SPLIT: f64 = 134217728. + 1.; // 0x1p27 + 1 === (2 ^ 27) + 1 + +fn sq(x: f64) -> (f64, f64) { + let xh: f64; + let xl: f64; + let xc: f64; + + xc = x * SPLIT; + xh = x - xc + xc; + xl = x - xh; + let hi = x*x; + let lo = xh*xh - hi + 2.*xh*xl + xl*xl; + (hi, lo) +} + +#[inline] +pub fn hypot(mut x: f64, mut y: f64) -> f64 { + let x1p700 = f64::from_bits(0x6bb0000000000000); // 0x1p700 === 2 ^ 700 + let x1p_700 = f64::from_bits(0x1430000000000000); // 0x1p-700 === 2 ^ -700 + + let mut uxi = x.to_bits(); + let mut uyi = y.to_bits(); + let uti; + let ex: i64; + let ey: i64; + let mut z: f64; + + /* arrange |x| >= |y| */ + uxi &= -1i64 as u64 >> 1; + uyi &= -1i64 as u64 >> 1; + if uxi < uyi { + uti = uxi; + uxi = uyi; + uyi = uti; + } + + /* special cases */ + ex = (uxi>>52) as i64; + ey = (uyi>>52) as i64; + x = f64::from_bits(uxi); + y = f64::from_bits(uyi); + /* note: hypot(inf,nan) == inf */ + if ey == 0x7ff { + return y; + } + if ex == 0x7ff || uyi == 0 { + return x; + } + /* note: hypot(x,y) ~= x + y*y/x/2 with inexact for small y/x */ + /* 64 difference is enough for ld80 double_t */ + if ex - ey > 64 { + return x + y; + } + + /* precise sqrt argument in nearest rounding mode without overflow */ + /* xh*xh must not overflow and xl*xl must not underflow in sq */ + z = 1.; + if ex > 0x3ff+510 { + z = x1p700; + x *= x1p_700; + y *= x1p_700; + } else if ey < 0x3ff-450 { + z = x1p_700; + x *= x1p700; + y *= x1p700; + } + let (hx, lx) = sq(x); + let (hy, ly) = sq(y); + return z*sqrt(ly+lx+hy+hx); +} diff --git a/library/compiler-builtins/libm/src/math/hypotf.rs b/library/compiler-builtins/libm/src/math/hypotf.rs new file mode 100644 index 000000000000..146aab4e4464 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/hypotf.rs @@ -0,0 +1,43 @@ +use core::f32; + +use super::sqrtf; + +#[inline] +pub fn hypotf(mut x: f32, mut y: f32) -> f32 { + let x1p90 = f32::from_bits(0x6c800000); // 0x1p90f === 2 ^ 90 + let x1p_90 = f32::from_bits(0x12800000); // 0x1p-90f === 2 ^ -90 + + let mut uxi = x.to_bits(); + let mut uyi = y.to_bits(); + let uti; + let mut z: f32; + + uxi &= -1i32 as u32 >> 1; + uyi &= -1i32 as u32 >> 1; + if uxi < uyi { + uti = uxi; + uxi = uyi; + uyi = uti; + } + + x = f32::from_bits(uxi); + y = f32::from_bits(uyi); + if uyi == 0xff<<23 { + return y; + } + if uxi >= 0xff<<23 || uyi == 0 || uxi - uyi >= 25<<23 { + return x + y; + } + + z = 1.; + if uxi >= (0x7f+60)<<23 { + z = x1p90; + x *= x1p_90; + y *= x1p_90; + } else if uyi < (0x7f-60)<<23 { + z = x1p_90; + x *= x1p90; + y *= x1p90; + } + z*sqrtf((x as f64 * x as f64 + y as f64 * y as f64) as f32) +} diff --git a/library/compiler-builtins/libm/src/math/logf.rs b/library/compiler-builtins/libm/src/math/logf.rs new file mode 100644 index 000000000000..76b4ede19514 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/logf.rs @@ -0,0 +1,49 @@ +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 */ + +#[inline] +pub fn logf(mut x: f32) -> f32 { + let x1p25 = f32::from_bits(0x4c000000); // 0x1p25f === 2 ^ 25 + + let mut ix = x.to_bits(); + let mut k = 0i32; + + if (ix < 0x00800000) || ((ix>>31) != 0) { /* x < 2**-126 */ + if ix<<1 == 0 { + return -1./(x*x); /* log(+-0)=-inf */ + } + if (ix>>31) != 0 { + return (x-x)/0.; /* log(-#) = NaN */ + } + /* subnormal number, scale up x */ + k -= 25; + x *= x1p25; + ix = x.to_bits(); + } else if ix >= 0x7f800000 { + return x; + } else if ix == 0x3f800000 { + return 0.; + } + + /* reduce x into [sqrt(2)/2, sqrt(2)] */ + ix += 0x3f800000 - 0x3f3504f3; + k += ((ix>>23) as i32) - 0x7f; + ix = (ix & 0x007fffff) + 0x3f3504f3; + x = f32::from_bits(ix); + + let f = x - 1.; + let s = f/(2. + f); + let z = s*s; + let w = z*z; + let t1 = w*(LG2+w*LG4); + let t2 = z*(LG1+w*LG3); + let r = t2 + t1; + let hfsq = 0.5*f*f; + let dk = k as f32; + s*(hfsq+r) + dk*LN2_LO - hfsq + f + dk*LN2_HI +} diff --git a/library/compiler-builtins/libm/src/math/mod.rs b/library/compiler-builtins/libm/src/math/mod.rs index 8032bdabf9ac..93a067102119 100644 --- a/library/compiler-builtins/libm/src/math/mod.rs +++ b/library/compiler-builtins/libm/src/math/mod.rs @@ -12,17 +12,39 @@ mod fmodf; mod powf; mod round; mod roundf; +mod scalbn; mod scalbnf; +mod sqrt; mod sqrtf; +mod logf; +mod expf; +mod floor; +mod trunc; +mod truncf; +mod hypot; +mod hypotf; -pub use self::fabs::fabs; -pub use self::fabsf::fabsf; -pub use self::fmodf::fmodf; -pub use self::powf::powf; -pub use self::round::round; -pub use self::roundf::roundf; -pub use self::scalbnf::scalbnf; -pub use self::sqrtf::sqrtf; +//mod service; + +pub use self::{ + fabs::fabs, + fabsf::fabsf, + fmodf::fmodf, + powf::powf, + round::round, + roundf::roundf, + scalbn::scalbn, + scalbnf::scalbnf, + sqrt::sqrt, + sqrtf::sqrtf, + logf::logf, + expf::expf, + floor::floor, + trunc::trunc, + truncf::truncf, + hypot::hypot, + hypotf::hypotf, +}; fn isnanf(x: f32) -> bool { x.to_bits() & 0x7fffffff > 0x7f800000 diff --git a/library/compiler-builtins/libm/src/math/scalbn.rs b/library/compiler-builtins/libm/src/math/scalbn.rs new file mode 100644 index 000000000000..76e06f03ece2 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/scalbn.rs @@ -0,0 +1,33 @@ +#[inline] +pub fn scalbn(x : f64, mut n: i32) -> f64 { + let x1p1023 = f64::from_bits(0x7fe0000000000000); // 0x1p1023 === 2 ^ 1023 + let x1p53 = f64::from_bits(0x4340000000000000); // 0x1p53 === 2 ^ 53 + let x1p_1022 = f64::from_bits(0x0010000000000000); // 0x1p-1022 === 2 ^ (-1022) + + let mut y = x; + + if n > 1023 { + y *= x1p1023; + n -= 1023; + if n > 1023 { + y *= x1p1023; + n -= 1023; + if n > 1023 { + n = 1023; + } + } + } else if n < -1022 { + /* make sure final n < -53 to avoid double + rounding in the subnormal range */ + y *= x1p_1022 * x1p53; + n += 1022 - 53; + if n < -1022 { + y *= x1p_1022 * x1p53; + n += 1022 - 53; + if n < -1022 { + n = -1022; + } + } + } + y*f64::from_bits(((0x3ff+n) as u64)<<52) +} diff --git a/library/compiler-builtins/libm/src/math/scalbnf.rs b/library/compiler-builtins/libm/src/math/scalbnf.rs index 2ae8bf31b772..31f93d323c5c 100644 --- a/library/compiler-builtins/libm/src/math/scalbnf.rs +++ b/library/compiler-builtins/libm/src/math/scalbnf.rs @@ -1,35 +1,29 @@ #[inline] -pub fn scalbnf(mut x: f32, mut n: i32) -> f32 { - let x1p127 = f32::from_bits(0x7f000000); // 0x1p127f === 2 ^ 127 - let x1p_126 = f32::from_bits(0x800000); // 0x1p-126f === 2 ^ -126 - let x1p24 = f32::from_bits(0x4b800000); // 0x1p24f === 2 ^ 24 - - let mut y: f32 = x; - +pub fn scalbnf(mut x: f32, mut n : i32) -> f32 { + let x1p127 = f32::from_bits(0x7f000000); // 0x1p127f === 2 ^ 127 + let x1p_126 = f32::from_bits(0x800000); // 0x1p-126f === 2 ^ -126 + let x1p24 = f32::from_bits(0x4b800000); // 0x1p24f === 2 ^ 24 + if n > 127 { - y *= x1p127; + x *= x1p127; n -= 127; if n > 127 { - y *= x1p127; + x *= x1p127; n -= 127; if n > 127 { n = 127; } } } else if n < -126 { - y *= x1p_126; - y *= x1p24; + x *= x1p_126 * x1p24; n += 126 - 24; if n < -126 { - y *= x1p_126; - y *= x1p24; + x *= x1p_126 * x1p24; n += 126 - 24; if n < -126 { n = -126; } } } - - x = y * f32::from_bits((0x7f + n as u32) << 23); - x + x * f32::from_bits(((0x7f+n) as u32)<<23) } diff --git a/library/compiler-builtins/libm/src/math/sqrt.rs b/library/compiler-builtins/libm/src/math/sqrt.rs new file mode 100644 index 000000000000..49fc58fff03c --- /dev/null +++ b/library/compiler-builtins/libm/src/math/sqrt.rs @@ -0,0 +1,129 @@ +use core::f64; + +const TINY: f64 = 1.0e-300; + +#[inline] +pub fn sqrt(x: f64) -> f64 { + let mut z: f64; + let sign: u32 = 0x80000000; + let mut ix0: i32; + let mut s0: i32; + let mut q: i32; + let mut m: i32; + let mut t: i32; + let mut i: i32; + let mut r: u32; + let mut t1: u32; + let mut s1: u32; + let mut ix1: u32; + let mut q1: u32; + + ix0 = (x.to_bits() >> 32) as i32; + ix1 = x.to_bits() as u32; + + /* take care of Inf and NaN */ + if (ix0&0x7ff00000) == 0x7ff00000 { + return x*x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */ + } + /* take care of zero */ + if ix0 <= 0 { + if ((ix0&!(sign as i32))|ix1 as i32) == 0 { + return x; /* sqrt(+-0) = +-0 */ + } + if ix0 < 0 { + return (x - x) / (x - x); /* sqrt(-ve) = sNaN */ + } + } + /* normalize x */ + m = ix0>>20; + if m == 0 { /* subnormal x */ + while ix0 == 0 { + m -= 21; + ix0 |= (ix1>>11) as i32; + ix1 <<= 21; + } + i=0; + while (ix0&0x00100000) == 0 { + i += 1; + ix0 <<= 1; + } + m -= i - 1; + ix0 |= (ix1>>(32-i)) as i32; + ix1 <<= i; + } + m -= 1023; /* unbias exponent */ + ix0 = (ix0&0x000fffff)|0x00100000; + if (m & 1) == 1 { /* odd m, double x to make it even */ + ix0 += ix0 + ((ix1&sign)>>31) as i32; + ix1 += ix1; + } + m >>= 1; /* m = [m/2] */ + + /* generate sqrt(x) bit by bit */ + ix0 += ix0 + ((ix1&sign)>>31) as i32; + ix1 += ix1; + q = 0; /* [q,q1] = sqrt(x) */ + q1 = 0; + s0 = 0; + s1 = 0; + r = 0x00200000; /* r = moving bit from right to left */ + + while r != 0 { + t = s0 + r as i32; + if t <= ix0 { + s0 = t + r as i32; + ix0 -= t; + q += r as i32; + } + ix0 += ix0 + ((ix1&sign)>>31) as i32; + ix1 += ix1; + r >>= 1; + } + + r = sign; + while r != 0 { + t1 = s1 + r; + t = s0; + if t < ix0 || (t == ix0 && t1 <= ix1) { + s1 = t1 + r; + if (t1&sign) == sign && (s1&sign) == 0 { + s0 += 1; + } + ix0 -= t; + if ix1 < t1 { + ix0 -= 1; + } + ix1 -= t1; + q1 += r; + } + ix0 += ix0 + ((ix1&sign)>>31) as i32; + ix1 += ix1; + r >>= 1; + } + + /* use floating add to find out rounding direction */ + if (ix0 as u32|ix1) != 0 { + z = 1.0 - TINY; /* raise inexact flag */ + if z >= 1.0 { + z = 1.0 + TINY; + if q1 == 0xffffffff { + q1 = 0; + q+=1; + } else if z > 1.0 { + if q1 == 0xfffffffe { + q += 1; + } + q1 += 2; + } else { + q1 += q1 & 1; + } + } + } + ix0 = (q>>1) + 0x3fe00000; + ix1 = q1>>1; + if (q&1) == 1 { + ix1 |= sign; + } + ix0 += m << 20; + f64::from_bits((ix0 as u64) << 32 | ix1 as u64) +} diff --git a/library/compiler-builtins/libm/src/math/trunc.rs b/library/compiler-builtins/libm/src/math/trunc.rs new file mode 100644 index 000000000000..6bea67cbc165 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/trunc.rs @@ -0,0 +1,32 @@ +use core::f64; + +#[inline] +pub fn trunc(x: f64) -> f64 { + let x1p120 = f64::from_bits(0x4770000000000000); // 0x1p120f === 2 ^ 120 + + let mut i: u64 = x.to_bits(); + let mut e: i64 = (i >> 52 & 0x7ff) as i64 - 0x3ff + 12; + let m: u64; + + if e >= 52 + 12 { + return x; + } + if e < 12 { + e = 1; + } + m = -1i64 as u64 >> e; + if (i & m) == 0 { + return x; + } + force_eval!(x + x1p120); + i &= !m; + f64::from_bits(i) +} + +#[cfg(test)] +mod tests { + #[test] + fn sanity_check() { + assert_eq!(super::trunc(1.1), 1.0); + } +} diff --git a/library/compiler-builtins/libm/src/math/truncf.rs b/library/compiler-builtins/libm/src/math/truncf.rs new file mode 100644 index 000000000000..9d42620d9666 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/truncf.rs @@ -0,0 +1,32 @@ +use core::f32; + +#[inline] +pub fn truncf(x: f32) -> f32 { + let x1p120 = f32::from_bits(0x7b800000); // 0x1p120f === 2 ^ 120 + + let mut i: u32 = x.to_bits(); + let mut e: i32 = (i >> 23 & 0xff) as i32 - 0x7f + 9; + let m: u32; + + if e >= 23 + 9 { + return x; + } + if e < 9 { + e = 1; + } + m = -1i32 as u32 >> e; + if (i & m) == 0 { + return x; + } + force_eval!(x + x1p120); + i &= !m; + f32::from_bits(i) +} + +#[cfg(test)] +mod tests { + #[test] + fn sanity_check() { + assert_eq!(super::truncf(1.1), 1.0); + } +} diff --git a/library/compiler-builtins/libm/test-generator/src/main.rs b/library/compiler-builtins/libm/test-generator/src/main.rs index 2f78a8ed38da..6474ab6f475c 100644 --- a/library/compiler-builtins/libm/test-generator/src/main.rs +++ b/library/compiler-builtins/libm/test-generator/src/main.rs @@ -62,6 +62,8 @@ macro_rules! f32_f32 { extern crate libm; + use std::panic; + #[test] fn {0}() {{ const CASES: &[(u32, u32)] = &[ @@ -71,16 +73,23 @@ macro_rules! f32_f32 { for case in CASES {{ let (inp, expected) = *case; - let outf = libm::{0}(f32::from_bits(inp)); - let outi = outf.to_bits(); + if let Ok(outf) = + panic::catch_unwind(|| libm::{0}(f32::from_bits(inp))) + {{ + let outi = outf.to_bits(); - if !((outf.is_nan() && f32::from_bits(expected).is_nan()) || - libm::_eqf(outi, expected)) {{ + if !((outf.is_nan() && f32::from_bits(expected).is_nan()) + || libm::_eqf(outi, expected)) + {{ + panic!( + \"input: {{}}, output: {{}}, expected: {{}}\", + inp, outi, expected, + ); + }} + }} else {{ panic!( - \"input: {{}}, output: {{}}, expected: {{}}\", - inp, - outi, - expected, + \"input: {{}}, output: PANIC, expected: {{}}\", + inp, expected, ); }} }} @@ -124,6 +133,8 @@ macro_rules! f32f32_f32 { extern crate libm; + use std::panic; + #[test] fn {0}() {{ const CASES: &[((u32, u32), u32)] = &[ @@ -133,15 +144,25 @@ macro_rules! f32f32_f32 { for case in CASES {{ let ((i1, i2), expected) = *case; - let outf = libm::{0}(f32::from_bits(i1), f32::from_bits(i2)); - let outi = outf.to_bits(); + if let Ok(outf) = panic::catch_unwind(|| {{ + libm::{0}(f32::from_bits(i1), f32::from_bits(i2)) + }}) {{ + let outi = outf.to_bits(); - if !((outf.is_nan() && f32::from_bits(expected).is_nan()) || - libm::_eqf(outi, expected)) {{ + if !((outf.is_nan() && f32::from_bits(expected).is_nan()) + || libm::_eqf(outi, expected)) + {{ + panic!( + \"input: {{:?}}, output: {{}}, expected: {{}}\", + (i1, i2), + outi, + expected, + ); + }} + }} else {{ panic!( - \"input: {{:?}}, output: {{}}, expected: {{}}\", + \"input: {{:?}}, output: PANIC, expected: {{}}\", (i1, i2), - outi, expected, ); }} @@ -188,6 +209,8 @@ macro_rules! f32f32f32_f32 { extern crate libm; + use std::panic; + #[test] fn {0}() {{ const CASES: &[((u32, u32, u32), u32)] = &[ @@ -197,19 +220,29 @@ macro_rules! f32f32f32_f32 { for case in CASES {{ let ((i1, i2, i3), expected) = *case; - let outf = libm::{0}( - f32::from_bits(i1), - f32::from_bits(i2), - f32::from_bits(i3), - ); - let outi = outf.to_bits(); + if let Ok(outf) = panic::catch_unwind(|| {{ + libm::{0}( + f32::from_bits(i1), + f32::from_bits(i2), + f32::from_bits(i3), + ) + }}) {{ + let outi = outf.to_bits(); - if !((outf.is_nan() && f32::from_bits(expected).is_nan()) || - libm::_eqf(outi, expected)) {{ + if !((outf.is_nan() && f32::from_bits(expected).is_nan()) + || libm::_eqf(outi, expected)) + {{ + panic!( + \"input: {{:?}}, output: {{}}, expected: {{}}\", + (i1, i2, i3), + outi, + expected, + ); + }} + }} else {{ panic!( - \"input: {{:?}}, output: {{}}, expected: {{}}\", + \"input: {{:?}}, output: PANIC, expected: {{}}\", (i1, i2, i3), - outi, expected, ); }} @@ -253,6 +286,8 @@ macro_rules! f32i32_f32 { extern crate libm; + use std::panic; + #[test] fn {0}() {{ const CASES: &[((u32, i16), u32)] = &[ @@ -262,15 +297,25 @@ macro_rules! f32i32_f32 { for case in CASES {{ let ((i1, i2), expected) = *case; - let outf = libm::{0}(f32::from_bits(i1), i2 as i32); - let outi = outf.to_bits(); + if let Ok(outf) = panic::catch_unwind(|| {{ + libm::{0}(f32::from_bits(i1), i2 as i32) + }}) {{ + let outi = outf.to_bits(); - if !((outf.is_nan() && f32::from_bits(expected).is_nan()) || - libm::_eqf(outi, expected)) {{ + if !((outf.is_nan() && f32::from_bits(expected).is_nan()) + || libm::_eqf(outi, expected)) + {{ + panic!( + \"input: {{:?}}, output: {{}}, expected: {{}}\", + (i1, i2), + outi, + expected, + ); + }} + }} else {{ panic!( - \"input: {{:?}}, output: {{}}, expected: {{}}\", + \"input: {{:?}}, output: PANIC, expected: {{}}\", (i1, i2), - outi, expected, ); }} @@ -314,6 +359,8 @@ macro_rules! f64_f64 { extern crate libm; + use std::panic; + #[test] fn {0}() {{ const CASES: &[(u64, u64)] = &[ @@ -323,15 +370,25 @@ macro_rules! f64_f64 { for case in CASES {{ let (inp, expected) = *case; - let outf = libm::{0}(f64::from_bits(inp)); - let outi = outf.to_bits(); + if let Ok(outf) = panic::catch_unwind(|| {{ + libm::{0}(f64::from_bits(inp)) + }}) {{ + let outi = outf.to_bits(); - if !((outf.is_nan() && f64::from_bits(expected).is_nan()) || - libm::_eq(outi, expected)) {{ + if !((outf.is_nan() && f64::from_bits(expected).is_nan()) + || libm::_eq(outi, expected)) + {{ + panic!( + \"input: {{}}, output: {{}}, expected: {{}}\", + inp, + outi, + expected, + ); + }} + }} else {{ panic!( - \"input: {{}}, output: {{}}, expected: {{}}\", + \"input: {{}}, output: PANIC, expected: {{}}\", inp, - outi, expected, ); }} @@ -376,6 +433,8 @@ macro_rules! f64f64_f64 { extern crate libm; + use std::panic; + #[test] fn {0}() {{ const CASES: &[((u64, u64), u64)] = &[ @@ -385,15 +444,24 @@ macro_rules! f64f64_f64 { for case in CASES {{ let ((i1, i2), expected) = *case; - let outf = libm::{0}(f64::from_bits(i1), f64::from_bits(i2)); - let outi = outf.to_bits(); + if let Ok(outf) = panic::catch_unwind(|| {{ + libm::{0}(f64::from_bits(i1), f64::from_bits(i2)) + }}) {{ + let outi = outf.to_bits(); - if !((outf.is_nan() && f64::from_bits(expected).is_nan()) || - libm::_eq(outi, expected)) {{ + if !((outf.is_nan() && f64::from_bits(expected).is_nan()) || + libm::_eq(outi, expected)) {{ + panic!( + \"input: {{:?}}, output: {{}}, expected: {{}}\", + (i1, i2), + outi, + expected, + ); + }} + }} else {{ panic!( - \"input: {{:?}}, output: {{}}, expected: {{}}\", + \"input: {{:?}}, output: PANIC, expected: {{}}\", (i1, i2), - outi, expected, ); }} @@ -440,6 +508,8 @@ macro_rules! f64f64f64_f64 { extern crate libm; + use std::panic; + #[test] fn {0}() {{ const CASES: &[((u64, u64, u64), u64)] = &[ @@ -449,19 +519,29 @@ macro_rules! f64f64f64_f64 { for case in CASES {{ let ((i1, i2, i3), expected) = *case; - let outf = libm::{0}( - f64::from_bits(i1), - f64::from_bits(i2), - f64::from_bits(i3), - ); - let outi = outf.to_bits(); + if let Ok(outf) = panic::catch_unwind(|| {{ + libm::{0}( + f64::from_bits(i1), + f64::from_bits(i2), + f64::from_bits(i3), + ) + }}) {{ + let outi = outf.to_bits(); - if !((outf.is_nan() && f64::from_bits(expected).is_nan()) || - libm::_eq(outi, expected)) {{ + if !((outf.is_nan() && f64::from_bits(expected).is_nan()) + || libm::_eq(outi, expected)) + {{ + panic!( + \"input: {{:?}}, output: {{}}, expected: {{}}\", + (i1, i2, i3), + outi, + expected, + ); + }} + }} else {{ panic!( - \"input: {{:?}}, output: {{}}, expected: {{}}\", + \"input: {{:?}}, output: PANIC, expected: {{}}\", (i1, i2, i3), - outi, expected, ); }} @@ -505,6 +585,8 @@ macro_rules! f64i32_f64 { extern crate libm; + use std::panic; + #[test] fn {0}() {{ const CASES: &[((u64, i16), u64)] = &[ @@ -514,15 +596,24 @@ macro_rules! f64i32_f64 { for case in CASES {{ let ((i1, i2), expected) = *case; - let outf = libm::{0}(f64::from_bits(i1), i2 as i32); - let outi = outf.to_bits(); + if let Ok(outf) = panic::catch_unwind(|| {{ + libm::{0}(f64::from_bits(i1), i2 as i32) + }}) {{ + let outi = outf.to_bits(); - if !((outf.is_nan() && f64::from_bits(expected).is_nan()) || - libm::_eq(outi, expected)) {{ + if !((outf.is_nan() && f64::from_bits(expected).is_nan()) || + libm::_eq(outi, expected)) {{ + panic!( + \"input: {{:?}}, output: {{}}, expected: {{}}\", + (i1, i2), + outi, + expected, + ); + }} + }} else {{ panic!( - \"input: {{:?}}, output: {{}}, expected: {{}}\", + \"input: {{:?}}, output: PANIC, expected: {{}}\", (i1, i2), - outi, expected, ); }} @@ -562,7 +653,7 @@ fn main() -> Result<(), Box> { f32_f32! { // acosf, // floorf, - // truncf + truncf, // asinf, // atanf, // cbrtf, @@ -570,10 +661,11 @@ f32_f32! { // cosf, // coshf, // exp2f, - // expf, + expf, // fdimf, // log10f, // log2f, + logf, roundf, // sinf, // sinhf, @@ -586,7 +678,7 @@ f32_f32! { // With signature `fn(f32, f32) -> f32` f32f32_f32! { // atan2f, - // hypotf, + hypotf, fmodf, powf, } @@ -613,7 +705,7 @@ f64_f64! { // exp, // exp2, // expm1, - // floor, + floor, // log, // log10, // log1p, @@ -621,10 +713,10 @@ f64_f64! { round, // sin, // sinh, - // sqrt, + sqrt, // tan, // tanh, - // trunc, + trunc, fabs, } @@ -633,7 +725,7 @@ f64f64_f64! { // atan2, // fdim, // fmod, - // hypot, + hypot, // pow, } @@ -644,5 +736,5 @@ f64f64f64_f64! { // With signature `fn(f64, i32) -> f64` f64i32_f64! { - // scalbn, + scalbn, }