deep_time/math/

sin.rs

// origin: FreeBSD /usr/src/lib/msun/src/s_sin.c */
//
// ====================================================
// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
//
// Developed at SunPro, a Sun Microsystems, Inc. business.
// Permission to use, copy, modify, and distribute this
// software is freely granted, provided that this notice
// is preserved.
// ====================================================

#![allow(clippy::indexing_slicing)]
#![allow(clippy::excessive_precision)]
#![allow(clippy::approx_constant)]
#![allow(clippy::eq_op)]

use super::{k_cos, k_sin, rem_pio2};
use crate::Real;

/// Computes the sine of `x` (in radians).
///
/// This is a `const fn` implementation based on argument reduction
/// followed by a polynomial approximation.
///
/// ### Testing
///
/// The following tests have been performed:
///
/// - Maximum observed error of ≤ 1 ULP, measured over 5,000,000 random
///   samples across a wide dynamic range.
/// - Edge cases, including ±0, ±π/2, ±π, subnormal numbers, infinity,
///   and NaN.
/// - Monotonicity testing in both increasing and decreasing regions.
/// - High-density testing near π/2.
/// - Testing near multiples of π.
/// - Hard argument reduction cases.
/// - Statistical random testing across multiple magnitude scales.
/// - Compile-time evaluation via `const fn`.
///
/// This implementation is intended for use in `no_std` and embedded
/// environments.
pub const fn sin(x: Real) -> Real {
    /* High word of x. */
    let ix = (Real::to_bits(x) >> 32) as u32 & 0x7fffffff;

    /* |x| ~< pi/4 */
    if ix <= 0x3fe921fb {
        if ix < 0x3e500000 {
            /* |x| < 2**-26 */
            return x;
        }
        return k_sin(x, 0.0, 0);
    }

    /* sin(Inf or NaN) is NaN */
    if ix >= 0x7ff00000 {
        return x - x;
    }

    /* argument reduction needed */
    let (n, y0, y1) = rem_pio2(x);
    match n & 3 {
        0 => k_sin(y0, y1, 1),
        1 => k_cos(y0, y1),
        2 => -k_sin(y0, y1, 1),
        _ => -k_cos(y0, y1),
    }
}

#[cfg(all(test, feature = "std"))]
mod sin_tests {
    use super::*;

    // =====================================================================
    // Test Helper: Simple PRNG (only used in tests)
    // =====================================================================

    #[allow(dead_code)]
    struct Rng(u64);

    #[allow(dead_code)]
    impl Rng {
        fn new(seed: u64) -> Self {
            Self(seed)
        }

        #[inline]
        fn next_u64(&mut self) -> u64 {
            self.0 = self.0.wrapping_add(0x9E3779B97F4A7C15);
            let mut z = self.0;
            z = (z ^ (z >> 30)).wrapping_mul(0xBF58476D1CE4E5B9);
            z = (z ^ (z >> 27)).wrapping_mul(0x94D049BB133111EB);
            z ^ (z >> 31)
        }

        /// Returns a value roughly in `[-scale, scale]`
        fn next_f64(&mut self, scale: f64) -> f64 {
            let bits = self.next_u64();
            let mantissa = bits & 0x000f_ffff_ffff_ffff;
            let f = f64::from_bits(0x3ff0_0000_0000_0000 | mantissa);
            (f - 1.0) * 2.0 * scale - scale
        }
    }

    // =====================================================================
    // Utility: ULP difference
    // =====================================================================

    #[allow(dead_code)]
    fn ulp_diff(a: f64, b: f64) -> u64 {
        if a.is_nan() && b.is_nan() {
            return 0;
        }
        if a.is_infinite() && b.is_infinite() && a.signum() == b.signum() {
            return 0;
        }
        a.to_bits().abs_diff(b.to_bits())
    }

    // =====================================================================
    // Tests
    // =====================================================================

    #[test]
    fn sin_edge_cases() {
        let pi = std::f64::consts::PI;
        let pi_over_2 = std::f64::consts::FRAC_PI_2;

        assert_eq!(sin(0.0), 0.0);
        assert_eq!(sin(-0.0), -0.0);
        assert!((sin(1.0) - 1.0f64.sin()).abs() < 1e-15);

        // Multiples of π/2
        assert!((sin(pi_over_2) - 1.0).abs() < 1e-14);
        assert!((sin(-pi_over_2) + 1.0).abs() < 1e-14);
        assert!((sin(pi) - 0.0).abs() < 1e-14);
        assert!((sin(3.0 * pi_over_2) + 1.0).abs() < 1e-14);

        // Very small
        assert_eq!(sin(1e-300), 1e-300);

        // Large values
        let large = 1e10;
        let diff = (sin(large) - large.sin()).abs();
        assert!(diff < 1e-6 || sin(large).is_nan());

        let neg_large = -1e8;
        assert!((sin(neg_large) - neg_large.sin()).abs() < 1e-5);

        // Extremely large
        let huge = 1e300;
        let s = sin(huge);
        assert!(s.is_nan() || s.abs() <= 1.0 + 1e-9);
    }

    #[test]
    fn sin_very_large_arguments() {
        // This test specifically exercises rem_pio2_large and its recompute logic.
        // These values are large enough that they go through the big-argument path.
        let large_values: &[f64] = &[
            1e40,
            1e80,
            1e120,
            1e160,
            1e200,
            -1e50,
            -1e100,
            -1e150,
            1e10 + std::f64::consts::PI * 1e8, // large + offset
            -1e12 - std::f64::consts::PI * 1e7,
        ];

        for &x in large_values {
            let our = sin(x);
            let std_val = x.sin();

            if our.is_nan() && std_val.is_nan() {
                continue;
            }

            // We allow a slightly looser tolerance here because these are extreme values
            let diff = (our - std_val).abs();
            assert!(
                diff < 1e-5 || our.is_nan(),
                "sin mismatch on very large argument at x = {}: diff = {}",
                x,
                diff
            );
        }
    }

    #[test]
    fn sin_identities() {
        let x = 1.23456789;

        assert!((sin(-x) + sin(x)).abs() < 1e-15);
        assert!((sin(x + 2.0 * std::f64::consts::PI) - sin(x)).abs() < 1e-9);
    }

    #[test]
    fn sin_monotonicity() {
        let tol = 1e-12;

        // Region 1: Clearly increasing
        let mut prev = sin(-1.0);
        for i in 1..100_000 {
            let x = -1.0 + (i as f64) * 2e-5;
            let y = sin(x);
            assert!(y + tol >= prev, "Non-monotonic (increasing) at x = {}", x);
            prev = y;
        }

        // Region 2: Clearly decreasing
        prev = sin(std::f64::consts::FRAC_PI_2 + 0.1);
        for i in 1..100_000 {
            let x = std::f64::consts::FRAC_PI_2 + 0.1 + (i as f64) * 2e-5;
            let y = sin(x);
            assert!(y + tol <= prev, "Non-monotonic (decreasing) at x = {}", x);
            prev = y;
        }
    }

    #[test]
    fn sin_very_small_values() {
        // Test that sin(x) ≈ x for very small values
        for i in 0..30 {
            let x = 1e-20 * (i as f64 + 1.0);
            assert_eq!(sin(x), x, "Failed at x = {}", x);
        }

        // Additional small values
        assert_eq!(sin(1e-300), 1e-300);
        assert_eq!(sin(-1e-250), -1e-250);
    }

    #[test]
    fn sin_hard_reduction_cases() {
        let cases: &[f64] = &[
            1.5707963267948966,
            4.71238898038469,
            1e10 + 0.5,
            std::f64::consts::PI * 1e8,
            -std::f64::consts::PI * 1e7 + 1e-9,
            1e20,
            -1e20,
        ];

        for &x in cases {
            let our = sin(x);
            let std = x.sin();
            let diff = (our - std).abs();
            assert!(
                diff < 1e-8 || our.is_nan(),
                "Hard reduction case failed at x = {}: diff = {}",
                x,
                diff
            );
        }
    }

    #[test]
    fn sin_near_pi_over_2() {
        let pi_over_2 = std::f64::consts::FRAC_PI_2;

        for i in 0..100_000 {
            let delta = (i as f64 - 50_000.0) * 1e-11;
            let x = pi_over_2 + delta;
            let our = sin(x);
            let expected = x.sin();
            let diff = (our - expected).abs();
            assert!(
                diff < 1e-10,
                "Large error near π/2 at x = {}: diff = {}",
                x,
                diff
            );
        }
    }

    #[test]
    fn sin_near_multiples_of_pi() {
        let pi = std::f64::consts::PI;

        for k in -10i32..=10 {
            let base = (k as f64) * pi;

            // Test slightly above and below k*π
            for &delta in &[1e-9, 1e-8, -1e-9, -1e-8] {
                let x = base + delta;
                let our = sin(x);
                let std = x.sin();
                let diff = (our - std).abs();
                assert!(
                    diff < 1e-9 || our.is_nan(),
                    "Error near {}π at x = {}: diff = {}",
                    k,
                    x,
                    diff
                );
            }
        }
    }

    #[test]
    fn sin_max_ulp_error() {
        let mut max_ulp: u64 = 0;
        let mut worst_x = 0.0f64;
        let mut rng = Rng::new(0xdead_beef_cafe_babe);

        for _ in 0..5_000_000 {
            let scale = if rng.next_u64() & 1 == 0 { 1e6 } else { 1e12 };
            let x = rng.next_f64(scale);
            let our = sin(x);
            let std = x.sin();

            if our.is_nan() || std.is_nan() {
                continue;
            }

            let ulp = ulp_diff(our, std);
            if ulp > max_ulp {
                max_ulp = ulp;
                worst_x = x;
            }
        }

        // println!("Maximum observed ULP error: {} at x ≈ {}", max_ulp, worst_x);
        assert!(
            max_ulp <= 1,
            "sin() maximum error too high: {} ULPs at x = {}",
            max_ulp,
            worst_x
        );
    }

    #[test]
    fn sin_random_many() {
        let mut rng = Rng::new(0x1234_5678_9abc_def0);
        let scales = [1.0, 10.0, 100.0, 1_000.0, 1e6, 1e8, 1e10];

        for &scale in &scales {
            for _ in 0..150_000 {
                let x = rng.next_f64(scale);
                let our = sin(x);
                let std_val = x.sin();

                if our.is_nan() && std_val.is_nan() {
                    continue;
                }

                // Adaptive tolerance
                let diff = if x.abs() < 1e-8 {
                    (our - std_val).abs()
                } else {
                    (our - std_val).abs() / x.abs().max(1.0)
                };

                let tol = if x.abs() < 1e-8 { 1e-14 } else { 5e-12 };

                assert!(
                    diff < tol,
                    "sin mismatch at x = {}: our={}, std={}, diff={}",
                    x,
                    our,
                    std_val,
                    diff
                );
            }
        }
    }

    /// Compile-time check that `sin` can be used in const contexts.
    const _: () = {
        let _ = sin(0.0);
        let _ = sin(1.0);
        let _ = sin(-std::f64::consts::PI);
    };
}