Add a generic rempio2

etc/generate-constants.jl

@@ -0,0 +1,121 @@
+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
+
+        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 = $(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;
+                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<Self> = $(bigf_hi(5 * big(pi) / 4));
+                const FRAC_3PI_4_HI: HalfRep<Self> = $(bigf_hi(3 * big(pi) / 4));
+                const FRAC_9PI_4_HI: HalfRep<Self> = $(bigf_hi(9 * big(pi) / 4));
+                const FRAC_7PI_4_HI: HalfRep<Self> = $(bigf_hi(7 * big(pi) / 4));
+                const FRAC_3PI_2_HI: HalfRep<Self> = $(bigf_hi(3 * big(pi) / 2));
+                const TAU_HI: HalfRep<Self> = $(bigf_hi(2 * big(pi)));
+                const FRAC_PI_2_HI: HalfRep<Self> = 0x921fb;
+                const FRAC_2_POW_20_PI_2: HalfRep<Self> = $(bigf_hi((big(2)^20) * pi / 2));
+
+                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)
+                }
+            }
+            """
+        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()

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);

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<F: Float, Args>(ctx: &CheckCtx, cases: impl Iterator<Item = Args>) {
+//     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::<Op>(&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::<<Op as MathOp>::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::<Op>(&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::<Op>(&ctx).0;
+//                 mp_runner_rem_pio2(&ctx, cases);
+//             }
+//         }
+//     };
+// }
+
+// mp_tests_rem_pio2! {
+//     fn_name: rem_pio2,
+//     attrs: [],
+// }

libm/src/math/generic/generated/mod.rs

@@ -0,0 +1,3 @@
+// autogenerated
+
+mod rem_pio2;

libm/src/math/generic/generated/rem_pio2.rs

@@ -0,0 +1,29 @@
+use super::super::rem_pio2::RemPio2Support;
+use crate::support::{HalfRep, f64_to_bits};
+
+impl RemPio2Support for f64 {
+    const TO_INT: Self = 1.5 / f64::EPSILON;
+    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;
+    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<Self> = (f64_to_bits(hf64!("0x3.ed4f452aa70bcp+0")) >> 32) as u32;
+    const FRAC_3PI_4_HI: HalfRep<Self> = (f64_to_bits(hf64!("0x2.5b2f8fe6643a4p+0")) >> 32) as u32;
+    const FRAC_9PI_4_HI: HalfRep<Self> = (f64_to_bits(hf64!("0x7.118eafb32caecp+0")) >> 32) as u32;
+    const FRAC_7PI_4_HI: HalfRep<Self> = (f64_to_bits(hf64!("0x5.7f6efa6ee9dd4p+0")) >> 32) as u32;
+    const FRAC_3PI_2_HI: HalfRep<Self> = (f64_to_bits(hf64!("0x4.b65f1fccc8748p+0")) >> 32) as u32;
+    const TAU_HI: HalfRep<Self> = (f64_to_bits(hf64!("0x6.487ed5110b46p+0")) >> 32) as u32;
+    const FRAC_PI_2_HI: HalfRep<Self> = 0x921fb;
+    const FRAC_2_POW_20_PI_2: HalfRep<Self> =
+        (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)
+    }
+}

libm/src/math/generic/mod.rs

@@ -13,6 +13,8 @@ mod fmin;
 mod fminimum;
 mod fminimum_num;
 mod fmod;
+mod generated;
+mod rem_pio2;
 mod rint;
 mod round;
 mod scalbn;
@@ -31,6 +33,7 @@ pub use fmin::fmin;
 pub use fminimum::fminimum;
 pub use fminimum_num::fminimum_num;
 pub use fmod::fmod;
+pub(crate) use rem_pio2::rem_pio2;
 pub use rint::rint_round;
 pub use round::round;
 pub use scalbn::scalbn;

libm/src/math/generic/rem_pio2.rs

@@ -0,0 +1,223 @@
+#![allow(unused)]
+
+use core::f64::consts;
+
+use crate::support::{CastFrom, CastInto, DInt, Float, HInt, HalfRep, Int, MinInt};
+
+pub(crate) trait RemPio2Support: Float
+where
+    Self::Int: DInt,
+    <Self::Int as DInt>::H: Int,
+{
+    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<Self>;
+    const FRAC_3PI_4_HI: HalfRep<Self>;
+    const FRAC_9PI_4_HI: HalfRep<Self>;
+    const FRAC_7PI_4_HI: HalfRep<Self>;
+    const FRAC_3PI_2_HI: HalfRep<Self>;
+    /// 2pi
+    const TAU_HI: HalfRep<Self>;
+    const FRAC_PI_2_HI: HalfRep<Self>;
+    /// (2^20)(pi/2)
+    const FRAC_2_POW_20_PI_2: HalfRep<Self>;
+
+    const MAGIC: u32 = 23;
+    const MAGIC_F: Self;
+
+    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<F>(x: F) -> (i32, F, F)
+where
+    F: RemPio2Support,
+    F: CastInto<i32>,
+    HalfRep<F>: Int + MinInt<Unsigned: Int<OtherSign: Int>>,
+    F::Int: DInt,
+    <F::Int as DInt>::H: Int,
+{
+    let ix: HalfRep<F> = 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 for f64
+                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;
+
+    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::<F>::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<F>(x: F, ix: HalfRep<F>) -> (i32, F, F)
+where
+    F: RemPio2Support,
+    F: CastInto<i32>,
+    HalfRep<F>: Int,
+    F::Int: DInt,
+    <F::Int as DInt>::H: 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::<F>::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)
+}

libm/src/math/rem_pio2.rs

@@ -1,3 +1,4 @@
+#![allow(unused)]
 // origin: FreeBSD /usr/src/lib/msun/src/e_rem_pio2.c
 //
 // ====================================================
@@ -12,181 +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) {
-    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_pio2(x);
 }
 
 #[cfg(test)]

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<F> = <<F as Float>::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);

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);

libm/src/math/support/mod.rs

@@ -10,8 +10,8 @@ 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(crate) use float_traits::{f32_from_bits, f64_from_bits};
+pub use float_traits::{DFloat, Float, HFloat, HalfRep, IntTy};
+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;