Extend the math expression parser's factorial past the whole numbers through the gamma function

libraries/math-parser/src/lib.rs

@@ -148,14 +148,36 @@ mod tests {
 			"gcd(10000000000000000000, 2)",
 			"mod(5, 0)",
 			"(-1)!",
-			"2.5!",
-			"i!",
-			"inf!",
+			"(-2)!",
+			"(-inf)!",
 		] {
 			assert!(evaluate(input).unwrap().is_err(), "expected `{input}` to be an evaluation error");
 		}
 	}
 
+	#[test]
+	fn factorial_extends_through_the_gamma_function() {
+		// Compared relatively, since the reference values span hundreds of orders of magnitude
+		for (input, expected) in [
+			("10.5!", Complex::from(11_899_423.083962247)),
+			("50.25!", Complex::from(8.112744267987253e64)),
+			("170.5!", Complex::from(9.483367566824801e307)),
+			("(-2.5)!", Complex::from(2.3632718012073544)),
+			("(-10.25)!", Complex::from(6.950338494377042e-6)),
+			("(-100.5)!", Complex::from(3.3704592739067173e-157)),
+			("(-170.25)!", Complex::from(2.8837604712815413e-305)),
+			("(-3.5+4i)!", Complex::new(-2.8327740563089983e-5, 5.018195008922803e-5)),
+			("(-20.5+i)!", Complex::new(-5.085633633020186e-19, 7.443957006169107e-20)),
+			("(300i)!", Complex::new(-2.1461927376275517e-204, -9.332698946010592e-204)),
+			("(-1+227i)!", Complex::new(-1.4098280082608946e-157, -2.309687847859193e-156)),
+			("(-1.5+300i)!", Complex::new(-9.760049091627542e-208, 1.5632983579858933e-207)),
+		] {
+			let Value::Number(actual) = evaluate(input).unwrap().unwrap();
+			let actual = actual.as_complex();
+			assert!((actual - expected).norm() / expected.norm() < 1e-12, "`{input}`: expected {expected}, got {actual}");
+		}
+	}
+
 	#[test]
 	fn host_values_are_admitted_where_read() {
 		struct Host;
@@ -432,6 +454,13 @@ mod tests {
 		factorial_nested: "(3 + 2)!" => 120.,
 		factorial_zero: "0!" => 1.,
 		factorial_chain: "3!!" => 720., // (3!)! = 6! = 720
+		factorial_half: "0.5!" => std::f64::consts::PI.sqrt() / 2.,
+		factorial_negative_half: "(-0.5)!" => std::f64::consts::PI.sqrt(),
+		factorial_fractional: "2.5!" => 3.323350970447842,
+		factorial_negative_fractional: "(-1.5)!" => -2. * std::f64::consts::PI.sqrt(),
+		factorial_imaginary: "i!" => Complex::new(0.498015668118356, -0.1549498283018107),
+		factorial_complex_left_half_plane: "(-1.5 + 2i)!" => Complex::new(-0.03903884916211552, -0.03516787606268694),
+		factorial_infinity: "inf!" => f64::INFINITY,
 
 		// Operations with negative values
 		negative_nested_parentheses: "-(5 + 3 * (2 - 1))" => -8.,
@@ -627,6 +656,11 @@ mod tests {
 		// Overflow-safe evaluation
 		factorial_overflows_to_infinity: "171!" => f64::INFINITY,
 		factorial_huge_input: "10000000000000000000000!" => f64::INFINITY,
+		factorial_fractional_overflows_to_infinity: "171.5!" => f64::INFINITY,
+		factorial_huge_fractional_input: "4503599627370495.5!" => f64::INFINITY,
+		factorial_complex_overflows_to_infinity: "|(171 + 0.5i)!|" => f64::INFINITY,
+		factorial_complex_underflows_to_zero: "(-190.5 + 0.5i)!" => 0.,
+		factorial_huge_negative_underflows_to_zero: "(-740.5)!" => 0.,
 		lcm_huge_no_overflow: "lcm(1099511627776, 1099511627775)" => 1099511627776. * 1099511627775.,
 		long_literal: "10000000000000000000000" => 1e22,
 		huge_exponent_saturates: "1e4294967296" => f64::INFINITY,

libraries/math-parser/src/value.rs

@@ -1,4 +1,5 @@
 use crate::ast::{BinaryOp, UnaryOp};
+use std::f64::consts::PI;
 
 pub type Complex = num_complex::Complex<f64>;
 
@@ -235,24 +236,88 @@ impl Number {
 				Number::Real(real) => real.abs(),
 				Number::Complex(complex) => complex.norm(),
 			})),
-			UnaryOp::Fac => {
-				// A factorial is defined for whole numbers at or above zero
-				let Number::Real(real) = self else { return None };
-				let whole = real.round();
-				if !real.is_finite() || whole < 0. || (real - whole).abs() > f64::EPSILON {
-					return None;
-				}
-
-				// Infinity above 170!, which overflows f64, also keeps huge inputs from spinning the loop
-				if whole > 170. {
-					return Some(Number::Real(f64::INFINITY));
-				}
-				Some(Number::Real((1..=whole as u64).fold(1., |accumulated, k| accumulated * k as f64)))
-			}
+			UnaryOp::Fac => Some(match self {
+				Number::Real(real) => Number::Real(real_factorial(real)?),
+				Number::Complex(complex) => Number::Complex(complex_gamma(complex + 1.)),
+			}),
 		}
 	}
 
 	pub fn from_f64(x: f64) -> Self {
 		Self::Real(x)
 	}
 }
+
+/// The factorial of a real number: the exact product for a whole number, and `x! = Γ(x + 1)` past the whole numbers, or
+/// `None` at the negative integers (the gamma function's poles, where no signed limit exists) and at -∞.
+fn real_factorial(x: f64) -> Option<f64> {
+	// The factorial overflows f64 from 171! on, and returning early keeps a huge whole number from spinning the product loop
+	if x > 171. {
+		return Some(f64::INFINITY);
+	}
+
+	if x.fract() == 0. {
+		(x >= 0.).then(|| (1..=x as u64).fold(1., |accumulated, k| accumulated * k as f64))
+	} else {
+		x.is_finite().then(|| real_gamma(x + 1.))
+	}
+}
+
+/// The Lanczos approximation's shift, for which [`LANCZOS_COEFFICIENTS`] give 15 significant digits near 1, falling to 13 by the f64 overflow limit.
+const LANCZOS_G: f64 = 7.;
+const LANCZOS_COEFFICIENTS: [f64; 9] = [
+	0.9999999999998099,
+	676.5203681218851,
+	-1259.1392167224028,
+	771.3234287776531,
+	-176.6150291621406,
+	12.507343278686905,
+	-0.13857109526572012,
+	9.984369578019572e-6,
+	1.5056327351493116e-7,
+];
+
+/// The gamma function on the reals by the Lanczos approximation. Below 1/2, where the series loses accuracy, it reflects
+/// through `Γ(x) Γ(1 - x) = π / sin(πx)`.
+fn real_gamma(x: f64) -> f64 {
+	if x < 0.5 {
+		return PI / (PI * x).sin() / real_gamma(1. - x);
+	}
+
+	let x = x - 1.;
+	let series = (1..LANCZOS_COEFFICIENTS.len()).fold(LANCZOS_COEFFICIENTS[0], |sum, k| sum + LANCZOS_COEFFICIENTS[k] / (x + k as f64));
+	let t = x + LANCZOS_G + 0.5;
+
+	// One exponential for the whole `t^(x + 1/2) e^-t` factor, so past f64's range it is ∞ rather than the NaN of ∞ × 0
+	(2. * PI).sqrt() * ((x + 0.5) * t.ln() - t).exp() * series
+}
+
+/// The gamma function over the complex plane, taken as one exponential of its logarithm so a magnitude past f64's range is ∞ or 0.
+fn complex_gamma(z: Complex) -> Complex {
+	complex_log_gamma(z).exp()
+}
+
+/// The natural logarithm of the gamma function over the complex plane, by the same approximation and reflection as [`real_gamma`].
+fn complex_log_gamma(z: Complex) -> Complex {
+	if z.re < 0.5 {
+		// Shifting the real part into `[0, 1)` keeps the sine exact, with a half turn (a sign) for each odd shift
+		let shift = z.re.floor();
+		let log_sin = complex_log_sin(PI * (z - shift)) + Complex::new(0., if shift.rem_euclid(2.) == 0. { 0. } else { PI });
+		return PI.ln() - log_sin - complex_log_gamma(1. - z);
+	}
+
+	let z = z - 1.;
+	let series = (1..LANCZOS_COEFFICIENTS.len()).fold(Complex::from(LANCZOS_COEFFICIENTS[0]), |sum, k| sum + LANCZOS_COEFFICIENTS[k] / (z + k as f64));
+	let t = z + LANCZOS_G + 0.5;
+
+	0.5 * (2. * PI).ln() + (z + 0.5) * t.ln() - t + series.ln()
+}
+
+/// `ln sin(w)` without forming `sin(w)`, which overflows past an imaginary part of about 710: with `s` the sign of that part,
+/// `sin(w) = e^(-siw) (e^(2siw) - 1) / (2si)`, and the growing exponential becomes a plain shift of the logarithm.
+fn complex_log_sin(w: Complex) -> Complex {
+	let s = if w.im < 0. { -1. } else { 1. };
+	let siw = Complex::new(0., s) * w;
+
+	-siw + (((2. * siw).exp() - 1.) / Complex::new(0., 2. * s)).ln()
+}