diff --git a/base/src/functions/engineering/mod.rs b/base/src/functions/engineering/mod.rs index 5a125ba..1f549cf 100644 --- a/base/src/functions/engineering/mod.rs +++ b/base/src/functions/engineering/mod.rs @@ -5,5 +5,3 @@ mod convert; mod misc; mod number_basis; mod transcendental; - -pub use transcendental::{fact, fact_double}; diff --git a/base/src/functions/engineering/transcendental/gamma.rs b/base/src/functions/engineering/transcendental/gamma.rs deleted file mode 100644 index cffec59..0000000 --- a/base/src/functions/engineering/transcendental/gamma.rs +++ /dev/null @@ -1,106 +0,0 @@ -use std::f64::consts::PI; - -// FIXME: Please do a better approximation - -/// Lanczos approximation for Gamma(z). -/// Valid for z > 0. For z <= 0 use `gamma` below (which does reflection). -fn gamma_lanczos(z: f64) -> f64 { - const COEFFS: [f64; 9] = [ - 0.999_999_999_999_809_9, - 676.5203681218851, - -1259.1392167224028, - 771.323_428_777_653_1, - -176.615_029_162_140_6, - 12.507343278686905, - -0.13857109526572012, - 9.984_369_578_019_572e-6, - 1.5056327351493116e-7, - ]; - - let g = 7.0; - let mut x = COEFFS[0]; - - // sum_{i=1}^{n} c_i / (z + i) - for (i, c) in COEFFS.iter().enumerate().skip(1) { - x += c / (z + i as f64); - } - - let t = z + g + 0.5; - // √(2π) * t^(z+0.5) * e^{-t} * x - (2.0 * PI).sqrt() * t.powf(z + 0.5) * (-t).exp() * x -} - -/// Gamma function for all real x (except poles), -/// using Lanczos + reflection formula. -fn gamma(x: f64) -> f64 { - if x.is_sign_negative() { - // Reflection formula: Γ(x) = π / (sin(πx) * Γ(1 - x)) - // Beware of sin(πx) near zeros. - let sin_pix = (PI * x).sin(); - if sin_pix == 0.0 { - return f64::NAN; - } - PI / (sin_pix * gamma_lanczos(1.0 - x)) - } else if x == 0.0 { - f64::INFINITY - } else { - gamma_lanczos(x) - } -} - -/// Factorial for real x: x! = Γ(x+1) -pub fn fact(x: f64) -> f64 { - if x < 0.0 { - return f64::NAN; - } - gamma(x + 1.0) -} - -/// Double factorial for real x. -/// -/// Strategy: -/// 1. If x is (almost) integer and >= 0 → do exact product (stable, no surprises) -/// 2. Else: -/// - If x is “even-ish”: x = 2t → x!! = 2^t * Γ(t+1) -/// - If x is “odd-ish”: x = 2t-1 → x!! = Γ(2t+1) / (2^t * Γ(t+1)) -pub fn fact_double(x: f64) -> f64 { - if x < 0.0 { - return f64::NAN; - } - - let xi = x.round(); - let is_int = (x - xi).abs() < 1e-9; - - if is_int { - // integer path - let mut acc = 1.0; - let mut k = xi; - while k > 1.0 { - acc *= k; - k -= 2.0; - } - return acc; - } - - // non-integer path: use continuous extension - let frac = x % 2.0; - // even-ish if close to 0 mod 2 - if frac.abs() < 1e-6 || (frac - 2.0).abs() < 1e-6 { - // x = 2t - let t = x / 2.0; - 2f64.powf(t) * gamma(t + 1.0) - } else if (frac - 1.0).abs() < 1e-6 || (frac + 1.0).abs() < 1e-6 { - // x = 2t - 1 => t = (x+1)/2 - let t = (x + 1.0) / 2.0; - // (2t - 1)!! = Γ(2t + 1) / (2^t Γ(t+1)) - // but 2t - 1 = x, so 2t + 1 = x + 2 - let num = gamma(x + 2.0); - let den = 2f64.powf(t) * gamma(t + 1.0); - num / den - } else { - // x is neither clearly even-ish nor odd-ish – fall back to a generic formula. - // A safe thing to do is to treat it as even-ish: - let t = x / 2.0; - 2f64.powf(t) * gamma(t + 1.0) - } -} diff --git a/base/src/functions/engineering/transcendental/mod.rs b/base/src/functions/engineering/transcendental/mod.rs index d4d301f..9a90a48 100644 --- a/base/src/functions/engineering/transcendental/mod.rs +++ b/base/src/functions/engineering/transcendental/mod.rs @@ -5,7 +5,6 @@ mod bessel_jn_yn; mod bessel_k; mod bessel_util; mod erf; -mod gamma; #[cfg(test)] mod test_bessel; @@ -15,4 +14,3 @@ pub(crate) use bessel_jn_yn::jn as bessel_j; pub(crate) use bessel_jn_yn::yn as bessel_y; pub(crate) use bessel_k::bessel_k; pub(crate) use erf::erf; -pub use gamma::{fact, fact_double}; diff --git a/base/src/functions/mathematical.rs b/base/src/functions/mathematical.rs index e1e522e..549d9c6 100644 --- a/base/src/functions/mathematical.rs +++ b/base/src/functions/mathematical.rs @@ -2,7 +2,6 @@ use crate::cast::NumberOrArray; use crate::constants::{LAST_COLUMN, LAST_ROW}; use crate::expressions::parser::ArrayNode; use crate::expressions::types::CellReferenceIndex; -use crate::functions::engineering::{fact, fact_double}; use crate::number_format::to_precision; use crate::single_number_fn; use crate::{ @@ -456,8 +455,35 @@ impl Model { Ok(1.0 / f64::cosh(f)) }); single_number_fn!(fn_exp, |f: f64| Ok(f64::exp(f))); - single_number_fn!(fn_fact, |f| Ok(fact(f))); - single_number_fn!(fn_factdouble, |f| Ok(fact_double(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 { + 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; + } + Ok(acc) + }); single_number_fn!(fn_sign, |f| Ok(f64::signum(f))); pub(crate) fn fn_pi(&mut self, args: &[Node], cell: CellReferenceIndex) -> CalcResult { diff --git a/base/src/functions/mod.rs b/base/src/functions/mod.rs index 0841fa7..fcc5f9f 100644 --- a/base/src/functions/mod.rs +++ b/base/src/functions/mod.rs @@ -602,6 +602,11 @@ impl Function { "SECH" => Some(Function::Sech), "ACOTH" => Some(Function::Acoth), + "FACT" => Some(Function::Fact), + "FACTDOUBLE" => Some(Function::Factdouble), + "EXP" => Some(Function::Exp), + "SIGN" => Some(Function::Sign), + "PI" => Some(Function::Pi), "ABS" => Some(Function::Abs), "SQRT" => Some(Function::Sqrt), diff --git a/webapp/IronCalc/package-lock.json b/webapp/IronCalc/package-lock.json index 8df1630..f782380 100644 --- a/webapp/IronCalc/package-lock.json +++ b/webapp/IronCalc/package-lock.json @@ -40,7 +40,7 @@ }, "../../bindings/wasm/pkg": { "name": "@ironcalc/wasm", - "version": "0.5.0", + "version": "0.6.0", "license": "MIT/Apache-2.0" }, "node_modules/@adobe/css-tools": { diff --git a/xlsx/tests/calc_tests/FACT_DOUBLEFACT.xlsx b/xlsx/tests/calc_tests/FACT_DOUBLEFACT.xlsx new file mode 100644 index 0000000..04e814f Binary files /dev/null and b/xlsx/tests/calc_tests/FACT_DOUBLEFACT.xlsx differ