diff --git a/library/compiler-builtins/libm/.travis.yml b/library/compiler-builtins/libm/.travis.yml index 0c0832793474..2bbb60fc64d1 100644 --- a/library/compiler-builtins/libm/.travis.yml +++ b/library/compiler-builtins/libm/.travis.yml @@ -15,6 +15,8 @@ matrix: - env: TARGET=powerpc64-unknown-linux-gnu - env: TARGET=powerpc64le-unknown-linux-gnu - env: TARGET=x86_64-unknown-linux-gnu + - env: TARGET=cargo-fmt + rust: beta before_install: set -e 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..6f8e984f3004 --- /dev/null +++ b/library/compiler-builtins/libm/CONTRIBUTING.md @@ -0,0 +1,88 @@ +# 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! Make sure to run `cargo fmt` on your code before sending the PR. Also + include "closes #42" in the PR description to close the corresponding issue. +- :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/ci/install.sh b/library/compiler-builtins/libm/ci/install.sh index efdbb06c206e..4d9552d231bc 100644 --- a/library/compiler-builtins/libm/ci/install.sh +++ b/library/compiler-builtins/libm/ci/install.sh @@ -1,6 +1,11 @@ set -euxo pipefail main() { + if [ $TARGET = cargo-fmt ]; then + rustup component add rustfmt-preview + return + fi + if ! hash cross >/dev/null 2>&1; then cargo install cross fi diff --git a/library/compiler-builtins/libm/ci/script.sh b/library/compiler-builtins/libm/ci/script.sh index 3db6bfeb12d9..f2a294b488f5 100644 --- a/library/compiler-builtins/libm/ci/script.sh +++ b/library/compiler-builtins/libm/ci/script.sh @@ -1,6 +1,11 @@ set -euxo pipefail main() { + if [ $TARGET = cargo-fmt ]; then + cargo fmt -- --check + return + fi + # quick check cargo check diff --git a/library/compiler-builtins/libm/src/lib.rs b/library/compiler-builtins/libm/src/lib.rs index 9f4365712d65..437a95fc3697 100644 --- a/library/compiler-builtins/libm/src/lib.rs +++ b/library/compiler-builtins/libm/src/lib.rs @@ -32,16 +32,12 @@ 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 { - #[cfg(todo)] +/// This trait is sealed and cannot be implemented outside of `libm`. +pub trait F32Ext: private::Sealed { fn floor(self) -> Self; - #[cfg(todo)] fn ceil(self) -> Self; - #[cfg(todo)] fn round(self) -> Self; fn trunc(self) -> Self; @@ -79,16 +75,13 @@ pub trait F32Ext { fn log(self, base: Self) -> Self; - #[cfg(todo)] fn log2(self) -> Self; - #[cfg(todo)] fn log10(self) -> Self; #[cfg(todo)] fn cbrt(self) -> Self; - #[cfg(todo)] fn hypot(self, other: Self) -> Self; #[cfg(todo)] @@ -143,19 +136,16 @@ pub trait F32Ext { } impl F32Ext for f32 { - #[cfg(todo)] #[inline] fn floor(self) -> Self { floorf(self) } - #[cfg(todo)] #[inline] fn ceil(self) -> Self { ceilf(self) } - #[cfg(todo)] #[inline] fn round(self) -> Self { roundf(self) @@ -235,13 +225,11 @@ impl F32Ext for f32 { self.ln() / base.ln() } - #[cfg(todo)] #[inline] fn log2(self) -> Self { log2f(self) } - #[cfg(todo)] #[inline] fn log10(self) -> Self { log10f(self) @@ -253,7 +241,6 @@ impl F32Ext for f32 { cbrtf(self) } - #[cfg(todo)] #[inline] fn hypot(self, other: Self) -> Self { hypotf(self, other) @@ -356,11 +343,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 { +/// This trait is sealed and cannot be implemented outside of `libm`. +pub trait F64Ext: private::Sealed { fn floor(self) -> Self; #[cfg(todo)] @@ -393,7 +379,6 @@ pub trait F64Ext { #[cfg(todo)] fn powf(self, n: Self) -> Self; - #[cfg(todo)] fn sqrt(self) -> Self; #[cfg(todo)] @@ -408,16 +393,13 @@ pub trait F64Ext { #[cfg(todo)] fn log(self, base: Self) -> Self; - #[cfg(todo)] fn log2(self) -> Self; - #[cfg(todo)] fn log10(self) -> Self; #[cfg(todo)] fn cbrt(self) -> Self; - #[cfg(todo)] fn hypot(self, other: Self) -> Self; #[cfg(todo)] @@ -538,7 +520,6 @@ impl F64Ext for f64 { pow(self, n) } - #[cfg(todo)] #[inline] fn sqrt(self) -> Self { sqrt(self) @@ -568,13 +549,11 @@ impl F64Ext for f64 { self.ln() / base.ln() } - #[cfg(todo)] #[inline] fn log2(self) -> Self { log2(self) } - #[cfg(todo)] #[inline] fn log10(self) -> Self { log10(self) @@ -586,7 +565,6 @@ impl F64Ext for f64 { cbrt(self) } - #[cfg(todo)] #[inline] fn hypot(self, other: Self) -> Self { hypot(self, other) @@ -689,3 +667,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/ceilf.rs b/library/compiler-builtins/libm/src/math/ceilf.rs new file mode 100644 index 000000000000..b4f58bfb86f3 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/ceilf.rs @@ -0,0 +1,29 @@ +use core::f32; + +pub fn ceilf(x: f32) -> f32 { + let mut ui = x.to_bits(); + let e = (((ui >> 23) & 0xff) - 0x7f) as i32; + + if e >= 23 { + return x; + } + if e >= 0 { + let m = 0x007fffff >> e; + if (ui & m) == 0 { + return x; + } + force_eval!(x + f32::from_bits(0x7b800000)); + if ui >> 31 == 0 { + ui += m; + } + ui &= !m; + } else { + force_eval!(x + f32::from_bits(0x7b800000)); + if ui >> 31 != 0 { + return -0.0; + } else if ui << 1 != 0 { + return 1.0; + } + } + return f32::from_bits(ui); +} diff --git a/library/compiler-builtins/libm/src/math/expf.rs b/library/compiler-builtins/libm/src/math/expf.rs index 1b645654e369..cffb557719f0 100644 --- a/library/compiler-builtins/libm/src/math/expf.rs +++ b/library/compiler-builtins/libm/src/math/expf.rs @@ -1,60 +1,67 @@ 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 */ +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 */ +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 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| */ - + 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 */ + if hx >= 0x42aeac50 { + /* if |x| >= -87.33655f or NaN */ + if hx > 0x7f800000 { + /* NaN */ return x; } - if (hx >= 0x42b17218) && (!signb) { /* x >= 88.722839f */ + 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. + 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; + 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; + hi = x - kf * LN2_HI; /* k*ln2hi is exact here */ + lo = kf * LN2_LO; x = hi - lo; - } else if hx > 0x39000000 { /* |x| > 2**-14 */ + } else if hx > 0x39000000 { + /* |x| > 2**-14 */ k = 0; hi = x; lo = 0.; @@ -63,11 +70,11 @@ pub fn expf(mut x: f32) -> f32 { 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); + 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 { diff --git a/library/compiler-builtins/libm/src/math/floor.rs b/library/compiler-builtins/libm/src/math/floor.rs index a5fb1757413b..997865d39e93 100644 --- a/library/compiler-builtins/libm/src/math/floor.rs +++ b/library/compiler-builtins/libm/src/math/floor.rs @@ -1,27 +1,27 @@ use core::f64; -const TOINT : f64 = 1. / f64::EPSILON; +const TOINT: f64 = 1. / f64::EPSILON; #[inline] -pub fn floor(x : f64) -> f64 { +pub fn floor(x: f64) -> f64 { let ui = x.to_bits(); - let e = ((ui >> 52) & 0x7ff) as i32; + let e = ((ui >> 52) & 0x7ff) as i32; - if (e >= 0x3ff+52) || (x == 0.) { - return x; + 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 + /* 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. { + /* 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/floorf.rs b/library/compiler-builtins/libm/src/math/floorf.rs new file mode 100644 index 000000000000..9c263b51828b --- /dev/null +++ b/library/compiler-builtins/libm/src/math/floorf.rs @@ -0,0 +1,30 @@ +use core::f32; + +#[inline] +pub fn floorf(x: f32) -> f32 { + let mut ui = x.to_bits(); + let e = (((ui >> 23) & 0xff) - 0x7f) as i32; + + if e >= 23 { + return x; + } + if e >= 0 { + let m: u32 = 0x007fffff >> e; + if (ui & m) == 0 { + return x; + } + force_eval!(x + f32::from_bits(0x7b800000)); + if ui >> 31 != 0 { + ui += m; + } + ui &= !m; + } else { + force_eval!(x + f32::from_bits(0x7b800000)); + if ui >> 31 == 0 { + ui = 0; + } else if ui << 1 != 0 { + return -1.0; + } + } + return f32::from_bits(ui); +} 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..7ad1baf796c8 --- /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..d59710ada395 --- /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/log10.rs b/library/compiler-builtins/libm/src/math/log10.rs new file mode 100644 index 000000000000..137d8966f521 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/log10.rs @@ -0,0 +1,98 @@ +use core::f64; + +const IVLN10HI: f64 = 4.34294481878168880939e-01; /* 0x3fdbcb7b, 0x15200000 */ +const IVLN10LO: f64 = 2.50829467116452752298e-11; /* 0x3dbb9438, 0xca9aadd5 */ +const LOG10_2HI: f64 = 3.01029995663611771306e-01; /* 0x3FD34413, 0x509F6000 */ +const LOG10_2LO: f64 = 3.69423907715893078616e-13; /* 0x3D59FEF3, 0x11F12B36 */ +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 log10(mut x: f64) -> f64 { + let x1p54 = f64::from_bits(0x4350000000000000); // 0x1p54 === 2 ^ 54 + + let mut ui: u64 = x.to_bits(); + let hfsq: f64; + let f: f64; + let s: f64; + let z: f64; + let r: f64; + let mut w: f64; + let t1: f64; + let t2: f64; + let dk: f64; + let y: f64; + let mut hi: f64; + let lo: f64; + let mut val_hi: f64; + let mut val_lo: f64; + let mut hx: u32; + let mut k: i32; + + hx = (ui >> 32) as u32; + k = 0; + if hx < 0x00100000 || (hx >> 31) > 0 { + 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); + + f = x - 1.0; + 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; + + /* See log2.c for details. */ + /* hi+lo = f - hfsq + s*(hfsq+R) ~ log(1+f) */ + hi = f - hfsq; + ui = hi.to_bits(); + ui &= (-1i64 as u64) << 32; + hi = f64::from_bits(ui); + lo = f - hi - hfsq + s * (hfsq + r); + + /* val_hi+val_lo ~ log10(1+f) + k*log10(2) */ + val_hi = hi * IVLN10HI; + dk = k as f64; + y = dk * LOG10_2HI; + val_lo = dk * LOG10_2LO + (lo + hi) * IVLN10LO + lo * IVLN10HI; + + /* + * Extra precision in for adding y is not strictly needed + * since there is no very large cancellation near x = sqrt(2) or + * x = 1/sqrt(2), but we do it anyway since it costs little on CPUs + * with some parallelism and it reduces the error for many args. + */ + w = y + val_hi; + val_lo += (y - w) + val_hi; + val_hi = w; + + return val_lo + val_hi; +} diff --git a/library/compiler-builtins/libm/src/math/log10f.rs b/library/compiler-builtins/libm/src/math/log10f.rs new file mode 100644 index 000000000000..58db093441f3 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/log10f.rs @@ -0,0 +1,76 @@ +use core::f32; + +const IVLN10HI: f32 = 4.3432617188e-01; /* 0x3ede6000 */ +const IVLN10LO: f32 = -3.1689971365e-05; /* 0xb804ead9 */ +const LOG10_2HI: f32 = 3.0102920532e-01; /* 0x3e9a2080 */ +const LOG10_2LO: f32 = 7.9034151668e-07; /* 0x355427db */ +/* |(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 log10f(mut x: f32) -> f32 { + let x1p25f = f32::from_bits(0x4c000000); // 0x1p25f === 2 ^ 25 + + let mut ui: u32 = x.to_bits(); + let hfsq: f32; + let f: f32; + let s: f32; + let z: f32; + let r: f32; + let w: f32; + let t1: f32; + let t2: f32; + let dk: f32; + let mut hi: f32; + let lo: f32; + let mut ix: u32; + let mut k: i32; + + ix = ui; + k = 0; + 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.0; /* log(-#) = NaN */ + } + /* subnormal number, scale up x */ + k -= 25; + x *= x1p25f; + ui = x.to_bits(); + ix = ui; + } 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; + ui = ix; + x = f32::from_bits(ui); + + f = x - 1.0; + 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; + + hi = f - hfsq; + ui = hi.to_bits(); + ui &= 0xfffff000; + hi = f32::from_bits(ui); + lo = f - hi - hfsq + s * (hfsq + r); + dk = k as f32; + return dk * LOG10_2LO + (lo + hi) * IVLN10LO + lo * IVLN10HI + hi * IVLN10HI + dk * LOG10_2HI; +} diff --git a/library/compiler-builtins/libm/src/math/log2.rs b/library/compiler-builtins/libm/src/math/log2.rs new file mode 100644 index 000000000000..c0d3263e3b46 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/log2.rs @@ -0,0 +1,86 @@ +use core::f64; + +const IVLN2HI: f64 = 1.44269504072144627571e+00; /* 0x3ff71547, 0x65200000 */ +const IVLN2LO: f64 = 1.67517131648865118353e-10; /* 0x3de705fc, 0x2eefa200 */ +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 log2(mut x: f64) -> f64 { + let x1p54 = f64::from_bits(0x4350000000000000); // 0x1p54 === 2 ^ 54 + + let mut ui: u64 = x.to_bits(); + let hfsq: f64; + let f: f64; + let s: f64; + let z: f64; + let r: f64; + let mut w: f64; + let t1: f64; + let t2: f64; + let y: f64; + let mut hi: f64; + let lo: f64; + let mut val_hi: f64; + let mut val_lo: f64; + let mut hx: u32; + let mut k: i32; + + hx = (ui >> 32) as u32; + k = 0; + if hx < 0x00100000 || (hx >> 31) > 0 { + 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); + + f = x - 1.0; + 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; + + /* hi+lo = f - hfsq + s*(hfsq+R) ~ log(1+f) */ + hi = f - hfsq; + ui = hi.to_bits(); + ui &= (-1i64 as u64) << 32; + hi = f64::from_bits(ui); + lo = f - hi - hfsq + s * (hfsq + r); + + val_hi = hi * IVLN2HI; + val_lo = (lo + hi) * IVLN2LO + lo * IVLN2HI; + + /* spadd(val_hi, val_lo, y), except for not using double_t: */ + y = k.into(); + w = y + val_hi; + val_lo += (y - w) + val_hi; + val_hi = w; + + return val_lo + val_hi; +} diff --git a/library/compiler-builtins/libm/src/math/log2f.rs b/library/compiler-builtins/libm/src/math/log2f.rs new file mode 100644 index 000000000000..47a917bdb8d2 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/log2f.rs @@ -0,0 +1,71 @@ +use core::f32; + +const IVLN2HI: f32 = 1.4428710938e+00; /* 0x3fb8b000 */ +const IVLN2LO: f32 = -1.7605285393e-04; /* 0xb9389ad4 */ +/* |(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 log2f(mut x: f32) -> f32 { + let x1p25f = f32::from_bits(0x4c000000); // 0x1p25f === 2 ^ 25 + + let mut ui: u32 = x.to_bits(); + let hfsq: f32; + let f: f32; + let s: f32; + let z: f32; + let r: f32; + let w: f32; + let t1: f32; + let t2: f32; + let mut hi: f32; + let lo: f32; + let mut ix: u32; + let mut k: i32; + + ix = ui; + k = 0; + 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.0; /* log(-#) = NaN */ + } + /* subnormal number, scale up x */ + k -= 25; + x *= x1p25f; + ui = x.to_bits(); + ix = ui; + } 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; + ui = ix; + x = f32::from_bits(ui); + + f = x - 1.0; + 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; + + hi = f - hfsq; + ui = hi.to_bits(); + ui &= 0xfffff000; + hi = f32::from_bits(ui); + lo = f - hi - hfsq + s * (hfsq + r); + return (lo + hi) * IVLN2LO + lo * IVLN2HI + hi * IVLN2HI + k as f32; +} diff --git a/library/compiler-builtins/libm/src/math/logf.rs b/library/compiler-builtins/libm/src/math/logf.rs index 76b4ede19514..78c5e94ad92a 100644 --- a/library/compiler-builtins/libm/src/math/logf.rs +++ b/library/compiler-builtins/libm/src/math/logf.rs @@ -1,24 +1,25 @@ -const LN2_HI : f32 = 6.9313812256e-01; /* 0x3f317180 */ -const LN2_LO : f32 = 9.0580006145e-06; /* 0x3717f7d1 */ +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 */ +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 < 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 */ + if (ix >> 31) != 0 { + return (x - x) / 0.; /* log(-#) = NaN */ } /* subnormal number, scale up x */ k -= 25; @@ -32,18 +33,18 @@ pub fn logf(mut x: f32) -> f32 { /* reduce x into [sqrt(2)/2, sqrt(2)] */ ix += 0x3f800000 - 0x3f3504f3; - k += ((ix>>23) as i32) - 0x7f; + k += ((ix >> 23) as i32) - 0x7f; ix = (ix & 0x007fffff) + 0x3f3504f3; - x = f32::from_bits(ix); + 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 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 hfsq = 0.5 * f * f; let dk = k as f32; - s*(hfsq+r) + dk*LN2_LO - hfsq + f + dk*LN2_HI + 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 d844ee153d64..30e7a2c592db 100644 --- a/library/compiler-builtins/libm/src/math/mod.rs +++ b/library/compiler-builtins/libm/src/math/mod.rs @@ -1,16 +1,31 @@ macro_rules! force_eval { ($e:expr) => { - unsafe { ::core::ptr::read_volatile(&$e); } - } + unsafe { + ::core::ptr::read_volatile(&$e); + } + }; } +mod ceilf; +mod expf; mod fabs; mod fabsf; +mod floor; +mod floorf; mod fmodf; +mod hypot; +mod hypotf; +mod log10; +mod log10f; +mod log2; +mod log2f; +mod logf; mod powf; mod round; +mod roundf; mod scalbn; mod scalbnf; +mod sqrt; mod sqrtf; mod logf; mod expf; @@ -20,20 +35,10 @@ mod trunc; mod truncf; pub use self::{ - fabs::fabs, - fabsf::fabsf, - fmodf::fmodf, - powf::powf, - round::round, - scalbn::scalbn, - scalbnf::scalbnf, - sqrtf::sqrtf, - logf::logf, - expf::expf, - floor::floor, - cosf::cosf, - trunc::trunc, - truncf::truncf, + ceilf::ceilf, expf::expf, fabs::fabs, fabsf::fabsf, floor::floor, floorf::floorf, 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, cosf::cosf, }; mod k_cosf; diff --git a/library/compiler-builtins/libm/src/math/round.rs b/library/compiler-builtins/libm/src/math/round.rs index 2a9f67c1d0fa..73d431c51839 100644 --- a/library/compiler-builtins/libm/src/math/round.rs +++ b/library/compiler-builtins/libm/src/math/round.rs @@ -26,7 +26,7 @@ pub fn round(mut x: f64) -> f64 { } else { y = y + x; } - + if i >> 63 != 0 { -y } else { diff --git a/library/compiler-builtins/libm/src/math/roundf.rs b/library/compiler-builtins/libm/src/math/roundf.rs new file mode 100644 index 000000000000..bd2488fa99c8 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/roundf.rs @@ -0,0 +1,33 @@ +use core::f32; + +const TOINT: f32 = 1.0 / f32::EPSILON; + +pub fn roundf(mut x: f32) -> f32 { + let i = x.to_bits(); + let e: u32 = i >> 23 & 0xff; + let mut y: f32; + + if e >= 0x7f + 23 { + return x; + } + if i >> 31 != 0 { + x = -x; + } + if e < 0x7f - 1 { + force_eval!(x + TOINT); + return 0.0 * x; + } + y = x + TOINT - TOINT - x; + if y > 0.5f32 { + y = y + x - 1.0; + } else if y <= -0.5f32 { + y = y + x + 1.0; + } else { + y = y + x; + } + if i >> 31 != 0 { + -y + } else { + y + } +} diff --git a/library/compiler-builtins/libm/src/math/scalbn.rs b/library/compiler-builtins/libm/src/math/scalbn.rs index 76e06f03ece2..ad81072ddc98 100644 --- a/library/compiler-builtins/libm/src/math/scalbn.rs +++ b/library/compiler-builtins/libm/src/math/scalbn.rs @@ -1,11 +1,11 @@ #[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 +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; @@ -29,5 +29,5 @@ pub fn scalbn(x : f64, mut n: i32) -> f64 { } } } - y*f64::from_bits(((0x3ff+n) as u64)<<52) + 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 31f93d323c5c..901497e5e51d 100644 --- a/library/compiler-builtins/libm/src/math/scalbnf.rs +++ b/library/compiler-builtins/libm/src/math/scalbnf.rs @@ -1,9 +1,9 @@ #[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 - +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 { x *= x1p127; n -= 127; @@ -25,5 +25,5 @@ pub fn scalbnf(mut x: f32, mut n : i32) -> f32 { } } } - x * f32::from_bits(((0x7f+n) as u32)<<23) + 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..17de5a2e0308 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/sqrt.rs @@ -0,0 +1,131 @@ +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 index b50ffd77156a..6bea67cbc165 100644 --- a/library/compiler-builtins/libm/src/math/trunc.rs +++ b/library/compiler-builtins/libm/src/math/trunc.rs @@ -22,3 +22,11 @@ pub fn trunc(x: f64) -> f64 { 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 index f7d7249e3519..9d42620d9666 100644 --- a/library/compiler-builtins/libm/src/math/truncf.rs +++ b/library/compiler-builtins/libm/src/math/truncf.rs @@ -22,3 +22,11 @@ pub fn truncf(x: f32) -> f32 { 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 40eac207c153..65b03f585e23 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, ); }} @@ -561,21 +652,21 @@ fn main() -> Result<(), Box> { // With signature `fn(f32) -> f32` f32_f32! { // acosf, - // floorf, + floorf, truncf, // asinf, // atanf, // cbrtf, - // ceilf, cosf, + ceilf, // coshf, // exp2f, expf, // fdimf, - // log10f, - // log2f, + log10f, + log2f, logf, - // roundf, + roundf, // sinf, // sinhf, // tanf, @@ -587,7 +678,7 @@ f32_f32! { // With signature `fn(f32, f32) -> f32` f32f32_f32! { // atan2f, - // hypotf, + hypotf, fmodf, powf, } @@ -616,13 +707,13 @@ f64_f64! { // expm1, floor, // log, - // log10, + log10, // log1p, - // log2, + log2, round, // sin, // sinh, - // sqrt, + sqrt, // tan, // tanh, trunc, @@ -634,7 +725,7 @@ f64f64_f64! { // atan2, // fdim, // fmod, - // hypot, + hypot, // pow, }