From 7c4aef41b11a4928a55e2506b3a9eb4aa1fbff91 Mon Sep 17 00:00:00 2001 From: Trevor Gross Date: Tue, 22 Apr 2025 22:53:17 -0400 Subject: [PATCH 1/5] Migrate to a generic rem_pio2 --- libm-test/src/f8_impl.rs | 3 + libm/src/math/generic/mod.rs | 2 + libm/src/math/generic/rem_frac_pi_2.rs | 237 +++++++++++++++++++++++++ libm/src/math/rem_pio2.rs | 146 +-------------- libm/src/math/support/float_traits.rs | 11 +- libm/src/math/support/int_traits.rs | 46 +++-- libm/src/math/support/mod.rs | 2 +- 7 files changed, 284 insertions(+), 163 deletions(-) create mode 100644 libm/src/math/generic/rem_frac_pi_2.rs diff --git a/libm-test/src/f8_impl.rs b/libm-test/src/f8_impl.rs index 905c7d7fd..898b9ca60 100644 --- a/libm-test/src/f8_impl.rs +++ b/libm-test/src/f8_impl.rs @@ -24,6 +24,9 @@ impl Float for f8 { const ZERO: Self = Self(0b0_0000_000); const NEG_ZERO: Self = Self(0b1_0000_000); const ONE: Self = Self(0b0_0111_000); + const TWO: Self = Self(0b0_1000_000); + const THREE: Self = Self(0b0_1000_100); + const FOUR: Self = Self(0b0_1001_000); const NEG_ONE: Self = Self(0b1_0111_000); const MAX: Self = Self(0b0_1110_111); const MIN: Self = Self(0b1_1110_111); diff --git a/libm/src/math/generic/mod.rs b/libm/src/math/generic/mod.rs index 35846351a..11ab645f1 100644 --- a/libm/src/math/generic/mod.rs +++ b/libm/src/math/generic/mod.rs @@ -13,6 +13,7 @@ mod fmin; mod fminimum; mod fminimum_num; mod fmod; +mod rem_frac_pi_2; mod rint; mod round; mod scalbn; @@ -31,6 +32,7 @@ pub use fmin::fmin; pub use fminimum::fminimum; pub use fminimum_num::fminimum_num; pub use fmod::fmod; +pub(crate) use rem_frac_pi_2::{medium, rem_frac_pi_2}; pub use rint::rint_round; pub use round::round; pub use scalbn::scalbn; diff --git a/libm/src/math/generic/rem_frac_pi_2.rs b/libm/src/math/generic/rem_frac_pi_2.rs new file mode 100644 index 000000000..6038a97b0 --- /dev/null +++ b/libm/src/math/generic/rem_frac_pi_2.rs @@ -0,0 +1,237 @@ +#![allow(unused)] + +use core::f64::consts; + +use crate::support::{CastFrom, CastInto, DInt, Float, HInt, HalfRep, Int, MinInt}; + +pub(crate) trait RemFracPi2Support: Float> { + const TO_INT: Self; + const INV_PIO2: Self; + const PIO2_1: Self; + const PIO2_1T: Self; + const PIO2_2: Self; + const PIO2_2T: Self; + const PIO2_3: Self; + const PIO2_3T: Self; + + const FRAC_5PI_4_HI: HalfRep; + const FRAC_3PI_4_HI: HalfRep; + const FRAC_9PI_4_HI: HalfRep; + const FRAC_7PI_4_HI: HalfRep; + const FRAC_3PI_2_HI: HalfRep; + /// 2pi + const TAU_HI: HalfRep; + const FRAC_PI_2_HI: HalfRep; + /// (2^20)(pi/2) + const FRAC_2_POW_20_PI_2: HalfRep; + + const MAGIC: u32 = 23; + const MAGIC_F: Self; + + fn large(x: &[Self], y: &mut [Self], e0: i32, prec: usize) -> i32; +} + +const EPS: f64 = 2.2204460492503131e-16; + +impl RemFracPi2Support for f64 { + const TO_INT: Self = 1.5 / EPS; + const INV_PIO2: Self = 6.36619772367581382433e-01; + const PIO2_1: Self = 1.57079632673412561417e+00; + const PIO2_1T: Self = 6.07710050650619224932e-11; + const PIO2_2: Self = 6.07710050630396597660e-11; + const PIO2_2T: Self = 2.02226624879595063154e-21; + const PIO2_3: Self = 2.02226624871116645580e-21; + const PIO2_3T: Self = 8.47842766036889956997e-32; + const FRAC_5PI_4_HI: HalfRep = 0x400f6a7a; + const FRAC_3PI_4_HI: HalfRep = 0x4002d97c; + const FRAC_9PI_4_HI: HalfRep = 0x401c463b; + const FRAC_7PI_4_HI: HalfRep = 0x4015fdbc; + const FRAC_3PI_2_HI: HalfRep = 0x4012d97c; + const TAU_HI: HalfRep = 0x401921fb; + const FRAC_PI_2_HI: HalfRep = 0x921fb; + const FRAC_2_POW_20_PI_2: HalfRep = 0x413921fb; + + const MAGIC_F: Self = hf64!("0x1p24"); + + fn large(x: &[Self], y: &mut [Self], e0: i32, prec: usize) -> i32 { + super::super::rem_pio2_large(x, y, e0, prec) + } +} + +pub(crate) fn rem_frac_pi_2(x: F) -> (i32, F, F) +where + F: RemFracPi2Support, + F: CastInto, + HalfRep: Int + MinInt>, + // as Int>::Unsigned: Int, +{ + // let sign = x.is_sign_positive() + + let ix: HalfRep = x.abs().to_bits().hi(); + let pos = x.is_sign_positive(); + + if ix <= F::FRAC_5PI_4_HI { + /* |x| ~<= 5pi/4 */ + if (ix & F::SIG_MASK.hi()) == F::FRAC_PI_2_HI { + /* |x| ~= pi/2 or 2pi/2 */ + return medium(x, ix); /* cancellation -- use medium case */ + } + + if ix <= F::FRAC_3PI_4_HI { + /* |x| ~<= 3pi/4 */ + if pos { + let z = x - F::PIO2_1; /* one round good to 85 bits */ + let y0 = z - F::PIO2_1T; + let y1 = (z - y0) - F::PIO2_1T; + return (1, y0, y1); + } else { + let z = x + F::PIO2_1; + let y0 = z + F::PIO2_1T; + let y1 = (z - y0) + F::PIO2_1T; + return (-1, y0, y1); + } + } else if pos { + let z = x - F::TWO * F::PIO2_1; + let y0 = z - F::TWO * F::PIO2_1T; + let y1 = (z - y0) - F::TWO * F::PIO2_1T; + return (2, y0, y1); + } else { + let z = x + F::TWO * F::PIO2_1; + let y0 = z + F::TWO * F::PIO2_1T; + let y1 = (z - y0) + F::TWO * F::PIO2_1T; + return (-2, y0, y1); + } + } + + if ix <= F::FRAC_9PI_4_HI { + /* |x| ~<= 9pi/4 */ + if ix <= F::FRAC_7PI_4_HI { + /* |x| ~<= 7pi/4 */ + if ix == F::FRAC_3PI_2_HI { + /* |x| ~= 3pi/2 */ + return medium(x, ix); + } + + if pos { + let z = x - F::THREE * F::PIO2_1; + let y0 = z - F::THREE * F::PIO2_1T; + let y1 = (z - y0) - F::THREE * F::PIO2_1T; + return (3, y0, y1); + } else { + let z = x + F::THREE * F::PIO2_1; + let y0 = z + F::THREE * F::PIO2_1T; + let y1 = (z - y0) + F::THREE * F::PIO2_1T; + return (-3, y0, y1); + } + } else { + if ix == F::TAU_HI { + /* |x| ~= 4pi/2 */ + return medium(x, ix); + } + + if pos { + let z = x - F::FOUR * F::PIO2_1; + let y0 = z - F::FOUR * F::PIO2_1T; + let y1 = (z - y0) - F::FOUR * F::PIO2_1T; + return (4, y0, y1); + } else { + let z = x + F::FOUR * F::PIO2_1; + let y0 = z + F::FOUR * F::PIO2_1T; + let y1 = (z - y0) + F::FOUR * F::PIO2_1T; + return (-4, y0, y1); + } + } + } + + if ix < F::FRAC_2_POW_20_PI_2 { + /* |x| ~< 2^20*(pi/2), medium size */ + return medium(x, ix); + } + /* + + * all other (large) arguments + */ + if ix >= F::EXP_MASK.hi() { + /* x is inf or NaN */ + let y0 = x - x; + let y1 = y0; + return (0, y0, y1); + } + + /* set z = scalbn(|x|,-ilogb(x)+23) */ + let mut ui = x.to_bits(); + ui &= F::SIG_MASK; + // ui |= F::Int::cast_from((F::EXP_BIAS + F::MAGIC) << F::SIG_BITS); + ui |= F::Int::cast_from(F::EXP_BIAS + F::MAGIC) << F::SIG_BITS; + + let mut z = F::from_bits(ui); + let mut tx = [F::ZERO; 3]; + + for i in 0..2 { + // ?? + i!(tx, i, =, super::trunc(z)); + z = (z - i!(tx, i)) * F::MAGIC_F; + } + + i!(tx,2, =, z); + + /* skip zero terms, first term is non-zero */ + let mut i = 2; + while i != 0 && i!(tx, i) == F::ZERO { + i -= 1; + } + + let mut ty = [F::ZERO; 3]; + + let ex: i32 = (ix >> (HalfRep::::BITS - F::EXP_BITS - 1)).cast(); + let n = F::large(&tx[..=i], &mut ty, ex - (F::EXP_BIAS + F::MAGIC) as i32, 1); + + if !pos { + return (-n, -i!(ty, 0), -i!(ty, 1)); + } else { + (n, i!(ty, 0), i!(ty, 1)) + } +} + +pub fn medium(x: F, ix: HalfRep) -> (i32, F, F) +where + F: RemFracPi2Support, + F: CastInto, + HalfRep: Int, +{ + /* rint(x/(pi/2)), Assume round-to-nearest. */ + let tmp = x * F::INV_PIO2 + F::TO_INT; + // force rounding of tmp to its storage format on x87 to avoid + // excess precision issues. + #[cfg(all(target_arch = "x86", not(target_feature = "sse2")))] + let tmp = force_eval!(tmp); + let f_n = tmp - F::TO_INT; + let n: i32 = f_n.cast(); + let mut r = x - f_n * F::PIO2_1; + let mut w = f_n * F::PIO2_1T; /* 1st round, good to 85 bits */ + let mut y0 = r - w; + let ui = y0.to_bits(); + let ey = y0.ex().signed(); + let ex: i32 = (ix >> (HalfRep::::BITS - F::EXP_BITS - 1)).cast(); + + // (ix >> 20) as i32; + if ex - ey > 16 { + /* 2nd round, good to 118 bits */ + let t = r; + w = f_n * F::PIO2_2; + r = t - w; + w = f_n * F::PIO2_2T - ((t - r) - w); + y0 = r - w; + let ey = y0.ex().signed(); + if ex - ey > 49 { + /* 3rd round, good to 151 bits, covers all cases */ + let t = r; + w = f_n * F::PIO2_3; + r = t - w; + w = f_n * F::PIO2_3T - ((t - r) - w); + y0 = r - w; + } + } + let y1 = (r - y0) - w; + (n, y0, y1) +} diff --git a/libm/src/math/rem_pio2.rs b/libm/src/math/rem_pio2.rs index d677fd9dc..d5708c68b 100644 --- a/libm/src/math/rem_pio2.rs +++ b/libm/src/math/rem_pio2.rs @@ -1,3 +1,4 @@ +#![allow(unused)] // origin: FreeBSD /usr/src/lib/msun/src/e_rem_pio2.c // // ==================================================== @@ -43,150 +44,7 @@ const PIO2_3T: f64 = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ // caller must handle the case when reduction is not needed: |x| ~<= pi/4 */ #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub(crate) 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; - - fn medium(x: f64, ix: u32) -> (i32, f64, f64) { - /* rint(x/(pi/2)), Assume round-to-nearest. */ - let tmp = x * INV_PIO2 + TO_INT; - // force rounding of tmp to it's storage format on x87 to avoid - // excess precision issues. - #[cfg(all(target_arch = "x86", not(target_feature = "sse2")))] - let tmp = force_eval!(tmp); - let f_n = tmp - TO_INT; - let n = f_n as i32; - let mut r = x - f_n * PIO2_1; - let mut w = f_n * PIO2_1T; /* 1st round, good to 85 bits */ - let mut y0 = r - w; - let ui = f64::to_bits(y0); - let ey = (ui >> 52) as i32 & 0x7ff; - let ex = (ix >> 20) as i32; - if ex - ey > 16 { - /* 2nd round, good to 118 bits */ - let t = r; - w = f_n * PIO2_2; - r = t - w; - w = f_n * PIO2_2T - ((t - r) - w); - y0 = r - w; - let ey = (f64::to_bits(y0) >> 52) as i32 & 0x7ff; - if ex - ey > 49 { - /* 3rd round, good to 151 bits, covers all cases */ - let t = r; - w = f_n * PIO2_3; - r = t - w; - w = f_n * PIO2_3T - ((t - r) - w); - y0 = r - w; - } - } - let y1 = (r - y0) - w; - (n, y0, y1) - } - - if ix <= 0x400f6a7a { - /* |x| ~<= 5pi/4 */ - if (ix & 0xfffff) == 0x921fb { - /* |x| ~= pi/2 or 2pi/2 */ - return medium(x, ix); /* cancellation -- use medium case */ - } - if ix <= 0x4002d97c { - /* |x| ~<= 3pi/4 */ - if sign == 0 { - let z = x - PIO2_1; /* one round good to 85 bits */ - let y0 = z - PIO2_1T; - let y1 = (z - y0) - PIO2_1T; - return (1, y0, y1); - } else { - let z = x + PIO2_1; - let y0 = z + PIO2_1T; - let y1 = (z - y0) + PIO2_1T; - return (-1, y0, y1); - } - } else if sign == 0 { - let z = x - 2.0 * PIO2_1; - let y0 = z - 2.0 * PIO2_1T; - let y1 = (z - y0) - 2.0 * PIO2_1T; - return (2, y0, y1); - } else { - let z = x + 2.0 * PIO2_1; - let y0 = z + 2.0 * PIO2_1T; - let y1 = (z - y0) + 2.0 * PIO2_1T; - return (-2, y0, y1); - } - } - if ix <= 0x401c463b { - /* |x| ~<= 9pi/4 */ - if ix <= 0x4015fdbc { - /* |x| ~<= 7pi/4 */ - if ix == 0x4012d97c { - /* |x| ~= 3pi/2 */ - return medium(x, ix); - } - if sign == 0 { - let z = x - 3.0 * PIO2_1; - let y0 = z - 3.0 * PIO2_1T; - let y1 = (z - y0) - 3.0 * PIO2_1T; - return (3, y0, y1); - } else { - let z = x + 3.0 * PIO2_1; - let y0 = z + 3.0 * PIO2_1T; - let y1 = (z - y0) + 3.0 * PIO2_1T; - return (-3, y0, y1); - } - } else { - if ix == 0x401921fb { - /* |x| ~= 4pi/2 */ - return medium(x, ix); - } - if sign == 0 { - let z = x - 4.0 * PIO2_1; - let y0 = z - 4.0 * PIO2_1T; - let y1 = (z - y0) - 4.0 * PIO2_1T; - return (4, y0, y1); - } else { - let z = x + 4.0 * PIO2_1; - let y0 = z + 4.0 * PIO2_1T; - let y1 = (z - y0) + 4.0 * PIO2_1T; - return (-4, y0, y1); - } - } - } - if ix < 0x413921fb { - /* |x| ~< 2^20*(pi/2), medium size */ - return medium(x, ix); - } - /* - * all other (large) arguments - */ - if ix >= 0x7ff00000 { - /* x is inf or NaN */ - let y0 = x - x; - let y1 = y0; - return (0, y0, y1); - } - /* set z = scalbn(|x|,-ilogb(x)+23) */ - let mut ui = f64::to_bits(x); - ui &= (!1) >> 12; - ui |= (0x3ff + 23) << 52; - let mut z = f64::from_bits(ui); - let mut tx = [0.0; 3]; - for i in 0..2 { - i!(tx,i, =, z as i32 as f64); - z = (z - i!(tx, i)) * x1p24; - } - i!(tx,2, =, z); - /* skip zero terms, first term is non-zero */ - let mut i = 2; - while i != 0 && i!(tx, i) == 0.0 { - i -= 1; - } - let mut ty = [0.0; 3]; - let n = rem_pio2_large(&tx[..=i], &mut ty, ((ix as i32) >> 20) - (0x3ff + 23), 1); - if sign != 0 { - return (-n, -i!(ty, 0), -i!(ty, 1)); - } - (n, i!(ty, 0), i!(ty, 1)) + return super::generic::rem_frac_pi_2(x); } #[cfg(test)] diff --git a/libm/src/math/support/float_traits.rs b/libm/src/math/support/float_traits.rs index 8094a7b84..663e40117 100644 --- a/libm/src/math/support/float_traits.rs +++ b/libm/src/math/support/float_traits.rs @@ -1,6 +1,9 @@ use core::{fmt, mem, ops}; -use super::int_traits::{CastFrom, Int, MinInt}; +use super::int_traits::{CastFrom, DInt, Int, MinInt}; + +/// Wrapper to extract the integer type half of the float's size +pub type HalfRep = <::Int as DInt>::H; /// Trait for some basic operations on floats // #[allow(dead_code)] @@ -38,6 +41,9 @@ pub trait Float: const MAX: Self; const MIN: Self; const EPSILON: Self; + const TWO: Self; + const THREE: Self; + const FOUR: Self; const PI: Self; const NEG_PI: Self; const FRAC_PI_2: Self; @@ -222,6 +228,9 @@ macro_rules! float_impl { // Sign bit set, saturated mantissa, saturated exponent with last bit zeroed const MIN: Self = $from_bits(Self::Int::MAX & !(1 << Self::SIG_BITS)); const EPSILON: Self = <$ty>::EPSILON; + const TWO: Self = 2.0; + const THREE: Self = 3.0; + const FOUR: Self = 4.0; // Exponent is a 1 in the LSB const MIN_POSITIVE_NORMAL: Self = $from_bits(1 << Self::SIG_BITS); diff --git a/libm/src/math/support/int_traits.rs b/libm/src/math/support/int_traits.rs index 3ec1faba1..0e10530ff 100644 --- a/libm/src/math/support/int_traits.rs +++ b/libm/src/math/support/int_traits.rs @@ -410,26 +410,38 @@ macro_rules! cast_into { )*}; } -macro_rules! cast_into_float { - ($ty:ty) => { +macro_rules! cast_int_float { + ($ity:ty) => { #[cfg(f16_enabled)] - cast_into_float!($ty; f16); + cast_int_float!($ity; f16); - cast_into_float!($ty; f32, f64); + cast_int_float!($ity; f32, f64); #[cfg(f128_enabled)] - cast_into_float!($ty; f128); + cast_int_float!($ity; f128); }; - ($ty:ty; $($into:ty),*) => {$( - impl CastInto<$into> for $ty { - fn cast(self) -> $into { + ($ity:ty; $($fty:ty),*) => {$( + impl CastInto<$ity> for $fty { + fn cast(self) -> $ity { #[cfg(not(feature = "compiler-builtins"))] - debug_assert_eq!(self as $into as $ty, self, "inexact float cast"); - self as $into + debug_assert_eq!(self as $ity as $fty, self, "inexact float->int cast"); + self as $ity } - fn cast_lossy(self) -> $into { - self as $into + fn cast_lossy(self) -> $ity { + self as $ity + } + } + + impl CastInto<$fty> for $ity { + fn cast(self) -> $fty { + #[cfg(not(feature = "compiler-builtins"))] + debug_assert_eq!(self as $fty as $ity, self, "inexact int->float cast"); + self as $fty + } + + fn cast_lossy(self) -> $fty { + self as $fty } } )*}; @@ -448,8 +460,8 @@ cast_into!(i64); cast_into!(u128); cast_into!(i128); -cast_into_float!(i8); -cast_into_float!(i16); -cast_into_float!(i32); -cast_into_float!(i64); -cast_into_float!(i128); +cast_int_float!(i8); +cast_int_float!(i16); +cast_int_float!(i32); +cast_int_float!(i64); +cast_int_float!(i128); diff --git a/libm/src/math/support/mod.rs b/libm/src/math/support/mod.rs index ee3f2bbdf..1f0b59720 100644 --- a/libm/src/math/support/mod.rs +++ b/libm/src/math/support/mod.rs @@ -10,7 +10,7 @@ mod int_traits; pub use big::{i256, u256}; pub use env::{FpResult, Round, Status}; #[allow(unused_imports)] -pub use float_traits::{DFloat, Float, HFloat, IntTy}; +pub use float_traits::{DFloat, Float, HFloat, HalfRep, IntTy}; pub(crate) use float_traits::{f32_from_bits, f64_from_bits}; #[cfg(f16_enabled)] #[allow(unused_imports)] From f356c0aac1e9b3b8d2b06630569f815aa7a81bae Mon Sep 17 00:00:00 2001 From: Trevor Gross Date: Wed, 23 Apr 2025 00:05:45 -0400 Subject: [PATCH 2/5] Generate constants --- etc/generate-constants.jl | 120 ++++++++++++++++++ libm/src/math/generic/generated/mod.rs | 3 + libm/src/math/generic/generated/rem_pio2.rs | 29 +++++ libm/src/math/generic/mod.rs | 5 +- .../generic/{rem_frac_pi_2.rs => rem_pio2.rs} | 51 ++------ libm/src/math/rem_pio2.rs | 27 +--- libm/src/math/support/mod.rs | 2 +- 7 files changed, 170 insertions(+), 67 deletions(-) create mode 100644 etc/generate-constants.jl create mode 100644 libm/src/math/generic/generated/mod.rs create mode 100644 libm/src/math/generic/generated/rem_pio2.rs rename libm/src/math/generic/{rem_frac_pi_2.rs => rem_pio2.rs} (79%) diff --git a/etc/generate-constants.jl b/etc/generate-constants.jl new file mode 100644 index 000000000..30876740b --- /dev/null +++ b/etc/generate-constants.jl @@ -0,0 +1,120 @@ +using Printf + +function main() + repo_root = dirname(@__DIR__) + out_dir = joinpath(repo_root, "libm/src/math/generic/generated") + + mod_contents = "// autogenerated\n\n" + + cfg = Config(out_dir) + mod_contents *= update_rem_pio2(cfg) + + path = joinpath(cfg.out_dir, "mod.rs") + write(path, mod_contents) + println("wrote $path") +end + +struct Config + out_dir::String +end + +@enum FloatTy F16 F32 F64 F128 + +struct TyInfo + bits::UInt32 + sig_bits::UInt32 + ity::Type + ty_name::String +end + +function ty_info(fty::FloatTy)::TyInfo + if fty == F16 + bits = 16 + sig_bits = 10 + ity = UInt16 + elseif fty == F32 + bits = 32 + sig_bits = 23 + ity = UInt32 + elseif fty == F64 + bits = 64 + sig_bits = 52 + ity = UInt64 + elseif fty == F128 + bits = 128 + sig_bits = 112 + ity = UInt128 + else + @assert(false) + end + + return TyInfo( + bits, + sig_bits, + ity, + lowercase(string(fty)), + ) +end + +function update_rem_pio2(cfg::Config)::String + prefix= """ + use crate::support::HalfRep; + use super::super::rem_pio2::RemPio2Support; + """ + ret = "" + + for fty in [F64] + info = ty_info(fty) + ty_name = info.ty_name + halfbits = info.bits / 2 + + if info.bits == 32 || info.bits == 64 + to_bits = "$(ty_name)_to_bits" + prefix *= "use crate::support::$to_bits;\n" + else + to_bits = "$ty_name::to_bits" + end + + hi_int(x::BigFloat) = @sprintf "(%s(hf%d!(\"%a\")) >> %d) as u%d" to_bits info.bits x halfbits halfbits + + setprecision(ty_info(fty).sig_bits + 1) + + ty_impl = """ + impl RemPio2Support for $ty_name { + const TO_INT: Self = 1.5 / $ty_name::EPSILON; + const INV_PIO2: Self = 6.36619772367581382433e-01; + const PIO2_1: Self = 1.57079632673412561417e+00; + const PIO2_1T: Self = 6.07710050650619224932e-11; + const PIO2_2: Self = 6.07710050630396597660e-11; + const PIO2_2T: Self = 2.02226624879595063154e-21; + const PIO2_3: Self = 2.02226624871116645580e-21; + const PIO2_3T: Self = 8.47842766036889956997e-32; + + const FRAC_5PI_4_HI: HalfRep = $(hi_int(big(pi)*5/4)); + const FRAC_3PI_4_HI: HalfRep = $(hi_int(big(pi)*3/4)); + const FRAC_9PI_4_HI: HalfRep = $(hi_int(big(pi)*9/4)); + const FRAC_7PI_4_HI: HalfRep = $(hi_int(big(pi)*7/4)); + const FRAC_3PI_2_HI: HalfRep = $(hi_int(big(pi)*3/2)); + const TAU_HI: HalfRep = $(hi_int(big(pi)*2)); + const FRAC_PI_2_HI: HalfRep = 0x921fb; + const FRAC_2_POW_20_PI_2: HalfRep = $(hi_int((big(2)^20) * pi / 2)); + + const MAGIC_F: Self = hf64!("0x1p24"); + + fn large(x: &[Self], y: &mut [Self], e0: i32, prec: usize) -> i32 { + super::super::super::rem_pio2_large(x, y, e0, prec) + } + } + """ + ret *= ty_impl + end + + ret = prefix * "\n" * ret + path = joinpath(cfg.out_dir, "rem_pio2.rs") + write(path, ret) + println("wrote $path") + + return "mod rem_pio2;" +end + +main() diff --git a/libm/src/math/generic/generated/mod.rs b/libm/src/math/generic/generated/mod.rs new file mode 100644 index 000000000..4f36e0561 --- /dev/null +++ b/libm/src/math/generic/generated/mod.rs @@ -0,0 +1,3 @@ +// autogenerated + +mod rem_pio2; \ No newline at end of file diff --git a/libm/src/math/generic/generated/rem_pio2.rs b/libm/src/math/generic/generated/rem_pio2.rs new file mode 100644 index 000000000..035f5c4ef --- /dev/null +++ b/libm/src/math/generic/generated/rem_pio2.rs @@ -0,0 +1,29 @@ +use crate::support::HalfRep; +use super::super::rem_pio2::RemPio2Support; +use crate::support::f64_to_bits; + +impl RemPio2Support for f64 { + const TO_INT: Self = 1.5 / f64::EPSILON; + const INV_PIO2: Self = 6.36619772367581382433e-01; + const PIO2_1: Self = 1.57079632673412561417e+00; + const PIO2_1T: Self = 6.07710050650619224932e-11; + const PIO2_2: Self = 6.07710050630396597660e-11; + const PIO2_2T: Self = 2.02226624879595063154e-21; + const PIO2_3: Self = 2.02226624871116645580e-21; + const PIO2_3T: Self = 8.47842766036889956997e-32; + + const FRAC_5PI_4_HI: HalfRep = (f64_to_bits(hf64!("0x3.ed4f452aa70bcp+0")) >> 32) as u32; + const FRAC_3PI_4_HI: HalfRep = (f64_to_bits(hf64!("0x2.5b2f8fe6643a4p+0")) >> 32) as u32; + const FRAC_9PI_4_HI: HalfRep = (f64_to_bits(hf64!("0x7.118eafb32caecp+0")) >> 32) as u32; + const FRAC_7PI_4_HI: HalfRep = (f64_to_bits(hf64!("0x5.7f6efa6ee9dd4p+0")) >> 32) as u32; + const FRAC_3PI_2_HI: HalfRep = (f64_to_bits(hf64!("0x4.b65f1fccc8748p+0")) >> 32) as u32; + const TAU_HI: HalfRep = (f64_to_bits(hf64!("0x6.487ed5110b46p+0")) >> 32) as u32; + const FRAC_PI_2_HI: HalfRep = 0x921fb; + const FRAC_2_POW_20_PI_2: HalfRep = (f64_to_bits(hf64!("0x1.921fb54442d18p+20")) >> 32) as u32; + + const MAGIC_F: Self = hf64!("0x1p24"); + + fn large(x: &[Self], y: &mut [Self], e0: i32, prec: usize) -> i32 { + super::super::super::rem_pio2_large(x, y, e0, prec) + } +} diff --git a/libm/src/math/generic/mod.rs b/libm/src/math/generic/mod.rs index 11ab645f1..cd75aadf5 100644 --- a/libm/src/math/generic/mod.rs +++ b/libm/src/math/generic/mod.rs @@ -13,7 +13,8 @@ mod fmin; mod fminimum; mod fminimum_num; mod fmod; -mod rem_frac_pi_2; +mod generated; +mod rem_pio2; mod rint; mod round; mod scalbn; @@ -32,7 +33,7 @@ pub use fmin::fmin; pub use fminimum::fminimum; pub use fminimum_num::fminimum_num; pub use fmod::fmod; -pub(crate) use rem_frac_pi_2::{medium, rem_frac_pi_2}; +pub(crate) use rem_pio2::rem_pio2; pub use rint::rint_round; pub use round::round; pub use scalbn::scalbn; diff --git a/libm/src/math/generic/rem_frac_pi_2.rs b/libm/src/math/generic/rem_pio2.rs similarity index 79% rename from libm/src/math/generic/rem_frac_pi_2.rs rename to libm/src/math/generic/rem_pio2.rs index 6038a97b0..7b3d9fa5c 100644 --- a/libm/src/math/generic/rem_frac_pi_2.rs +++ b/libm/src/math/generic/rem_pio2.rs @@ -4,14 +4,21 @@ use core::f64::consts; use crate::support::{CastFrom, CastInto, DInt, Float, HInt, HalfRep, Int, MinInt}; -pub(crate) trait RemFracPi2Support: Float> { +pub(crate) trait RemPio2Support: Float> { const TO_INT: Self; + /// 53 bits of 2/pi const INV_PIO2: Self; + /// first 33 bits of pi/2 const PIO2_1: Self; + /// pi/2 - PIO2_1 const PIO2_1T: Self; + /// second 33 bits of pi/2 const PIO2_2: Self; + /// pi/2 - (PIO2_1+PIO2_2) const PIO2_2T: Self; + /// third 33 bits of pi/2 const PIO2_3: Self; + /// pi/2 - (PIO2_1+PIO2_2+PIO2_3) const PIO2_3T: Self; const FRAC_5PI_4_HI: HalfRep; @@ -31,42 +38,12 @@ pub(crate) trait RemFracPi2Support: Float> { fn large(x: &[Self], y: &mut [Self], e0: i32, prec: usize) -> i32; } -const EPS: f64 = 2.2204460492503131e-16; - -impl RemFracPi2Support for f64 { - const TO_INT: Self = 1.5 / EPS; - const INV_PIO2: Self = 6.36619772367581382433e-01; - const PIO2_1: Self = 1.57079632673412561417e+00; - const PIO2_1T: Self = 6.07710050650619224932e-11; - const PIO2_2: Self = 6.07710050630396597660e-11; - const PIO2_2T: Self = 2.02226624879595063154e-21; - const PIO2_3: Self = 2.02226624871116645580e-21; - const PIO2_3T: Self = 8.47842766036889956997e-32; - const FRAC_5PI_4_HI: HalfRep = 0x400f6a7a; - const FRAC_3PI_4_HI: HalfRep = 0x4002d97c; - const FRAC_9PI_4_HI: HalfRep = 0x401c463b; - const FRAC_7PI_4_HI: HalfRep = 0x4015fdbc; - const FRAC_3PI_2_HI: HalfRep = 0x4012d97c; - const TAU_HI: HalfRep = 0x401921fb; - const FRAC_PI_2_HI: HalfRep = 0x921fb; - const FRAC_2_POW_20_PI_2: HalfRep = 0x413921fb; - - const MAGIC_F: Self = hf64!("0x1p24"); - - fn large(x: &[Self], y: &mut [Self], e0: i32, prec: usize) -> i32 { - super::super::rem_pio2_large(x, y, e0, prec) - } -} - -pub(crate) fn rem_frac_pi_2(x: F) -> (i32, F, F) +pub(crate) fn rem_pio2(x: F) -> (i32, F, F) where - F: RemFracPi2Support, + F: RemPio2Support, F: CastInto, HalfRep: Int + MinInt>, - // as Int>::Unsigned: Int, { - // let sign = x.is_sign_positive() - let ix: HalfRep = x.abs().to_bits().hi(); let pos = x.is_sign_positive(); @@ -148,9 +125,8 @@ where return medium(x, ix); } /* - - * all other (large) arguments - */ + * all other (large) arguments + */ if ix >= F::EXP_MASK.hi() { /* x is inf or NaN */ let y0 = x - x; @@ -161,7 +137,6 @@ where /* set z = scalbn(|x|,-ilogb(x)+23) */ let mut ui = x.to_bits(); ui &= F::SIG_MASK; - // ui |= F::Int::cast_from((F::EXP_BIAS + F::MAGIC) << F::SIG_BITS); ui |= F::Int::cast_from(F::EXP_BIAS + F::MAGIC) << F::SIG_BITS; let mut z = F::from_bits(ui); @@ -195,7 +170,7 @@ where pub fn medium(x: F, ix: HalfRep) -> (i32, F, F) where - F: RemFracPi2Support, + F: RemPio2Support, F: CastInto, HalfRep: Int, { diff --git a/libm/src/math/rem_pio2.rs b/libm/src/math/rem_pio2.rs index d5708c68b..92dd6fb34 100644 --- a/libm/src/math/rem_pio2.rs +++ b/libm/src/math/rem_pio2.rs @@ -13,38 +13,13 @@ // Optimized by Bruce D. Evans. */ use super::rem_pio2_large; -// #if FLT_EVAL_METHOD==0 || FLT_EVAL_METHOD==1 -// #define EPS DBL_EPSILON -const EPS: f64 = 2.2204460492503131e-16; -// #elif FLT_EVAL_METHOD==2 -// #define EPS LDBL_EPSILON -// #endif - -// TODO: Support FLT_EVAL_METHOD? - -const TO_INT: f64 = 1.5 / EPS; -/// 53 bits of 2/pi -const INV_PIO2: f64 = 6.36619772367581382433e-01; /* 0x3FE45F30, 0x6DC9C883 */ -/// first 33 bits of pi/2 -const PIO2_1: f64 = 1.57079632673412561417e+00; /* 0x3FF921FB, 0x54400000 */ -/// pi/2 - PIO2_1 -const PIO2_1T: f64 = 6.07710050650619224932e-11; /* 0x3DD0B461, 0x1A626331 */ -/// second 33 bits of pi/2 -const PIO2_2: f64 = 6.07710050630396597660e-11; /* 0x3DD0B461, 0x1A600000 */ -/// pi/2 - (PIO2_1+PIO2_2) -const PIO2_2T: f64 = 2.02226624879595063154e-21; /* 0x3BA3198A, 0x2E037073 */ -/// third 33 bits of pi/2 -const PIO2_3: f64 = 2.02226624871116645580e-21; /* 0x3BA3198A, 0x2E000000 */ -/// pi/2 - (PIO2_1+PIO2_2+PIO2_3) -const PIO2_3T: f64 = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */ - // return the remainder of x rem pi/2 in y[0]+y[1] // use rem_pio2_large() for large x // // caller must handle the case when reduction is not needed: |x| ~<= pi/4 */ #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub(crate) fn rem_pio2(x: f64) -> (i32, f64, f64) { - return super::generic::rem_frac_pi_2(x); + return super::generic::rem_pio2(x); } #[cfg(test)] diff --git a/libm/src/math/support/mod.rs b/libm/src/math/support/mod.rs index 1f0b59720..8db3a6d39 100644 --- a/libm/src/math/support/mod.rs +++ b/libm/src/math/support/mod.rs @@ -11,7 +11,7 @@ pub use big::{i256, u256}; pub use env::{FpResult, Round, Status}; #[allow(unused_imports)] pub use float_traits::{DFloat, Float, HFloat, HalfRep, IntTy}; -pub(crate) use float_traits::{f32_from_bits, f64_from_bits}; +pub(crate) use float_traits::{f32_from_bits, f64_from_bits, f64_to_bits}; #[cfg(f16_enabled)] #[allow(unused_imports)] pub use hex_float::hf16; From 927d430fa22f3cc19019638ee8914c97c1612c99 Mon Sep 17 00:00:00 2001 From: Trevor Gross Date: Wed, 23 Apr 2025 05:17:33 +0000 Subject: [PATCH 3/5] msrv and format --- libm/src/math/generic/generated/mod.rs | 2 +- libm/src/math/generic/generated/rem_pio2.rs | 6 +++--- libm/src/math/generic/rem_pio2.rs | 10 +++++++++- 3 files changed, 13 insertions(+), 5 deletions(-) diff --git a/libm/src/math/generic/generated/mod.rs b/libm/src/math/generic/generated/mod.rs index 4f36e0561..d5116219c 100644 --- a/libm/src/math/generic/generated/mod.rs +++ b/libm/src/math/generic/generated/mod.rs @@ -1,3 +1,3 @@ // autogenerated -mod rem_pio2; \ No newline at end of file +mod rem_pio2; diff --git a/libm/src/math/generic/generated/rem_pio2.rs b/libm/src/math/generic/generated/rem_pio2.rs index 035f5c4ef..8549a1568 100644 --- a/libm/src/math/generic/generated/rem_pio2.rs +++ b/libm/src/math/generic/generated/rem_pio2.rs @@ -1,6 +1,5 @@ -use crate::support::HalfRep; use super::super::rem_pio2::RemPio2Support; -use crate::support::f64_to_bits; +use crate::support::{HalfRep, f64_to_bits}; impl RemPio2Support for f64 { const TO_INT: Self = 1.5 / f64::EPSILON; @@ -19,7 +18,8 @@ impl RemPio2Support for f64 { const FRAC_3PI_2_HI: HalfRep = (f64_to_bits(hf64!("0x4.b65f1fccc8748p+0")) >> 32) as u32; const TAU_HI: HalfRep = (f64_to_bits(hf64!("0x6.487ed5110b46p+0")) >> 32) as u32; const FRAC_PI_2_HI: HalfRep = 0x921fb; - const FRAC_2_POW_20_PI_2: HalfRep = (f64_to_bits(hf64!("0x1.921fb54442d18p+20")) >> 32) as u32; + const FRAC_2_POW_20_PI_2: HalfRep = + (f64_to_bits(hf64!("0x1.921fb54442d18p+20")) >> 32) as u32; const MAGIC_F: Self = hf64!("0x1p24"); diff --git a/libm/src/math/generic/rem_pio2.rs b/libm/src/math/generic/rem_pio2.rs index 7b3d9fa5c..d6c6c444a 100644 --- a/libm/src/math/generic/rem_pio2.rs +++ b/libm/src/math/generic/rem_pio2.rs @@ -4,7 +4,11 @@ use core::f64::consts; use crate::support::{CastFrom, CastInto, DInt, Float, HInt, HalfRep, Int, MinInt}; -pub(crate) trait RemPio2Support: Float> { +pub(crate) trait RemPio2Support: Float +where + Self::Int: DInt, + ::H: Int, +{ const TO_INT: Self; /// 53 bits of 2/pi const INV_PIO2: Self; @@ -43,6 +47,8 @@ where F: RemPio2Support, F: CastInto, HalfRep: Int + MinInt>, + F::Int: DInt, + ::H: Int, { let ix: HalfRep = x.abs().to_bits().hi(); let pos = x.is_sign_positive(); @@ -173,6 +179,8 @@ where F: RemPio2Support, F: CastInto, HalfRep: Int, + F::Int: DInt, + ::H: Int, { /* rint(x/(pi/2)), Assume round-to-nearest. */ let tmp = x * F::INV_PIO2 + F::TO_INT; From 975a47d0ed5750efcf905ea85b86540b97b02cc3 Mon Sep 17 00:00:00 2001 From: Trevor Gross Date: Wed, 23 Apr 2025 05:55:17 +0000 Subject: [PATCH 4/5] WIP --- libm-test/tests/multiprecision.rs | 64 ++++++++++++++++++++++++++++++- libm/src/math/generic/rem_pio2.rs | 11 ++++-- 2 files changed, 70 insertions(+), 5 deletions(-) diff --git a/libm-test/tests/multiprecision.rs b/libm-test/tests/multiprecision.rs index 80b2c7868..6f039f2ed 100644 --- a/libm-test/tests/multiprecision.rs +++ b/libm-test/tests/multiprecision.rs @@ -3,7 +3,7 @@ #![cfg(feature = "build-mpfr")] use libm_test::generate::{case_list, edge_cases, random, spaced}; -use libm_test::mpfloat::MpOp; +use libm_test::mpfloat::{MpFloat, MpOp}; use libm_test::{CheckBasis, CheckCtx, CheckOutput, GeneratorKind, MathOp, TupleCall}; const BASIS: CheckBasis = CheckBasis::Mpfr; @@ -77,3 +77,65 @@ libm_macros::for_each_function! { nextafterf, ], } + +// fn mp_runner_rem_pio2(ctx: &CheckCtx, cases: impl Iterator) { +// let x = MpFloat::new(prec) + +// // let mut mp_vals = Op::new_mp(); +// for input in cases { +// // let mp_res = Op::run(&mut mp_vals, input); +// let crate_res = input.call_intercept_panics(Op::ROUTINE); + +// crate_res.validate(mp_res, input, ctx).unwrap(); +// } +// } + +// macro_rules! mp_tests_rem_pio2 { +// ( +// fn_name: $fn_name:ident, +// attrs: [$($attr:meta),*], +// ) => { +// paste::paste! { +// #[test] +// $(#[$attr])* +// fn [< mp_case_list_ $fn_name >]() { +// type Op = libm_test::op::sin::Routine; +// let ctx = CheckCtx::new(Op::IDENTIFIER, BASIS, GeneratorKind::List); +// let cases = case_list::get_test_cases_basis::(&ctx).0; +// mp_runner_rem_pio2(&ctx, cases); +// } + +// #[test] +// $(#[$attr])* +// fn [< mp_random_ $fn_name >]() { +// type Op = libm_test::op::sin::Routine; +// let ctx = CheckCtx::new(Op::IDENTIFIER, BASIS, GeneratorKind::Random); +// let cases = random::get_test_cases::<::RustArgs>(&ctx).0; +// mp_runner_rem_pio2(&ctx, cases); +// } + +// #[test] +// $(#[$attr])* +// fn [< mp_edge_case_ $fn_name >]() { +// type Op = libm_test::op::sin::Routine; +// let ctx = CheckCtx::new(Op::IDENTIFIER, BASIS, GeneratorKind::EdgeCases); +// let cases = edge_cases::get_test_cases::(&ctx).0; +// mp_runner_rem_pio2(&ctx, cases); +// } + +// #[test] +// $(#[$attr])* +// fn [< mp_quickspace_ $fn_name >]() { +// type Op = libm_test::op::sin::Routine; +// let ctx = CheckCtx::new(Op::IDENTIFIER, BASIS, GeneratorKind::QuickSpaced); +// let cases = spaced::get_test_cases::(&ctx).0; +// mp_runner_rem_pio2(&ctx, cases); +// } +// } +// }; +// } + +// mp_tests_rem_pio2! { +// fn_name: rem_pio2, +// attrs: [], +// } diff --git a/libm/src/math/generic/rem_pio2.rs b/libm/src/math/generic/rem_pio2.rs index d6c6c444a..d82dd4f35 100644 --- a/libm/src/math/generic/rem_pio2.rs +++ b/libm/src/math/generic/rem_pio2.rs @@ -42,6 +42,9 @@ where fn large(x: &[Self], y: &mut [Self], e0: i32, prec: usize) -> i32; } +/// Return `x % pi/2` as `y[0] + y[1]`. +/// +/// Caller must handle the case when reduction is not needed: `|x| ~<= pi/4`. pub(crate) fn rem_pio2(x: F) -> (i32, F, F) where F: RemPio2Support, @@ -56,14 +59,14 @@ where if ix <= F::FRAC_5PI_4_HI { /* |x| ~<= 5pi/4 */ if (ix & F::SIG_MASK.hi()) == F::FRAC_PI_2_HI { - /* |x| ~= pi/2 or 2pi/2 */ - return medium(x, ix); /* cancellation -- use medium case */ + // |x| ~= pi/2 or 2pi/2 + return medium(x, ix); // cancellation -- use medium case } if ix <= F::FRAC_3PI_4_HI { - /* |x| ~<= 3pi/4 */ + // |x| ~<= 3pi/4 if pos { - let z = x - F::PIO2_1; /* one round good to 85 bits */ + let z = x - F::PIO2_1; // one round good to 85 bits for f64 let y0 = z - F::PIO2_1T; let y1 = (z - y0) - F::PIO2_1T; return (1, y0, y1); From 25e21bf46c3f51b35cd629a3b8c15a2d2b90d3d1 Mon Sep 17 00:00:00 2001 From: Trevor Gross Date: Wed, 23 Apr 2025 15:05:17 -0400 Subject: [PATCH 5/5] Update script --- etc/generate-constants.jl | 21 +++++++++++---------- libm/src/math/generic/generated/rem_pio2.rs | 2 +- 2 files changed, 12 insertions(+), 11 deletions(-) diff --git a/etc/generate-constants.jl b/etc/generate-constants.jl index 30876740b..2be4871ce 100644 --- a/etc/generate-constants.jl +++ b/etc/generate-constants.jl @@ -75,14 +75,15 @@ function update_rem_pio2(cfg::Config)::String to_bits = "$ty_name::to_bits" end - hi_int(x::BigFloat) = @sprintf "(%s(hf%d!(\"%a\")) >> %d) as u%d" to_bits info.bits x halfbits halfbits + bigf_hi(x::BigFloat) = @sprintf "(%s(hf%d!(\"%a\")) >> %d) as u%d" to_bits info.bits x halfbits halfbits + bigf_hex(x::BigFloat) = @sprintf "hf%d!(\"%a\")" info.bits x setprecision(ty_info(fty).sig_bits + 1) ty_impl = """ impl RemPio2Support for $ty_name { const TO_INT: Self = 1.5 / $ty_name::EPSILON; - const INV_PIO2: Self = 6.36619772367581382433e-01; + const INV_PIO2: Self = $(bigf_hex(2 / big(pi))); const PIO2_1: Self = 1.57079632673412561417e+00; const PIO2_1T: Self = 6.07710050650619224932e-11; const PIO2_2: Self = 6.07710050630396597660e-11; @@ -90,16 +91,16 @@ function update_rem_pio2(cfg::Config)::String const PIO2_3: Self = 2.02226624871116645580e-21; const PIO2_3T: Self = 8.47842766036889956997e-32; - const FRAC_5PI_4_HI: HalfRep = $(hi_int(big(pi)*5/4)); - const FRAC_3PI_4_HI: HalfRep = $(hi_int(big(pi)*3/4)); - const FRAC_9PI_4_HI: HalfRep = $(hi_int(big(pi)*9/4)); - const FRAC_7PI_4_HI: HalfRep = $(hi_int(big(pi)*7/4)); - const FRAC_3PI_2_HI: HalfRep = $(hi_int(big(pi)*3/2)); - const TAU_HI: HalfRep = $(hi_int(big(pi)*2)); + const FRAC_5PI_4_HI: HalfRep = $(bigf_hi(5 * big(pi) / 4)); + const FRAC_3PI_4_HI: HalfRep = $(bigf_hi(3 * big(pi) / 4)); + const FRAC_9PI_4_HI: HalfRep = $(bigf_hi(9 * big(pi) / 4)); + const FRAC_7PI_4_HI: HalfRep = $(bigf_hi(7 * big(pi) / 4)); + const FRAC_3PI_2_HI: HalfRep = $(bigf_hi(3 * big(pi) / 2)); + const TAU_HI: HalfRep = $(bigf_hi(2 * big(pi))); const FRAC_PI_2_HI: HalfRep = 0x921fb; - const FRAC_2_POW_20_PI_2: HalfRep = $(hi_int((big(2)^20) * pi / 2)); + const FRAC_2_POW_20_PI_2: HalfRep = $(bigf_hi((big(2)^20) * pi / 2)); - const MAGIC_F: Self = hf64!("0x1p24"); + const MAGIC_F: Self = hf$(info.bits)!("0x1p24"); fn large(x: &[Self], y: &mut [Self], e0: i32, prec: usize) -> i32 { super::super::super::rem_pio2_large(x, y, e0, prec) diff --git a/libm/src/math/generic/generated/rem_pio2.rs b/libm/src/math/generic/generated/rem_pio2.rs index 8549a1568..db3566b8f 100644 --- a/libm/src/math/generic/generated/rem_pio2.rs +++ b/libm/src/math/generic/generated/rem_pio2.rs @@ -3,7 +3,7 @@ use crate::support::{HalfRep, f64_to_bits}; impl RemPio2Support for f64 { const TO_INT: Self = 1.5 / f64::EPSILON; - const INV_PIO2: Self = 6.36619772367581382433e-01; + const INV_PIO2: Self = hf64!("0xa.2f9836e4e4418p-4"); const PIO2_1: Self = 1.57079632673412561417e+00; const PIO2_1T: Self = 6.07710050650619224932e-11; const PIO2_2: Self = 6.07710050630396597660e-11;