From 7e379e24e7a5e2119e686a06265a17b91d87ea9c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=A1s=20Hatcher=20Andr=C3=A9s?= Date: Thu, 30 Oct 2025 18:28:07 +0100 Subject: [PATCH] UPDATE: Adds simple functions (#489) Exp, Fact, Factdouble and sign --- .../src/expressions/parser/static_analysis.rs | 8 ++ base/src/functions/engineering/mod.rs | 2 + .../engineering/transcendental/gamma.rs | 106 ++++++++++++++++++ .../engineering/transcendental/mod.rs | 2 + base/src/functions/mathematical.rs | 5 + base/src/functions/mod.rs | 39 +++---- 6 files changed, 141 insertions(+), 21 deletions(-) create mode 100644 base/src/functions/engineering/transcendental/gamma.rs diff --git a/base/src/expressions/parser/static_analysis.rs b/base/src/expressions/parser/static_analysis.rs index e50e023..50e9d82 100644 --- a/base/src/expressions/parser/static_analysis.rs +++ b/base/src/expressions/parser/static_analysis.rs @@ -842,6 +842,10 @@ fn get_function_args_signature(kind: &Function, arg_count: usize) -> Vec args_signature_scalars(arg_count, 1, 0), Function::Sec => args_signature_scalars(arg_count, 1, 0), Function::Sech => args_signature_scalars(arg_count, 1, 0), + Function::Exp => args_signature_scalars(arg_count, 1, 0), + Function::Fact => args_signature_scalars(arg_count, 1, 0), + Function::Factdouble => args_signature_scalars(arg_count, 1, 0), + Function::Sign => args_signature_scalars(arg_count, 1, 0), } } @@ -1072,5 +1076,9 @@ fn static_analysis_on_function(kind: &Function, args: &[Node]) -> StaticResult { Function::Csch => scalar_arguments(args), Function::Sec => scalar_arguments(args), Function::Sech => scalar_arguments(args), + Function::Exp => scalar_arguments(args), + Function::Fact => scalar_arguments(args), + Function::Factdouble => scalar_arguments(args), + Function::Sign => scalar_arguments(args), } } diff --git a/base/src/functions/engineering/mod.rs b/base/src/functions/engineering/mod.rs index 1f549cf..5a125ba 100644 --- a/base/src/functions/engineering/mod.rs +++ b/base/src/functions/engineering/mod.rs @@ -5,3 +5,5 @@ 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 new file mode 100644 index 0000000..cffec59 --- /dev/null +++ b/base/src/functions/engineering/transcendental/gamma.rs @@ -0,0 +1,106 @@ +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 9a90a48..d4d301f 100644 --- a/base/src/functions/engineering/transcendental/mod.rs +++ b/base/src/functions/engineering/transcendental/mod.rs @@ -5,6 +5,7 @@ mod bessel_jn_yn; mod bessel_k; mod bessel_util; mod erf; +mod gamma; #[cfg(test)] mod test_bessel; @@ -14,3 +15,4 @@ 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 290c0cf..e1e522e 100644 --- a/base/src/functions/mathematical.rs +++ b/base/src/functions/mathematical.rs @@ -2,6 +2,7 @@ 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::{ @@ -454,6 +455,10 @@ impl Model { } else { 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_sign, |f| Ok(f64::signum(f))); pub(crate) fn fn_pi(&mut self, args: &[Node], cell: CellReferenceIndex) -> CalcResult { if !args.is_empty() { diff --git a/base/src/functions/mod.rs b/base/src/functions/mod.rs index 95c7e84..0841fa7 100644 --- a/base/src/functions/mod.rs +++ b/base/src/functions/mod.rs @@ -85,6 +85,11 @@ pub enum Function { Sec, Sech, + Exp, + Fact, + Factdouble, + Sign, + // Information ErrorType, Formulatext, @@ -278,7 +283,7 @@ pub enum Function { } impl Function { - pub fn into_iter() -> IntoIter { + pub fn into_iter() -> IntoIter { [ Function::And, Function::False, @@ -320,6 +325,10 @@ impl Function { Function::Sec, Function::Sech, Function::Power, + Function::Exp, + Function::Fact, + Function::Factdouble, + Function::Sign, Function::Max, Function::Min, Function::Product, @@ -836,7 +845,6 @@ impl fmt::Display for Function { Function::Asinh => write!(f, "ASINH"), Function::Acosh => write!(f, "ACOSH"), Function::Atanh => write!(f, "ATANH"), - Function::Acot => write!(f, "ACOT"), Function::Acoth => write!(f, "ACOTH"), Function::Cot => write!(f, "COT"), @@ -845,7 +853,6 @@ impl fmt::Display for Function { Function::Csch => write!(f, "CSCH"), Function::Sec => write!(f, "SEC"), Function::Sech => write!(f, "SECH"), - Function::Abs => write!(f, "ABS"), Function::Pi => write!(f, "PI"), Function::Sqrt => write!(f, "SQRT"), @@ -910,7 +917,6 @@ impl fmt::Display for Function { Function::Isformula => write!(f, "ISFORMULA"), Function::Type => write!(f, "TYPE"), Function::Sheet => write!(f, "SHEET"), - Function::Average => write!(f, "AVERAGE"), Function::Averagea => write!(f, "AVERAGEA"), Function::Averageif => write!(f, "AVERAGEIF"), @@ -1035,8 +1041,11 @@ impl fmt::Display for Function { Function::Convert => write!(f, "CONVERT"), Function::Delta => write!(f, "DELTA"), Function::Gestep => write!(f, "GESTEP"), - Function::Subtotal => write!(f, "SUBTOTAL"), + Function::Exp => write!(f, "EXP"), + Function::Fact => write!(f, "FACT"), + Function::Factdouble => write!(f, "FACTDOUBLE"), + Function::Sign => write!(f, "SIGN"), } } } @@ -1065,7 +1074,6 @@ impl Model { cell: CellReferenceIndex, ) -> CalcResult { match kind { - // Logical Function::And => self.fn_and(args, cell), Function::False => self.fn_false(args, cell), Function::If => self.fn_if(args, cell), @@ -1077,34 +1085,27 @@ impl Model { Function::Switch => self.fn_switch(args, cell), Function::True => self.fn_true(args, cell), Function::Xor => self.fn_xor(args, cell), - // Math and trigonometry Function::Log => self.fn_log(args, cell), Function::Log10 => self.fn_log10(args, cell), Function::Ln => self.fn_ln(args, cell), Function::Sin => self.fn_sin(args, cell), Function::Cos => self.fn_cos(args, cell), Function::Tan => self.fn_tan(args, cell), - Function::Asin => self.fn_asin(args, cell), Function::Acos => self.fn_acos(args, cell), Function::Atan => self.fn_atan(args, cell), - Function::Sinh => self.fn_sinh(args, cell), Function::Cosh => self.fn_cosh(args, cell), Function::Tanh => self.fn_tanh(args, cell), - Function::Asinh => self.fn_asinh(args, cell), Function::Acosh => self.fn_acosh(args, cell), Function::Atanh => self.fn_atanh(args, cell), - Function::Pi => self.fn_pi(args, cell), Function::Abs => self.fn_abs(args, cell), - Function::Sqrt => self.fn_sqrt(args, cell), Function::Sqrtpi => self.fn_sqrtpi(args, cell), Function::Atan2 => self.fn_atan2(args, cell), Function::Power => self.fn_power(args, cell), - Function::Max => self.fn_max(args, cell), Function::Min => self.fn_min(args, cell), Function::Product => self.fn_product(args, cell), @@ -1116,8 +1117,6 @@ impl Model { Function::Sum => self.fn_sum(args, cell), Function::Sumif => self.fn_sumif(args, cell), Function::Sumifs => self.fn_sumifs(args, cell), - - // Lookup and Reference Function::Choose => self.fn_choose(args, cell), Function::Column => self.fn_column(args, cell), Function::Columns => self.fn_columns(args, cell), @@ -1131,7 +1130,6 @@ impl Model { Function::Rows => self.fn_rows(args, cell), Function::Vlookup => self.fn_vlookup(args, cell), Function::Xlookup => self.fn_xlookup(args, cell), - // Text Function::Concatenate => self.fn_concatenate(args, cell), Function::Exact => self.fn_exact(args, cell), Function::Value => self.fn_value(args, cell), @@ -1149,7 +1147,6 @@ impl Model { Function::Trim => self.fn_trim(args, cell), Function::Unicode => self.fn_unicode(args, cell), Function::Upper => self.fn_upper(args, cell), - // Information Function::Isnumber => self.fn_isnumber(args, cell), Function::Isnontext => self.fn_isnontext(args, cell), Function::Istext => self.fn_istext(args, cell), @@ -1167,7 +1164,6 @@ impl Model { Function::Isformula => self.fn_isformula(args, cell), Function::Type => self.fn_type(args, cell), Function::Sheet => self.fn_sheet(args, cell), - // Statistical Function::Average => self.fn_average(args, cell), Function::Averagea => self.fn_averagea(args, cell), Function::Averageif => self.fn_averageif(args, cell), @@ -1180,7 +1176,6 @@ impl Model { Function::Maxifs => self.fn_maxifs(args, cell), Function::Minifs => self.fn_minifs(args, cell), Function::Geomean => self.fn_geomean(args, cell), - // Date and Time Function::Year => self.fn_year(args, cell), Function::Day => self.fn_day(args, cell), Function::Eomonth => self.fn_eomonth(args, cell), @@ -1206,7 +1201,6 @@ impl Model { Function::WorkdayIntl => self.fn_workday_intl(args, cell), Function::Yearfrac => self.fn_yearfrac(args, cell), Function::Isoweeknum => self.fn_isoweeknum(args, cell), - // Financial Function::Pmt => self.fn_pmt(args, cell), Function::Pv => self.fn_pv(args, cell), Function::Rate => self.fn_rate(args, cell), @@ -1240,7 +1234,6 @@ impl Model { Function::Db => self.fn_db(args, cell), Function::Cumprinc => self.fn_cumprinc(args, cell), Function::Cumipmt => self.fn_cumipmt(args, cell), - // Engineering Function::Besseli => self.fn_besseli(args, cell), Function::Besselj => self.fn_besselj(args, cell), Function::Besselk => self.fn_besselk(args, cell), @@ -1304,6 +1297,10 @@ impl Model { Function::Csch => self.fn_csch(args, cell), Function::Sec => self.fn_sec(args, cell), Function::Sech => self.fn_sech(args, cell), + Function::Exp => self.fn_exp(args, cell), + Function::Fact => self.fn_fact(args, cell), + Function::Factdouble => self.fn_factdouble(args, cell), + Function::Sign => self.fn_sign(args, cell), } } }