From a00ecfdadeb7697bfbef3f3f57d7615bd3153ae0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=A1s=20Hatcher?= Date: Tue, 25 Nov 2025 18:14:17 +0100 Subject: [PATCH] FIX: Uses lookup tables and lngamma for large numbers Fixes #574 --- base/src/functions/mathematical.rs | 149 +++++++++++++++++++++++++---- base/src/test/mod.rs | 1 + base/src/test/test_fn_fact.rs | 29 ++++++ 3 files changed, 159 insertions(+), 20 deletions(-) create mode 100644 base/src/test/test_fn_fact.rs diff --git a/base/src/functions/mathematical.rs b/base/src/functions/mathematical.rs index 2b088af..4fccfb6 100644 --- a/base/src/functions/mathematical.rs +++ b/base/src/functions/mathematical.rs @@ -8,8 +8,74 @@ use crate::single_number_fn; use crate::{ calc_result::CalcResult, expressions::parser::Node, expressions::token::Error, model::Model, }; +use statrs::function::gamma::ln_gamma; +use std::f64::consts::LN_2; use std::f64::consts::PI; +const FACT_TABLE: [f64; 26] = [ + 1.0, // 0! + 1.0, // 1! + 2.0, // 2! + 6.0, // 3! + 24.0, // 4! + 120.0, // 5! + 720.0, // 6! + 5040.0, // 7! + 40320.0, // 8! + 362880.0, // 9! + 3628800.0, // 10! + 39916800.0, // 11! + 479001600.0, // 12! + 6227020800.0, // 13! + 87178291200.0, // 14! + 1307674368000.0, // 15! + 20922789888000.0, // 16! + 355687428096000.0, // 17! + 6402373705728000.0, // 18! + 121645100408832000.0, // 19! + 2432902008176640000.0, // 20! + 51090942171709440000.0, // 21! + 1124000727777607680000.0, // 22! + 25852016738884976640000.0, // 23! + 620448401733239439360000.0, // 24! + 15511210043330985984000000.0, // 25! +]; + +const FACTDOUBLE_TABLE: [f64; 32] = [ + 1.0, + 1.0, + 2.0, + 3.0, + 8.0, + 15.0, + 48.0, + 105.0, + 384.0, + 945.0, + 3840.0, + 10395.0, + 46080.0, + 135135.0, + 645120.0, + 2027025.0, + 10321920.0, + 34459425.0, + 185794560.0, + 654729075.0, + 3715891200.0, + 13749310575.0, + 81749606400.0, + 316234143225.0, + 1961990553600.0, + 7905853580625.0, + 51011754393600.0, + 213458046676875.0, + 1428329123020800.0, + 6190283353629370.0, + 42849873690624000.0, + 191898783962511000.0, +]; + #[cfg(not(target_arch = "wasm32"))] pub fn random() -> f64 { rand::random() @@ -1319,35 +1385,78 @@ impl Model { single_number_fn!(fn_sec, |f| Ok(1.0 / f64::cos(f))); single_number_fn!(fn_sech, |f| Ok(1.0 / f64::cosh(f))); single_number_fn!(fn_exp, |f: f64| Ok(f64::exp(f))); + single_number_fn!(fn_fact, |x: f64| { let x = x.floor(); + if x < 0.0 { return Err(Error::NUM); } - let mut acc = 1.0; - let mut k = 2.0; - while k <= x { - acc *= k; - k += 1.0; - } - Ok(acc) - }); - single_number_fn!(fn_factdouble, |x: f64| { - let x = x.floor(); - if x < -1.0 { - return Err(Error::NUM); - } - if x < 0.0 { + + if x == 0.0 { return Ok(1.0); } - let mut acc = 1.0; - let mut k = if x % 2.0 == 0.0 { 2.0 } else { 1.0 }; - while k <= x { - acc *= k; - k += 2.0; + if x < FACT_TABLE.len() as f64 { + return Ok(FACT_TABLE[x as usize]); } - Ok(acc) + + // Use ln Γ(x+1) to avoid overflow while deciding. + let ln_val = ln_gamma(x + 1.0); + + // If gamma overflows or is invalid, map to NUM + if !ln_val.is_finite() || ln_val > f64::MAX.ln() { + return Err(Error::NUM); + } + + Ok(ln_val.exp()) }); + + single_number_fn!(fn_factdouble, |x: f64| { + let x = x.floor(); + if x <= -1.0 { + return Err(Error::NUM); + } + + if x <= 1.0 { + return Ok(1.0); + } + + // From here x > 1 and integer + let n = x as i64; + + if n < FACTDOUBLE_TABLE.len() as i64 { + return Ok(FACTDOUBLE_TABLE[n as usize]); + } + + // n!! grows very fast, so we compute it using gamma in log-space: + // + // If n is even, n = 2k: + // n!! = 2^k * k! + // If n is odd, n = 2k + 1: + // n!! = (2k+1)! / (2^k * k!) + // + // and we use ln_gamma for factorials. + + let ln_val = if n % 2 == 0 { + // even n = 2k + let k = (n / 2) as f64; + // ln(n!!) = k * ln(2) + ln(k!) + k * LN_2 + ln_gamma(k + 1.0) + } else { + // odd n = 2k + 1 + let k = ((n - 1) / 2) as f64; + let nn = n as f64; + // ln(n!!) = ln((2k+1)!) - (k * ln(2) + ln(k!)) + ln_gamma(nn + 1.0) - (k * LN_2 + ln_gamma(k + 1.0)) + }; + + if !ln_val.is_finite() || ln_val > f64::MAX.ln() { + return Err(Error::NUM); + } + + Ok(ln_val.exp()) + }); + single_number_fn!(fn_sign, |f| { if f == 0.0 { Ok(0.0) diff --git a/base/src/test/mod.rs b/base/src/test/mod.rs index 3e80ace..95816d2 100644 --- a/base/src/test/mod.rs +++ b/base/src/test/mod.rs @@ -18,6 +18,7 @@ mod test_fn_concatenate; mod test_fn_count; mod test_fn_day; mod test_fn_exact; +mod test_fn_fact; mod test_fn_financial; mod test_fn_formulatext; mod test_fn_if; diff --git a/base/src/test/test_fn_fact.rs b/base/src/test/test_fn_fact.rs new file mode 100644 index 0000000..4e0b8f0 --- /dev/null +++ b/base/src/test/test_fn_fact.rs @@ -0,0 +1,29 @@ +#![allow(clippy::unwrap_used)] + +use crate::test::util::new_empty_model; + +#[test] +fn large_numbers() { + let mut model = new_empty_model(); + model._set("A1", "=FACT(170)"); + model._set("A2", "=FACTDOUBLE(36)"); + + model._set("B1", "=FACT(6)"); + model._set("B2", "=FACTDOUBLE(6)"); + + model._set("C3", "=FACTDOUBLE(15)"); + + model._set("F3", "=FACT(-0.1)"); + model._set("F4", "=FACT(0)"); + + model.evaluate(); + + assert_eq!(model._get_text("A1"), *"7.25742E+306"); + assert_eq!(model._get_text("A2"), *"1.67834E+21"); + assert_eq!(model._get_text("B1"), *"720"); + assert_eq!(model._get_text("B2"), *"48"); + assert_eq!(model._get_text("C3"), *"2027025"); + + assert_eq!(model._get_text("F3"), *"#NUM!"); + assert_eq!(model._get_text("F4"), *"1"); +}