From 8313a32b8152165dc38e91786635e41face232f8 Mon Sep 17 00:00:00 2001 From: Jeong YunWon Date: Tue, 13 Jan 2026 22:04:18 +0900 Subject: [PATCH 1/2] more edge values --- .gitignore | 1 + src/math/integer.rs | 33 +++++++++++++++++++---- src/test.rs | 64 ++++++++++++++++++++++++++++++++++++++------- 3 files changed, 83 insertions(+), 15 deletions(-) diff --git a/.gitignore b/.gitignore index 96ef6c0..66113fa 100644 --- a/.gitignore +++ b/.gitignore @@ -1,2 +1,3 @@ /target Cargo.lock +cpython diff --git a/src/math/integer.rs b/src/math/integer.rs index 3638ac9..3c11080 100644 --- a/src/math/integer.rs +++ b/src/math/integer.rs @@ -715,7 +715,7 @@ mod tests { use pyo3::prelude::*; /// Edge i64 values for testing integer math functions (gcd, lcm, isqrt, factorial, comb, perm) - const EDGE_I64: [i64; 24] = [ + const EDGE_I64: [i64; 44] = [ // Zero and small values 0, 1, @@ -726,27 +726,50 @@ mod tests { 7, 13, 97, - // Powers of 2 + 127, // table boundary in comb/perm + 128, // Table boundary + 1 + // Powers of 2 and boundaries 64, + 63, // 2^6 - 1 + 65, // 2^6 + 1 1024, - 65536, + 65535, // 2^16 - 1 + 65536, // 2^16 + 65537, // 2^16 + 1 (Fermat prime) // Factorial-relevant + 12, // 12! = 479001600 fits in u32 + 13, // 13! overflows u32 20, // 20! fits in u64 21, // 21! overflows u64 + 170, // factorial(170) is the largest that fits in f64 + 171, // factorial(171) overflows f64 + // Comb/perm algorithm switching points + 34, // FAST_COMB_LIMITS1 boundary + 35, // Large values 1_000_000, -1_000_000, i32::MAX as i64, + i32::MAX as i64 + 1, i32::MIN as i64, + i32::MIN as i64 - 1, // Near i64 bounds i64::MAX, i64::MIN, i64::MAX - 1, i64::MIN + 1, // Square root boundaries - (1i64 << 31) - 1, // sqrt fits in u32 - 1i64 << 32, // sqrt boundary + (1i64 << 15) - 1, // 32767, sqrt = 181 + 1i64 << 16, // 65536, sqrt = 256 (exact) + (1i64 << 31) - 1, // sqrt fits in u16 + 1i64 << 32, // sqrt boundary (exact power of 2) + (1i64 << 32) - 1, // near sqrt boundary + (1i64 << 32) + 1, // just above sqrt boundary (1i64 << 62) - 1, // large but valid for isqrt + 1i64 << 62, // exact power of 2 + // Near perfect squares + 99, // sqrt(99) = 9.949... + 101, // sqrt(101) = 10.049... ]; fn test_gcd_impl(args: &[i64]) { diff --git a/src/test.rs b/src/test.rs index cf9ebf7..d33f0c6 100644 --- a/src/test.rs +++ b/src/test.rs @@ -3,7 +3,7 @@ use pyo3::{Python, prelude::*}; /// Edge values for testing floating-point functions. /// Includes: zeros, infinities, various NaNs, subnormals, and values at different scales. -pub(crate) const EDGE_VALUES: [f64; 30] = [ +pub(crate) const EDGE_VALUES: [f64; 64] = [ // Zeros 0.0, -0.0, @@ -15,35 +15,79 @@ pub(crate) const EDGE_VALUES: [f64; 30] = [ -f64::NAN, // Additional NaN with different payload (quiet NaN with payload 1) f64::from_bits(0x7FF8_0000_0000_0001_u64), + // Signaling NaN (sNaN) - may trigger FP exceptions on some platforms + f64::from_bits(0x7FF0_0000_0000_0001_u64), // Subnormal (denormalized) values - f64::MIN_POSITIVE * 0.5, // smallest subnormal + f64::MIN_POSITIVE * 0.5, -f64::MIN_POSITIVE * 0.5, + 5e-324, + -5e-324, // Boundary values - f64::MIN_POSITIVE, // smallest positive normal - f64::MAX, // largest finite - f64::MIN, // most negative finite (not smallest!) + f64::MIN_POSITIVE, + f64::MAX, + f64::MIN, // Near-infinity large values f64::MAX * 0.5, -f64::MAX * 0.5, 1e308, -1e308, + // Overflow/underflow thresholds for exp + 710.0, + -745.0, // Small scale 1e-10, -1e-10, 1e-300, + -1e-300, // Normal scale 1.0, -1.0, 0.5, -0.5, 2.0, - // Trigonometric special values (where sin/cos/tan have exact or near-zero results) - std::f64::consts::PI, // sin(PI) ≈ 0 + -2.0, + 3.0, // for cbrt + -3.0, + // Values near 1.0 (log, expm1, log1p, acosh boundary) + 1.0 - 1e-15, + 1.0 + 1e-15, + f64::EPSILON, + 1.0 - f64::EPSILON, + 1.0 + f64::EPSILON, + // asin/acos domain boundaries [-1, 1] + 1.0000000000000002, // just outside domain (1 + eps) + -1.0000000000000002, + // atanh domain boundaries (-1, 1) + 0.9999999999999999, // just inside domain + -0.9999999999999999, + // log1p domain boundary (> -1) + -0.9999999999999999, // just above -1 + -1.0 + 1e-15, // very close to -1 + // gamma/lgamma poles (negative integers) + -1.0, + -2.0, + -3.0, + -0.5, // gamma(-0.5) = -2*sqrt(pi) + // Mathematical constants + std::f64::consts::E, + std::f64::consts::LN_2, + std::f64::consts::LOG10_E, + // Trigonometric special values + std::f64::consts::PI, -std::f64::consts::PI, - std::f64::consts::FRAC_PI_2, // cos(PI/2) ≈ 0 + std::f64::consts::FRAC_PI_2, -std::f64::consts::FRAC_PI_2, - std::f64::consts::FRAC_PI_4, // tan(PI/4) = 1 - std::f64::consts::TAU, // sin(2*PI) ≈ 0, cos(2*PI) = 1 + std::f64::consts::FRAC_PI_4, + std::f64::consts::TAU, + 1.5 * std::f64::consts::PI, // 3π/2 + // Large values for trig (precision loss) + 1e15, + -1e15, + // Near-integer values (ceil, floor, trunc, round) + 0.49999999999999994, + 0.50000000000000006, + -0.49999999999999994, + -0.50000000000000006, ]; pub(crate) fn unwrap<'py>( From 0c83ff1012066d568b877a0887622a396d063218 Mon Sep 17 00:00:00 2001 From: Jeong YunWon Date: Thu, 15 Jan 2026 09:25:10 +0900 Subject: [PATCH 2/2] try --- src/cmath/exponential.rs | 5 +++-- src/cmath/misc.rs | 37 ++++++++++++++++++++++++++++++++----- src/math/integer.rs | 14 +++++++------- 3 files changed, 42 insertions(+), 14 deletions(-) diff --git a/src/cmath/exponential.rs b/src/cmath/exponential.rs index f306065..9b3e793 100644 --- a/src/cmath/exponential.rs +++ b/src/cmath/exponential.rs @@ -4,7 +4,7 @@ use super::{ CM_LARGE_DOUBLE, CM_LOG_LARGE_DOUBLE, INF, M_LN2, N, P, P12, P14, P34, U, c, special_type, special_value, }; -use crate::{Error, Result, m}; +use crate::{Error, Result, m, mul_add}; use num_complex::Complex64; // Local constants @@ -160,7 +160,8 @@ pub(crate) fn ln(z: Complex64) -> Result { if (0.71..=1.73).contains(&h) { let am = if ax > ay { ax } else { ay }; // max(ax, ay) let an = if ax > ay { ay } else { ax }; // min(ax, ay) - m::log1p((am - 1.0) * (am + 1.0) + an * an) / 2.0 + let log1p_arg = mul_add(am - 1.0, am + 1.0, an * an); + m::log1p(log1p_arg) / 2.0 } else { m::log(h) } diff --git a/src/cmath/misc.rs b/src/cmath/misc.rs index 69514e2..f7136fc 100644 --- a/src/cmath/misc.rs +++ b/src/cmath/misc.rs @@ -28,14 +28,41 @@ pub fn phase(z: Complex64) -> Result { } } +#[inline] +fn c_abs_raw(z: Complex64) -> f64 { + if !z.re.is_finite() || !z.im.is_finite() { + // C99 rules: if either part is infinite, return infinity, + // even if the other part is NaN. + if z.re.is_infinite() { + return m::fabs(z.re); + } + if z.im.is_infinite() { + return m::fabs(z.im); + } + return f64::NAN; + } + m::hypot(z.re, z.im) +} + +#[inline] +fn c_abs_checked(z: Complex64) -> Result { + if !z.re.is_finite() || !z.im.is_finite() { + return Ok(c_abs_raw(z)); + } + crate::err::set_errno(0); + let r = m::hypot(z.re, z.im); + if r.is_infinite() { + Err(Error::ERANGE) + } else { + Ok(r) + } +} + /// Convert z to polar coordinates (r, phi). #[inline] pub fn polar(z: Complex64) -> Result<(f64, f64)> { let phi = m::atan2(z.im, z.re); - let r = m::hypot(z.re, z.im); - if r.is_infinite() && z.re.is_finite() && z.im.is_finite() { - return Err(Error::ERANGE); - } + let r = c_abs_checked(z)?; Ok((r, phi)) } @@ -94,7 +121,7 @@ pub fn isinf(z: Complex64) -> bool { /// Complex absolute value (magnitude). #[inline] pub fn abs(z: Complex64) -> f64 { - m::hypot(z.re, z.im) + c_abs_raw(z) } /// Determine whether two complex numbers are close in value. diff --git a/src/math/integer.rs b/src/math/integer.rs index 3c11080..2ef9331 100644 --- a/src/math/integer.rs +++ b/src/math/integer.rs @@ -715,7 +715,7 @@ mod tests { use pyo3::prelude::*; /// Edge i64 values for testing integer math functions (gcd, lcm, isqrt, factorial, comb, perm) - const EDGE_I64: [i64; 44] = [ + const EDGE_I64: &[i64] = &[ // Zero and small values 0, 1, @@ -982,9 +982,9 @@ mod tests { #[test] fn edgetest_gcd() { // Test all edge values - gcd handles arbitrary large integers - for &a in &EDGE_I64 { + for &a in EDGE_I64 { test_gcd_impl(&[a]); - for &b in &EDGE_I64 { + for &b in EDGE_I64 { test_gcd_impl(&[a, b]); } } @@ -997,9 +997,9 @@ mod tests { #[test] fn edgetest_lcm() { // Test all edge values - lcm handles arbitrary large integers - for &a in &EDGE_I64 { + for &a in EDGE_I64 { test_lcm_impl(&[a]); - for &b in &EDGE_I64 { + for &b in EDGE_I64 { test_lcm_impl(&[a, b]); } } @@ -1011,7 +1011,7 @@ mod tests { #[test] fn edgetest_isqrt() { - for &n in &EDGE_I64 { + for &n in EDGE_I64 { test_isqrt_impl(n); } // Additional boundary cases @@ -1031,7 +1031,7 @@ mod tests { #[test] fn edgetest_factorial() { - for &n in &EDGE_I64 { + for &n in EDGE_I64 { // factorial only makes sense for reasonable n values if n >= -10 && n <= 170 { test_factorial_impl(n);