Bugfix/nicolas bufixes (#491)
* UPDATE: package lock * FIX: Add function definitions * FIX: Small fix to get FACT working * FIX: We only need integer FACT and FACTDOUBLE * FIX: Make clippy happy
This commit is contained in:
Nicolás Hatcher Andrés
committed by
GitHub
GitHub
parent
7e379e24e7
commit
a768bc5974
@@ -5,5 +5,3 @@ mod convert;
|
||||
mod misc;
|
||||
mod number_basis;
|
||||
mod transcendental;
|
||||
|
||||
pub use transcendental::{fact, fact_double};
|
||||
|
||||
@@ -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)
|
||||
}
|
||||
}
|
||||
@@ -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};
|
||||
|
||||
@@ -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 {
|
||||
|
||||
@@ -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),
|
||||
|
||||
Reference in New Issue
Block a user