From 5cec1adf3d5a01e05dbaf7466b895550cc6e4e46 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Fri, 31 Oct 2025 20:04:30 +0100
Subject: [PATCH 01/35] Implement the Gamma function by Toshio Fukushima that
 has very low relative error.

This commit also extends that function to handle more special cases correctly.
---
 scripts/gamma_table.py |   5 +-
 src/gamma.rs           | 108 ++++++++++++++++++++++++++++++++++++-----
 src/utils.rs           |  12 +++++
 tests/gamma_test.rs    |  52 +++++++++++---------
 4 files changed, 139 insertions(+), 38 deletions(-)

diff --git a/scripts/gamma_table.py b/scripts/gamma_table.py
index 3ce4fb3..a9db97b 100644
--- a/scripts/gamma_table.py
+++ b/scripts/gamma_table.py
@@ -29,7 +29,8 @@
 # └──────────────────────────────────────────────────────────┘
 # x values to test (excluding inputs that would result in complex outputs)
 gamma_x_values = [
-    -1.5, -0.5,  # Negative non-integer values that result in real outputs
+    -1.5, -0.5, -100.5,  # Negative non-integer values that result in real outputs
+    -4503599627370495.5,  # Large negative non-integer value that underflows to zero.
     0.1, 0.2, 0.5,  # Values between 0 and 1
     1.0, 1.5, 2.0, 2.5,  # Values between 1 and 3
     3.0, 4.0, 5.0,  # Small integers and values between them
@@ -45,7 +46,7 @@
 # Print in Rust code format
 print("const GAMMA_TABLE: [(f64, f64); {}] = [".format(len(table)))
 for x, y in table:
-    print(f"    ({x:.14e}, {y:.14e}),")
+    print(f"    ({x:.16e}, {y:.16e}),")
 print("];")
 
 # a values to test for gammp, gammq, and invgammp
diff --git a/src/gamma.rs b/src/gamma.rs
index 7d8ad28..f8a5fce 100644
--- a/src/gamma.rs
+++ b/src/gamma.rs
@@ -7,8 +7,8 @@
 //! - `gammq`: Calculates the regularized upper incomplete gamma function.
 //! - `invgammp`: Calculates the inverse of the regularized lower incomplete gamma function.
 
-use crate::utils::factorial;
-use crate::{EPS, FPMIN, W, Y};
+use crate::{utils::polynomial, EPS, FPMIN, W, Y};
+use core::f64;
 use std::f64::consts::PI;
 const ASWITCH: usize = 100;
 const NGAU: usize = 18;
@@ -54,19 +54,103 @@ pub fn ln_gamma(z: f64) -> f64 {
 /// # Returns
 ///
 /// The value of the gamma function at `z`
-pub fn gamma(z: f64) -> f64 {
-    if z >= 1f64 {
-        let z_int = z as usize;
-        if (z - (z_int as f64)).abs() < EPS {
-            return factorial(z_int - 1);
-        }
+///
+// Based on the Fortran implementation by Toshio Fukushima: <https://www.researchgate.net/publication/336578125_xgamtxt_Fortran_90_test_program_package_of_qgam_dgam_and_sgam_fast_computation_of_Gamma_function_with_quadruple_double_and_single_precision_accuracy_respectively>
+pub fn gamma(mut x: f64) -> f64 {
+    /// This input gives an output of [`f64::MAX`].
+    /// Computed with WolframAlpha: https://www.wolframalpha.com/input?i=solve+gamma%28x%29+%3D+%281+%E2%88%92+2%5E%28%E2%88%9253%29%29*2%5E1024+for+x+%3E+1.
+    const MAX_INPUT: f64 = 171.624376956302720;
+
+    // Special cases.
+    if x > MAX_INPUT {
+        // Output is too large to represent.
+        return f64::INFINITY;
+    } else if x == 0.0 {
+        // The gamma function diverges for an input of zero.
+        // It diverges to positive or negative infinity depending on which direction zero is approached from.
+        // According to entry F.9.5.4 of the ISO C 99 standard (https://www.open-std.org/jtc1/sc22/wg14/www/docs/n1256.pdf)
+        // gamma(+/- 0) = +/- infinity.
+        // This is the standard followed by `scipy.special.gamma` as well as the `tgamma` function in `libc`.
+        // Since the Gamma function in (nightly) Rust's standard library currently corresponds to `tgamma`
+        // and thus follows the standard we do the same here.
+        return f64::INFINITY.copysign(x);
+    } else if x <= 0.0 && x.fract() == 0.0 {
+        // The Gamma function diverges for non-positive integers.
+        // There is however no clear way for the caller to signal whether the input approached the integer from above or below.
+        // According to the same standard as above the gamma function should return NaN for these inputs.
+        return f64::NAN;
     }
 
-    if z < 0.5 {
-        PI / ((PI * z).sin() * gamma(1f64 - z))
+    let f = if x > 3.5 {
+        let mut f = 1.0;
+        // Decrement x by 1 until less than 3.5.
+        // This will not be a massively long loop since we checked for too large input above.
+        while x >= 3.5 {
+            x -= 1.0;
+            f *= x;
+        }
+        f
+    } else if x < 2.5 {
+        let mut f = 1.0;
+        // Increment x by 1 until larger than 2.5.
+        while x <= 2.5 {
+            f *= x;
+            if f.is_infinite() {
+                // There is no point in continuing the calculation,
+                // as `f` will remain infinite when multiplied by any finite number.
+                // The final result will thus be `g` / inf,
+                // which underflows to 0 since `g` (defined below) is always finite.
+                // This early return stops the loop from running too long for very large negative inputs.
+                return 0.0;
+            }
+            x += 1.0;
+        }
+        f.recip()
+    } else {
+        1.0
+    };
+
+    let g = if x > 3.0 {
+        polynomial(
+            x - 3.25,
+            &[
+                2.5492569667185291891,
+                2.5925711651299805640,
+                1.7769198047099783654,
+                0.85885382288064637892,
+                0.34626855336056872905,
+                0.11526048279921955157,
+                0.034350198197417378854,
+                0.0088826032863543439287,
+                0.0021498495721144269788,
+                0.00045576674678513314905,
+                0.000095940184093793519324,
+                0.000016738398919923317512,
+            ],
+        )
+    } else if x < 3.0 {
+        polynomial(
+            x - 2.75,
+            &[
+                1.6083594219855455740,
+                1.3170871791928586304,
+                0.89116794802647661936,
+                0.38477912014818524560,
+                0.15595585685363435567,
+                0.046735160424814165697,
+                0.014576803693379529708,
+                0.0032284372223208491985,
+                0.00091425286320138797078,
+                0.00013276357775573942004,
+                0.000048420884483658205918,
+                7.3776577743984342365e-7,
+            ],
+        )
     } else {
-        ln_gamma(z).exp()
-    }
+        2.0
+    };
+
+    g * f
 }
 
 // =============================================================================
diff --git a/src/utils.rs b/src/utils.rs
index a67fe34..a5cf61b 100644
--- a/src/utils.rs
+++ b/src/utils.rs
@@ -103,3 +103,15 @@ pub fn sign(a: f64, b: f64) -> f64 {
         -(a.abs())
     }
 }
+
+/// Evaluate a polynomial at x using [Horner's method](https://en.wikipedia.org/wiki/Horner%27s_method).
+///
+/// The coefficients are assumed to be sorted in increasing order by the degree of the associated x-term.
+///
+/// For example, to evaluate 4x^2 + 3x + 2, use `coefficients = [2.0, 3.0, 4.0]`.
+pub fn polynomial<const N: usize>(x: f64, coefficients: &[f64; N]) -> f64 {
+    coefficients
+        .iter()
+        .rev()
+        .fold(0.0, |acc, &coeff| acc * x + coeff)
+}
diff --git a/tests/gamma_test.rs b/tests/gamma_test.rs
index a3f37e3..983e137 100644
--- a/tests/gamma_test.rs
+++ b/tests/gamma_test.rs
@@ -22,9 +22,11 @@ fn test_gamma() {
     for (x, y) in GAMMA_TABLE {
         let result = gamma(x);
         let abs_eps = f64::EPSILON;
-        let rel_eps = 1e-10;
+        let rel_eps = 1e-14;
         assert_relative_eq!(result, y, epsilon = abs_eps, max_relative = rel_eps);
     }
+
+    assert_eq!(gamma(-4503599627370495.5), 0.0);
 }
 
 #[test]
@@ -63,7 +65,7 @@ proptest! {
         let result = gamma(x);
         let expected = unsafe { tgamma(x) };
         let abs_eps = f64::EPSILON;
-        let rel_eps = 1e-9;
+        let rel_eps = 1e-14;
         assert_relative_eq!(result, expected, epsilon = abs_eps, max_relative = rel_eps);
     }
 
@@ -102,28 +104,30 @@ const LN_GAMMA_TABLE: [(f64, f64); 19] = [
     (1.00000000000000e+10, 2.20258509288811e+11),
 ];
 
-const GAMMA_TABLE: [(f64, f64); 21] = [
-    (-1.50000000000000e+00, 2.36327180120735e+00),
-    (-5.00000000000000e-01, -3.54490770181103e+00),
-    (1.00000000000000e-01, 9.51350769866873e+00),
-    (2.00000000000000e-01, 4.59084371199880e+00),
-    (5.00000000000000e-01, 1.77245385090552e+00),
-    (1.00000000000000e+00, 1.00000000000000e+00),
-    (1.50000000000000e+00, 8.86226925452758e-01),
-    (2.00000000000000e+00, 1.00000000000000e+00),
-    (2.50000000000000e+00, 1.32934038817914e+00),
-    (3.00000000000000e+00, 2.00000000000000e+00),
-    (4.00000000000000e+00, 6.00000000000000e+00),
-    (5.00000000000000e+00, 2.40000000000000e+01),
-    (1.00000000000000e+01, 3.62880000000000e+05),
-    (2.00000000000000e+01, 1.21645100408832e+17),
-    (5.00000000000000e+01, 6.08281864034268e+62),
-    (2.50000000000000e-01, 3.62560990822191e+00),
-    (7.50000000000000e-01, 1.22541670246518e+00),
-    (1.00000000000000e-05, 9.99994227942255e+04),
-    (1.00000000000000e-10, 9.99999999942278e+09),
-    (1.00000000000000e+02, 9.33262154439442e+155),
-    (1.70000000000000e+02, 4.26906800900471e+304),
+const GAMMA_TABLE: [(f64, f64); 23] = [
+    (-1.5000000000000000e+00, 2.3632718012073548e+00),
+    (-5.0000000000000000e-01, -3.5449077018110318e+00),
+    (-1.0050000000000000e+02, -3.3536908198076757e-159),
+    (-4.5035996273704955e+15, -0.0000000000000000e+00),
+    (1.0000000000000001e-01, 9.5135076986687324e+00),
+    (2.0000000000000001e-01, 4.5908437119988035e+00),
+    (5.0000000000000000e-01, 1.7724538509055159e+00),
+    (1.0000000000000000e+00, 1.0000000000000000e+00),
+    (1.5000000000000000e+00, 8.8622692545275794e-01),
+    (2.0000000000000000e+00, 1.0000000000000000e+00),
+    (2.5000000000000000e+00, 1.3293403881791370e+00),
+    (3.0000000000000000e+00, 2.0000000000000000e+00),
+    (4.0000000000000000e+00, 6.0000000000000000e+00),
+    (5.0000000000000000e+00, 2.4000000000000000e+01),
+    (1.0000000000000000e+01, 3.6288000000000000e+05),
+    (2.0000000000000000e+01, 1.2164510040883200e+17),
+    (5.0000000000000000e+01, 6.0828186403426756e+62),
+    (2.5000000000000000e-01, 3.6256099082219082e+00),
+    (7.5000000000000000e-01, 1.2254167024651774e+00),
+    (1.0000000000000001e-05, 9.9999422794225538e+04),
+    (1.0000000000000000e-10, 9.9999999994227848e+09),
+    (1.0000000000000000e+02, 9.3326215443944153e+155),
+    (1.7000000000000000e+02, 4.2690680090047053e+304),
 ];
 
 const GAMMP_TABLE: [(f64, f64, f64); 42] = [

From 2a73150f513794e95afbbd7540c6ff9c94f2a891 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Fri, 31 Oct 2025 20:37:44 +0100
Subject: [PATCH 02/35] Check strictly < 0 since we have already checked that x
 is not 0

---
 src/gamma.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/gamma.rs b/src/gamma.rs
index eb7314f..10c0c98 100644
--- a/src/gamma.rs
+++ b/src/gamma.rs
@@ -74,7 +74,7 @@ pub fn gamma(mut x: f64) -> f64 {
         // Since the Gamma function in (nightly) Rust's standard library currently corresponds to `tgamma`
         // and thus follows the standard we do the same here.
         return f64::INFINITY.copysign(x);
-    } else if x <= 0.0 && x.fract() == 0.0 {
+    } else if x < 0.0 && x.fract() == 0.0 {
         // The Gamma function diverges for non-positive integers.
         // There is however no clear way for the caller to signal whether the input approached the integer from above or below.
         // According to the same standard as above the gamma function should return NaN for these inputs.

From 77af3dd6610bc0bedf1c45b824234f9a47776354 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Fri, 31 Oct 2025 20:37:51 +0100
Subject: [PATCH 03/35] Make the tests stricter

---
 tests/gamma_test.rs | 39 ++++++++++++++++++++++-----------------
 1 file changed, 22 insertions(+), 17 deletions(-)

diff --git a/tests/gamma_test.rs b/tests/gamma_test.rs
index b41df38..e903b94 100644
--- a/tests/gamma_test.rs
+++ b/tests/gamma_test.rs
@@ -66,44 +66,49 @@ fn test_invgammp() {
 #[test]
 fn test_gamma_edge_cases() {
     // Test gamma at special values
+
+    // Γ(±0) = ±∞
+    assert_eq!(gamma(0.0), f64::INFINITY);
+    assert_eq!(gamma(-0.0), -f64::INFINITY);
+
     // Γ(1) = 1
-    assert_relative_eq!(gamma(1.0), 1.0, epsilon = 1e-15);
+    assert_eq!(gamma(1.0), 1.0);
 
     // Γ(2) = 1
-    assert_relative_eq!(gamma(2.0), 1.0, epsilon = 1e-15);
+    assert_eq!(gamma(2.0), 1.0);
 
     // Γ(1/2) = √π
-    let sqrt_pi = std::f64::consts::PI.sqrt();
-    assert_relative_eq!(gamma(0.5), sqrt_pi, epsilon = 1e-10);
+    let sqrt_pi = core::f64::consts::PI.sqrt();
+    assert_relative_eq!(gamma(0.5), sqrt_pi);
 
     // Γ(3/2) = √π/2
-    assert_relative_eq!(gamma(1.5), sqrt_pi / 2.0, epsilon = 1e-10);
+    assert_relative_eq!(gamma(1.5), sqrt_pi / 2.0);
 
     // Γ(5/2) = 3√π/4
-    assert_relative_eq!(gamma(2.5), 3.0 * sqrt_pi / 4.0, epsilon = 1e-10);
+    assert_relative_eq!(gamma(2.5), 3.0 * sqrt_pi / 4.0);
 
-    // Test very small positive values
-    let result_small = gamma(1e-8);
-    assert!(result_small > 1e7 && result_small < 1e9);
+    // Test very small positive values.
+    // Computed with WolframAlpha: https://www.wolframalpha.com/input?i=Gamma%281e-8%29.
+    assert_relative_eq!(gamma(1e-8), 99999999.42278434, max_relative = 1.5*f64::EPSILON);
 
-    // Test negative values (using reflection formula)
+    // Test negative values
     // Γ(-0.5) = -2√π
-    assert_relative_eq!(gamma(-0.5), -2.0 * sqrt_pi, epsilon = 1e-10);
+    assert_relative_eq!(gamma(-0.5), -2.0 * sqrt_pi);
 
     // Test near-integer negative values
-    let result_neg = gamma(-1.5);
-    assert!(result_neg > 2.0 && result_neg < 3.0);
+    // Computed with WolframAlpha: https://www.wolframalpha.com/input?i=Gamma%28-1.5%29
+    assert_relative_eq!(gamma(-1.5), 2.363271801207355);
 
     // Test value very close to 0 but positive
-    let result_tiny = gamma(1e-15);
-    assert!(result_tiny > 1e14);
+    // Computed with WolframAlpha: https://www.wolframalpha.com/input?i=Gamma%281e-15%29
+    assert_relative_eq!(gamma(1e-15), 999999999999999.4, max_relative = 2.0*f64::EPSILON);
 
     // Test integer factorials
     // Γ(6) = 5! = 120
-    assert_relative_eq!(gamma(6.0), 120.0, epsilon = 1e-13);
+    assert_eq!(gamma(6.0), 120.0);
 
     // Γ(11) = 10! = 3628800
-    assert_relative_eq!(gamma(11.0), 3628800.0, epsilon = 1e-9);
+    assert_eq!(gamma(11.0), 3628800.0);
 }
 
 #[test]

From d035b0fae9867964ff935d54bdca9c2aac1d6ddb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Fri, 31 Oct 2025 20:41:08 +0100
Subject: [PATCH 04/35] Remove unnecessary line in doc-comment

---
 src/gamma.rs | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/gamma.rs b/src/gamma.rs
index 10c0c98..6216adf 100644
--- a/src/gamma.rs
+++ b/src/gamma.rs
@@ -54,7 +54,6 @@ pub fn ln_gamma(z: f64) -> f64 {
 /// # Returns
 ///
 /// The value of the gamma function at `z`
-///
 // Based on the Fortran implementation by Toshio Fukushima: <https://www.researchgate.net/publication/336578125_xgamtxt_Fortran_90_test_program_package_of_qgam_dgam_and_sgam_fast_computation_of_Gamma_function_with_quadruple_double_and_single_precision_accuracy_respectively>
 pub fn gamma(mut x: f64) -> f64 {
     /// This input gives an output of [`f64::MAX`].

From 5ff810c07c060fde75ec8ffb51e3df53d887fb9b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Fri, 31 Oct 2025 20:41:13 +0100
Subject: [PATCH 05/35] `cargo fmt`

---
 tests/gamma_test.rs | 12 ++++++++++--
 1 file changed, 10 insertions(+), 2 deletions(-)

diff --git a/tests/gamma_test.rs b/tests/gamma_test.rs
index e903b94..919bb91 100644
--- a/tests/gamma_test.rs
+++ b/tests/gamma_test.rs
@@ -89,7 +89,11 @@ fn test_gamma_edge_cases() {
 
     // Test very small positive values.
     // Computed with WolframAlpha: https://www.wolframalpha.com/input?i=Gamma%281e-8%29.
-    assert_relative_eq!(gamma(1e-8), 99999999.42278434, max_relative = 1.5*f64::EPSILON);
+    assert_relative_eq!(
+        gamma(1e-8),
+        99999999.42278434,
+        max_relative = 1.5 * f64::EPSILON
+    );
 
     // Test negative values
     // Γ(-0.5) = -2√π
@@ -101,7 +105,11 @@ fn test_gamma_edge_cases() {
 
     // Test value very close to 0 but positive
     // Computed with WolframAlpha: https://www.wolframalpha.com/input?i=Gamma%281e-15%29
-    assert_relative_eq!(gamma(1e-15), 999999999999999.4, max_relative = 2.0*f64::EPSILON);
+    assert_relative_eq!(
+        gamma(1e-15),
+        999999999999999.4,
+        max_relative = 2.0 * f64::EPSILON
+    );
 
     // Test integer factorials
     // Γ(6) = 5! = 120

From 1f5b3848e0dfb4945fde1628994b10640f21b546 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Fri, 31 Oct 2025 20:45:10 +0100
Subject: [PATCH 06/35] Explain handling of x -negative int

---
 src/gamma.rs | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/gamma.rs b/src/gamma.rs
index 6216adf..b310a4c 100644
--- a/src/gamma.rs
+++ b/src/gamma.rs
@@ -97,6 +97,7 @@ pub fn gamma(mut x: f64) -> f64 {
             if f.is_infinite() {
                 // There is no point in continuing the calculation,
                 // as `f` will remain infinite when multiplied by any finite number.
+                // (except zero, but that can only happen if x is an integer, which we have handled above.)
                 // The final result will thus be `g` / inf,
                 // which underflows to 0 since `g` (defined below) is always finite.
                 // This early return stops the loop from running too long for very large negative inputs.

From 45246486c270920ff6cc160008d64ec1f966866e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Fri, 31 Oct 2025 20:49:32 +0100
Subject: [PATCH 07/35] Fix clippy lints

---
 src/gamma.rs | 50 +++++++++++++++++++++++++-------------------------
 1 file changed, 25 insertions(+), 25 deletions(-)

diff --git a/src/gamma.rs b/src/gamma.rs
index b310a4c..897210d 100644
--- a/src/gamma.rs
+++ b/src/gamma.rs
@@ -58,7 +58,7 @@ pub fn ln_gamma(z: f64) -> f64 {
 pub fn gamma(mut x: f64) -> f64 {
     /// This input gives an output of [`f64::MAX`].
     /// Computed with WolframAlpha: https://www.wolframalpha.com/input?i=solve+gamma%28x%29+%3D+%281+%E2%88%92+2%5E%28%E2%88%9253%29%29*2%5E1024+for+x+%3E+1.
-    const MAX_INPUT: f64 = 171.624376956302720;
+    const MAX_INPUT: f64 = 171.624_376_956_302_7;
 
     // Special cases.
     if x > MAX_INPUT {
@@ -114,36 +114,36 @@ pub fn gamma(mut x: f64) -> f64 {
         polynomial(
             x - 3.25,
             &[
-                2.5492569667185291891,
-                2.5925711651299805640,
-                1.7769198047099783654,
-                0.85885382288064637892,
-                0.34626855336056872905,
-                0.11526048279921955157,
-                0.034350198197417378854,
-                0.0088826032863543439287,
-                0.0021498495721144269788,
-                0.00045576674678513314905,
-                0.000095940184093793519324,
-                0.000016738398919923317512,
+                2.549_256_966_718_529,
+                2.592_571_165_129_980_8,
+                1.776_919_804_709_978_3,
+                0.858_853_822_880_646_4,
+                0.346_268_553_360_568_7,
+                0.115_260_482_799_219_55,
+                0.034_350_198_197_417_38,
+                0.008_882_603_286_354_344,
+                0.002_149_849_572_114_427,
+                0.000_455_766_746_785_133_16,
+                0.000_095_940_184_093_793_52,
+                0.000_016_738_398_919_923_317,
             ],
         )
     } else if x < 3.0 {
         polynomial(
             x - 2.75,
             &[
-                1.6083594219855455740,
-                1.3170871791928586304,
-                0.89116794802647661936,
-                0.38477912014818524560,
-                0.15595585685363435567,
-                0.046735160424814165697,
-                0.014576803693379529708,
-                0.0032284372223208491985,
-                0.00091425286320138797078,
-                0.00013276357775573942004,
-                0.000048420884483658205918,
-                7.3776577743984342365e-7,
+                1.608_359_421_985_545_5,
+                1.317_087_179_192_858_7,
+                0.891_167_948_026_476_6,
+                0.384_779_120_148_185_27,
+                0.155_955_856_853_634_36,
+                0.046_735_160_424_814_17,
+                0.014_576_803_693_379_53,
+                0.003_228_437_222_320_849,
+                0.000_914_252_863_201_387_9,
+                0.000_132_763_577_755_739_43,
+                0.000_048_420_884_483_658_204,
+                7.377_657_774_398_435e-7,
             ],
         )
     } else {

From de8c98b7508de06f5af7282ea404571b4707cbdf Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Fri, 31 Oct 2025 21:49:23 +0100
Subject: [PATCH 08/35] Set the expected answer in the test og Gamma(-1.5) to
 the exact correct answer

---
 tests/gamma_test.rs | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/tests/gamma_test.rs b/tests/gamma_test.rs
index 919bb91..9bde628 100644
--- a/tests/gamma_test.rs
+++ b/tests/gamma_test.rs
@@ -100,8 +100,7 @@ fn test_gamma_edge_cases() {
     assert_relative_eq!(gamma(-0.5), -2.0 * sqrt_pi);
 
     // Test near-integer negative values
-    // Computed with WolframAlpha: https://www.wolframalpha.com/input?i=Gamma%28-1.5%29
-    assert_relative_eq!(gamma(-1.5), 2.363271801207355);
+    assert_relative_eq!(gamma(-1.5), 4.0*sqrt_pi/3.0);
 
     // Test value very close to 0 but positive
     // Computed with WolframAlpha: https://www.wolframalpha.com/input?i=Gamma%281e-15%29

From 1937fc851a724db14254c749b2f41bdbd20a4232 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Fri, 31 Oct 2025 21:51:27 +0100
Subject: [PATCH 09/35] cargo fmt

---
 tests/gamma_test.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/gamma_test.rs b/tests/gamma_test.rs
index 9bde628..2fb6cd5 100644
--- a/tests/gamma_test.rs
+++ b/tests/gamma_test.rs
@@ -100,7 +100,7 @@ fn test_gamma_edge_cases() {
     assert_relative_eq!(gamma(-0.5), -2.0 * sqrt_pi);
 
     // Test near-integer negative values
-    assert_relative_eq!(gamma(-1.5), 4.0*sqrt_pi/3.0);
+    assert_relative_eq!(gamma(-1.5), 4.0 * sqrt_pi / 3.0);
 
     // Test value very close to 0 but positive
     // Computed with WolframAlpha: https://www.wolframalpha.com/input?i=Gamma%281e-15%29

From bc244d8b76dc8fc256c750ab765af8f4d4dfef9e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Fri, 31 Oct 2025 23:31:05 +0100
Subject: [PATCH 10/35] Implement the `polynomial` function in terms of
 `IntoIterator` instead of a const N array.

According to godbolt.org both of these implementations result in the same assembly, but since this one is more flexible for the caller we should use this one.
---
 src/gamma.rs | 48 ++++++++++++++++++++++++------------------------
 src/utils.rs | 14 ++++++++------
 2 files changed, 32 insertions(+), 30 deletions(-)

diff --git a/src/gamma.rs b/src/gamma.rs
index 897210d..8bcb445 100644
--- a/src/gamma.rs
+++ b/src/gamma.rs
@@ -113,37 +113,37 @@ pub fn gamma(mut x: f64) -> f64 {
     let g = if x > 3.0 {
         polynomial(
             x - 3.25,
-            &[
-                2.549_256_966_718_529,
-                2.592_571_165_129_980_8,
-                1.776_919_804_709_978_3,
-                0.858_853_822_880_646_4,
-                0.346_268_553_360_568_7,
-                0.115_260_482_799_219_55,
-                0.034_350_198_197_417_38,
-                0.008_882_603_286_354_344,
-                0.002_149_849_572_114_427,
-                0.000_455_766_746_785_133_16,
-                0.000_095_940_184_093_793_52,
+            [
                 0.000_016_738_398_919_923_317,
+                0.000_095_940_184_093_793_52,
+                0.000_455_766_746_785_133_16,
+                0.002_149_849_572_114_427,
+                0.008_882_603_286_354_344,
+                0.034_350_198_197_417_38,
+                0.115_260_482_799_219_55,
+                0.346_268_553_360_568_7,
+                0.858_853_822_880_646_4,
+                1.776_919_804_709_978_3,
+                2.592_571_165_129_980_8,
+                2.549_256_966_718_529,
             ],
         )
     } else if x < 3.0 {
         polynomial(
             x - 2.75,
-            &[
-                1.608_359_421_985_545_5,
-                1.317_087_179_192_858_7,
-                0.891_167_948_026_476_6,
-                0.384_779_120_148_185_27,
-                0.155_955_856_853_634_36,
-                0.046_735_160_424_814_17,
-                0.014_576_803_693_379_53,
-                0.003_228_437_222_320_849,
-                0.000_914_252_863_201_387_9,
-                0.000_132_763_577_755_739_43,
-                0.000_048_420_884_483_658_204,
+            [
                 7.377_657_774_398_435e-7,
+                0.000_048_420_884_483_658_204,
+                0.000_132_763_577_755_739_43,
+                0.000_914_252_863_201_387_9,
+                0.003_228_437_222_320_849,
+                0.014_576_803_693_379_53,
+                0.046_735_160_424_814_17,
+                0.155_955_856_853_634_36,
+                0.384_779_120_148_185_27,
+                0.891_167_948_026_476_6,
+                1.317_087_179_192_858_7,
+                1.608_359_421_985_545_5,
             ],
         )
     } else {
diff --git a/src/utils.rs b/src/utils.rs
index a5cf61b..58c3583 100644
--- a/src/utils.rs
+++ b/src/utils.rs
@@ -106,12 +106,14 @@ pub fn sign(a: f64, b: f64) -> f64 {
 
 /// Evaluate a polynomial at x using [Horner's method](https://en.wikipedia.org/wiki/Horner%27s_method).
 ///
-/// The coefficients are assumed to be sorted in increasing order by the degree of the associated x-term.
+/// The iterator is assumed to give the coefficients in decreasing order by the degree of the associated x-term.
 ///
-/// For example, to evaluate 4x^2 + 3x + 2, use `coefficients = [2.0, 3.0, 4.0]`.
-pub fn polynomial<const N: usize>(x: f64, coefficients: &[f64; N]) -> f64 {
+/// For example, to evaluate 4x^2 + 3x + 2, the iterator should give 4.0 then 3.0, and finally 2.0.
+pub fn polynomial<I>(x: f64, coefficients: I) -> f64
+where
+    I: IntoIterator<Item = f64>,
+{
     coefficients
-        .iter()
-        .rev()
-        .fold(0.0, |acc, &coeff| acc * x + coeff)
+        .into_iter()
+        .fold(0.0, |acc, coeff| acc * x + coeff)
 }

From f293682919aba4898930b69693fd0723c1427529 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Fri, 31 Oct 2025 23:43:27 +0100
Subject: [PATCH 11/35] Rename input to `z`

---
 src/gamma.rs | 40 ++++++++++++++++++++--------------------
 1 file changed, 20 insertions(+), 20 deletions(-)

diff --git a/src/gamma.rs b/src/gamma.rs
index 8bcb445..8014b6a 100644
--- a/src/gamma.rs
+++ b/src/gamma.rs
@@ -55,16 +55,16 @@ pub fn ln_gamma(z: f64) -> f64 {
 ///
 /// The value of the gamma function at `z`
 // Based on the Fortran implementation by Toshio Fukushima: <https://www.researchgate.net/publication/336578125_xgamtxt_Fortran_90_test_program_package_of_qgam_dgam_and_sgam_fast_computation_of_Gamma_function_with_quadruple_double_and_single_precision_accuracy_respectively>
-pub fn gamma(mut x: f64) -> f64 {
+pub fn gamma(mut z: f64) -> f64 {
     /// This input gives an output of [`f64::MAX`].
     /// Computed with WolframAlpha: https://www.wolframalpha.com/input?i=solve+gamma%28x%29+%3D+%281+%E2%88%92+2%5E%28%E2%88%9253%29%29*2%5E1024+for+x+%3E+1.
     const MAX_INPUT: f64 = 171.624_376_956_302_7;
 
     // Special cases.
-    if x > MAX_INPUT {
+    if z > MAX_INPUT {
         // Output is too large to represent.
         return f64::INFINITY;
-    } else if x == 0.0 {
+    } else if z == 0.0 {
         // The gamma function diverges for an input of zero.
         // It diverges to positive or negative infinity depending on which direction zero is approached from.
         // According to entry F.9.5.4 of the ISO C 99 standard (https://www.open-std.org/jtc1/sc22/wg14/www/docs/n1256.pdf)
@@ -72,47 +72,47 @@ pub fn gamma(mut x: f64) -> f64 {
         // This is the standard followed by `scipy.special.gamma` as well as the `tgamma` function in `libc`.
         // Since the Gamma function in (nightly) Rust's standard library currently corresponds to `tgamma`
         // and thus follows the standard we do the same here.
-        return f64::INFINITY.copysign(x);
-    } else if x < 0.0 && x.fract() == 0.0 {
+        return f64::INFINITY.copysign(z);
+    } else if z < 0.0 && z.fract() == 0.0 {
         // The Gamma function diverges for non-positive integers.
         // There is however no clear way for the caller to signal whether the input approached the integer from above or below.
         // According to the same standard as above the gamma function should return NaN for these inputs.
         return f64::NAN;
     }
 
-    let f = if x > 3.5 {
+    let f = if z > 3.5 {
         let mut f = 1.0;
-        // Decrement x by 1 until less than 3.5.
+        // Decrement z by 1 until less than 3.5.
         // This will not be a massively long loop since we checked for too large input above.
-        while x >= 3.5 {
-            x -= 1.0;
-            f *= x;
+        while z >= 3.5 {
+            z -= 1.0;
+            f *= z;
         }
         f
-    } else if x < 2.5 {
+    } else if z < 2.5 {
         let mut f = 1.0;
-        // Increment x by 1 until larger than 2.5.
-        while x <= 2.5 {
-            f *= x;
+        // Increment z by 1 until larger than 2.5.
+        while z <= 2.5 {
+            f *= z;
             if f.is_infinite() {
                 // There is no point in continuing the calculation,
                 // as `f` will remain infinite when multiplied by any finite number.
-                // (except zero, but that can only happen if x is an integer, which we have handled above.)
+                // (except zero, but that can only happen if z is an integer, which we have handled above.)
                 // The final result will thus be `g` / inf,
                 // which underflows to 0 since `g` (defined below) is always finite.
                 // This early return stops the loop from running too long for very large negative inputs.
                 return 0.0;
             }
-            x += 1.0;
+            z += 1.0;
         }
         f.recip()
     } else {
         1.0
     };
 
-    let g = if x > 3.0 {
+    let g = if z > 3.0 {
         polynomial(
-            x - 3.25,
+            z - 3.25,
             [
                 0.000_016_738_398_919_923_317,
                 0.000_095_940_184_093_793_52,
@@ -128,9 +128,9 @@ pub fn gamma(mut x: f64) -> f64 {
                 2.549_256_966_718_529,
             ],
         )
-    } else if x < 3.0 {
+    } else if z < 3.0 {
         polynomial(
-            x - 2.75,
+            z - 2.75,
             [
                 7.377_657_774_398_435e-7,
                 0.000_048_420_884_483_658_204,

From 14e99bdfc0039afd3ac8d7c8fa18e18fd0d3c69f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Fri, 31 Oct 2025 23:46:08 +0100
Subject: [PATCH 12/35] Make the Gamma function give NaN for NaN inputs

---
 src/gamma.rs        | 2 ++
 tests/gamma_test.rs | 4 ++++
 2 files changed, 6 insertions(+)

diff --git a/src/gamma.rs b/src/gamma.rs
index 8014b6a..b024e49 100644
--- a/src/gamma.rs
+++ b/src/gamma.rs
@@ -78,6 +78,8 @@ pub fn gamma(mut z: f64) -> f64 {
         // There is however no clear way for the caller to signal whether the input approached the integer from above or below.
         // According to the same standard as above the gamma function should return NaN for these inputs.
         return f64::NAN;
+    } else if z.is_nan() {
+        return f64::NAN;
     }
 
     let f = if z > 3.5 {
diff --git a/tests/gamma_test.rs b/tests/gamma_test.rs
index 2fb6cd5..27e0c15 100644
--- a/tests/gamma_test.rs
+++ b/tests/gamma_test.rs
@@ -1,3 +1,5 @@
+use core::f64;
+
 use approx::assert_relative_eq;
 use proptest::prelude::*;
 use puruspe::{gamma, gammp, gammq, invgammp, ln_gamma};
@@ -67,6 +69,8 @@ fn test_invgammp() {
 fn test_gamma_edge_cases() {
     // Test gamma at special values
 
+    assert!(gamma(f64::NAN).is_nan());
+
     // Γ(±0) = ±∞
     assert_eq!(gamma(0.0), f64::INFINITY);
     assert_eq!(gamma(-0.0), -f64::INFINITY);

From 480ae02ebaa33bd7968c57fc8430b615d190dedf Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Fri, 31 Oct 2025 23:49:33 +0100
Subject: [PATCH 13/35] Also test an input of -inf

---
 src/gamma.rs        | 3 +++
 tests/gamma_test.rs | 2 ++
 2 files changed, 5 insertions(+)

diff --git a/src/gamma.rs b/src/gamma.rs
index b024e49..ad9a98e 100644
--- a/src/gamma.rs
+++ b/src/gamma.rs
@@ -78,6 +78,9 @@ pub fn gamma(mut z: f64) -> f64 {
         // There is however no clear way for the caller to signal whether the input approached the integer from above or below.
         // According to the same standard as above the gamma function should return NaN for these inputs.
         return f64::NAN;
+    } else if z == f64::NEG_INFINITY {
+        // There is no well defined limit for z --> -inf.
+        return f64::NAN;
     } else if z.is_nan() {
         return f64::NAN;
     }
diff --git a/tests/gamma_test.rs b/tests/gamma_test.rs
index 27e0c15..2f91349 100644
--- a/tests/gamma_test.rs
+++ b/tests/gamma_test.rs
@@ -70,6 +70,8 @@ fn test_gamma_edge_cases() {
     // Test gamma at special values
 
     assert!(gamma(f64::NAN).is_nan());
+    assert!(gamma(f64::NEG_INFINITY).is_nan());
+    assert_eq!(gamma(f64::INFINITY), f64::INFINITY);
 
     // Γ(±0) = ±∞
     assert_eq!(gamma(0.0), f64::INFINITY);

From 409cd8e0b44ef19ada72485e1edbca66a010d7b5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Fri, 31 Oct 2025 23:51:28 +0100
Subject: [PATCH 14/35] Comment the special cases also in the tests

---
 tests/gamma_test.rs | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/tests/gamma_test.rs b/tests/gamma_test.rs
index 2f91349..1c099d1 100644
--- a/tests/gamma_test.rs
+++ b/tests/gamma_test.rs
@@ -69,8 +69,12 @@ fn test_invgammp() {
 fn test_gamma_edge_cases() {
     // Test gamma at special values
 
+    // Can't do anything with NaNs. Just propagate it.
     assert!(gamma(f64::NAN).is_nan());
+
+    // No limit for Γ(x) as x → -∞
     assert!(gamma(f64::NEG_INFINITY).is_nan());
+    // Γ(+∞) = +∞
     assert_eq!(gamma(f64::INFINITY), f64::INFINITY);
 
     // Γ(±0) = ±∞

From b5bbe01678814037e77b10308b43dd30055fa0ab Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Fri, 31 Oct 2025 23:54:18 +0100
Subject: [PATCH 15/35] Better grammar in comment

---
 src/gamma.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/gamma.rs b/src/gamma.rs
index ad9a98e..4a3ad71 100644
--- a/src/gamma.rs
+++ b/src/gamma.rs
@@ -74,7 +74,7 @@ pub fn gamma(mut z: f64) -> f64 {
         // and thus follows the standard we do the same here.
         return f64::INFINITY.copysign(z);
     } else if z < 0.0 && z.fract() == 0.0 {
-        // The Gamma function diverges for non-positive integers.
+        // The Gamma function also diverges for negative integers.
         // There is however no clear way for the caller to signal whether the input approached the integer from above or below.
         // According to the same standard as above the gamma function should return NaN for these inputs.
         return f64::NAN;

From 8af2e3751ce0ad688655fe07d8a73ba979086f4b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Fri, 31 Oct 2025 23:56:51 +0100
Subject: [PATCH 16/35] Make link in comment easier to just copy paste

---
 src/gamma.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/gamma.rs b/src/gamma.rs
index 4a3ad71..871d519 100644
--- a/src/gamma.rs
+++ b/src/gamma.rs
@@ -54,7 +54,7 @@ pub fn ln_gamma(z: f64) -> f64 {
 /// # Returns
 ///
 /// The value of the gamma function at `z`
-// Based on the Fortran implementation by Toshio Fukushima: <https://www.researchgate.net/publication/336578125_xgamtxt_Fortran_90_test_program_package_of_qgam_dgam_and_sgam_fast_computation_of_Gamma_function_with_quadruple_double_and_single_precision_accuracy_respectively>
+// Based on the Fortran implementation by Toshio Fukushima: https://www.researchgate.net/publication/336578125_xgamtxt_Fortran_90_test_program_package_of_qgam_dgam_and_sgam_fast_computation_of_Gamma_function_with_quadruple_double_and_single_precision_accuracy_respectively
 pub fn gamma(mut z: f64) -> f64 {
     /// This input gives an output of [`f64::MAX`].
     /// Computed with WolframAlpha: https://www.wolframalpha.com/input?i=solve+gamma%28x%29+%3D+%281+%E2%88%92+2%5E%28%E2%88%9253%29%29*2%5E1024+for+x+%3E+1.

From 8b56f640349c0c845cfcd7e36a7c457923f844f1 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Sat, 1 Nov 2025 00:00:45 +0100
Subject: [PATCH 17/35] Use correct name of input in comment

---
 tests/gamma_test.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/gamma_test.rs b/tests/gamma_test.rs
index 1c099d1..925d72b 100644
--- a/tests/gamma_test.rs
+++ b/tests/gamma_test.rs
@@ -72,7 +72,7 @@ fn test_gamma_edge_cases() {
     // Can't do anything with NaNs. Just propagate it.
     assert!(gamma(f64::NAN).is_nan());
 
-    // No limit for Γ(x) as x → -∞
+    // No limit for Γ(z) as z → -∞
     assert!(gamma(f64::NEG_INFINITY).is_nan());
     // Γ(+∞) = +∞
     assert_eq!(gamma(f64::INFINITY), f64::INFINITY);

From 7b823ce8b3747f7e30db20d4d5935ce1cf68eabb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Mon, 3 Nov 2025 19:45:22 +0100
Subject: [PATCH 18/35] Add a test at `f64::MIN_POSITIVE`

---
 tests/gamma_test.rs | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/tests/gamma_test.rs b/tests/gamma_test.rs
index 925d72b..f58b52e 100644
--- a/tests/gamma_test.rs
+++ b/tests/gamma_test.rs
@@ -104,6 +104,10 @@ fn test_gamma_edge_cases() {
         99999999.42278434,
         max_relative = 1.5 * f64::EPSILON
     );
+    assert_relative_eq!(
+        gamma(f64::MIN_POSITIVE),
+        4.4942328371557897e307,
+    );
 
     // Test negative values
     // Γ(-0.5) = -2√π

From 9cbc0ac6f8f65bf0e1e641d83f0088b8b85b157b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Mon, 3 Nov 2025 19:45:40 +0100
Subject: [PATCH 19/35] `cargo fmt`

---
 tests/gamma_test.rs | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/tests/gamma_test.rs b/tests/gamma_test.rs
index f58b52e..482bc3f 100644
--- a/tests/gamma_test.rs
+++ b/tests/gamma_test.rs
@@ -104,10 +104,7 @@ fn test_gamma_edge_cases() {
         99999999.42278434,
         max_relative = 1.5 * f64::EPSILON
     );
-    assert_relative_eq!(
-        gamma(f64::MIN_POSITIVE),
-        4.4942328371557897e307,
-    );
+    assert_relative_eq!(gamma(f64::MIN_POSITIVE), 4.4942328371557897e307,);
 
     // Test negative values
     // Γ(-0.5) = -2√π

From 33ff53c0725667d04db0df3078647a14e8cc580a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Mon, 3 Nov 2025 19:45:52 +0100
Subject: [PATCH 20/35] Remove a comma

---
 tests/gamma_test.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/gamma_test.rs b/tests/gamma_test.rs
index 482bc3f..3f3d6b5 100644
--- a/tests/gamma_test.rs
+++ b/tests/gamma_test.rs
@@ -104,7 +104,7 @@ fn test_gamma_edge_cases() {
         99999999.42278434,
         max_relative = 1.5 * f64::EPSILON
     );
-    assert_relative_eq!(gamma(f64::MIN_POSITIVE), 4.4942328371557897e307,);
+    assert_relative_eq!(gamma(f64::MIN_POSITIVE), 4.4942328371557897e307);
 
     // Test negative values
     // Γ(-0.5) = -2√π

From 6302264c5d1ddcded360459b5873011d6e2aa795 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Mon, 3 Nov 2025 19:46:39 +0100
Subject: [PATCH 21/35] It even gets it exactly right!

---
 tests/gamma_test.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/gamma_test.rs b/tests/gamma_test.rs
index 3f3d6b5..74af9c8 100644
--- a/tests/gamma_test.rs
+++ b/tests/gamma_test.rs
@@ -104,7 +104,7 @@ fn test_gamma_edge_cases() {
         99999999.42278434,
         max_relative = 1.5 * f64::EPSILON
     );
-    assert_relative_eq!(gamma(f64::MIN_POSITIVE), 4.4942328371557897e307);
+    assert_eq!(gamma(f64::MIN_POSITIVE), 4.4942328371557897e307);
 
     // Test negative values
     // Γ(-0.5) = -2√π

From 7feae764863e33d377287db032f3603b0949b2f5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Mon, 3 Nov 2025 19:47:17 +0100
Subject: [PATCH 22/35] Add comment with source

---
 tests/gamma_test.rs | 1 +
 1 file changed, 1 insertion(+)

diff --git a/tests/gamma_test.rs b/tests/gamma_test.rs
index 74af9c8..8f79279 100644
--- a/tests/gamma_test.rs
+++ b/tests/gamma_test.rs
@@ -104,6 +104,7 @@ fn test_gamma_edge_cases() {
         99999999.42278434,
         max_relative = 1.5 * f64::EPSILON
     );
+    // Computed with WolframAlpha: https://www.wolframalpha.com/input?i=Gamma%5B2.2250738585072014E-308%5D.
     assert_eq!(gamma(f64::MIN_POSITIVE), 4.4942328371557897e307);
 
     // Test negative values

From 185f2b2e702c8fb44c4e1649ae6c4057db3859a1 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Mon, 3 Nov 2025 19:48:19 +0100
Subject: [PATCH 23/35] Add back the relative equality

I do not have a guarantee that wolfram alpha has it exactly right to the exact floating point number
---
 tests/gamma_test.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/gamma_test.rs b/tests/gamma_test.rs
index 8f79279..a947115 100644
--- a/tests/gamma_test.rs
+++ b/tests/gamma_test.rs
@@ -105,7 +105,7 @@ fn test_gamma_edge_cases() {
         max_relative = 1.5 * f64::EPSILON
     );
     // Computed with WolframAlpha: https://www.wolframalpha.com/input?i=Gamma%5B2.2250738585072014E-308%5D.
-    assert_eq!(gamma(f64::MIN_POSITIVE), 4.4942328371557897e307);
+    assert_relative_eq!(gamma(f64::MIN_POSITIVE), 4.4942328371557897e307);
 
     // Test negative values
     // Γ(-0.5) = -2√π

From 0f7c9b26fbe7a5dc2ee77eb588aa9ce8c1f9d797 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Mon, 3 Nov 2025 19:49:23 +0100
Subject: [PATCH 24/35] Move test to be in a groups with all other similar
 tests

---
 tests/gamma_test.rs | 14 ++++++--------
 1 file changed, 6 insertions(+), 8 deletions(-)

diff --git a/tests/gamma_test.rs b/tests/gamma_test.rs
index a947115..ddab998 100644
--- a/tests/gamma_test.rs
+++ b/tests/gamma_test.rs
@@ -104,6 +104,12 @@ fn test_gamma_edge_cases() {
         99999999.42278434,
         max_relative = 1.5 * f64::EPSILON
     );
+    // Computed with WolframAlpha: https://www.wolframalpha.com/input?i=Gamma%281e-15%29
+    assert_relative_eq!(
+        gamma(1e-15),
+        999999999999999.4,
+        max_relative = 2.0 * f64::EPSILON
+    );
     // Computed with WolframAlpha: https://www.wolframalpha.com/input?i=Gamma%5B2.2250738585072014E-308%5D.
     assert_relative_eq!(gamma(f64::MIN_POSITIVE), 4.4942328371557897e307);
 
@@ -114,14 +120,6 @@ fn test_gamma_edge_cases() {
     // Test near-integer negative values
     assert_relative_eq!(gamma(-1.5), 4.0 * sqrt_pi / 3.0);
 
-    // Test value very close to 0 but positive
-    // Computed with WolframAlpha: https://www.wolframalpha.com/input?i=Gamma%281e-15%29
-    assert_relative_eq!(
-        gamma(1e-15),
-        999999999999999.4,
-        max_relative = 2.0 * f64::EPSILON
-    );
-
     // Test integer factorials
     // Γ(6) = 5! = 120
     assert_eq!(gamma(6.0), 120.0);

From 372961bc87dafba88bcff3dbf259a1372083235f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Mon, 3 Nov 2025 19:52:43 +0100
Subject: [PATCH 25/35] Test some more integer factorials

---
 tests/gamma_test.rs | 12 ++++--------
 1 file changed, 4 insertions(+), 8 deletions(-)

diff --git a/tests/gamma_test.rs b/tests/gamma_test.rs
index ddab998..1b178db 100644
--- a/tests/gamma_test.rs
+++ b/tests/gamma_test.rs
@@ -1,8 +1,6 @@
-use core::f64;
-
 use approx::assert_relative_eq;
 use proptest::prelude::*;
-use puruspe::{gamma, gammp, gammq, invgammp, ln_gamma};
+use puruspe::{gamma, gammp, gammq, invgammp, ln_gamma, utils::factorial};
 
 unsafe extern "C" {
     fn tgamma(x: f64) -> f64;
@@ -121,11 +119,9 @@ fn test_gamma_edge_cases() {
     assert_relative_eq!(gamma(-1.5), 4.0 * sqrt_pi / 3.0);
 
     // Test integer factorials
-    // Γ(6) = 5! = 120
-    assert_eq!(gamma(6.0), 120.0);
-
-    // Γ(11) = 10! = 3628800
-    assert_eq!(gamma(11.0), 3628800.0);
+    for i in 1..=25 {
+        assert_eq!(gamma(i as f64), factorial(i - 1));
+    }
 }
 
 #[test]

From 2439859971ea92c1948f43f29ea4803da35109ab Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Mon, 3 Nov 2025 21:42:28 +0100
Subject: [PATCH 26/35] Specialize on integers and calculate Gamma(integer)
 with factorial

---
 src/gamma.rs        | 19 +++++++++++++------
 tests/gamma_test.rs |  2 +-
 2 files changed, 14 insertions(+), 7 deletions(-)

diff --git a/src/gamma.rs b/src/gamma.rs
index 871d519..4f65982 100644
--- a/src/gamma.rs
+++ b/src/gamma.rs
@@ -7,7 +7,7 @@
 //! - `gammq`: Calculates the regularized upper incomplete gamma function.
 //! - `invgammp`: Calculates the inverse of the regularized lower incomplete gamma function.
 
-use crate::utils::polynomial;
+use crate::utils::{factorial, polynomial};
 use crate::{EPS, FPMIN, W, Y};
 use core::f64::consts::PI;
 const ASWITCH: usize = 100;
@@ -73,11 +73,18 @@ pub fn gamma(mut z: f64) -> f64 {
         // Since the Gamma function in (nightly) Rust's standard library currently corresponds to `tgamma`
         // and thus follows the standard we do the same here.
         return f64::INFINITY.copysign(z);
-    } else if z < 0.0 && z.fract() == 0.0 {
-        // The Gamma function also diverges for negative integers.
-        // There is however no clear way for the caller to signal whether the input approached the integer from above or below.
-        // According to the same standard as above the gamma function should return NaN for these inputs.
-        return f64::NAN;
+    } else if z.fract() == 0.0 {
+        // z is a non-zero integer less than or equal to 171.
+
+        if z < 0.0 {
+            // The Gamma function also diverges for negative integers.
+            // There is however no clear way for the caller to signal whether the input approached the integer from above or below.
+            // According to the same standard as above the gamma function should return NaN for these inputs.
+            return f64::NAN;
+        }
+
+        // 0 < z <= 171, we can just use the factorial.
+        return factorial(z as usize - 1);
     } else if z == f64::NEG_INFINITY {
         // There is no well defined limit for z --> -inf.
         return f64::NAN;
diff --git a/tests/gamma_test.rs b/tests/gamma_test.rs
index 1b178db..d0aefef 100644
--- a/tests/gamma_test.rs
+++ b/tests/gamma_test.rs
@@ -119,7 +119,7 @@ fn test_gamma_edge_cases() {
     assert_relative_eq!(gamma(-1.5), 4.0 * sqrt_pi / 3.0);
 
     // Test integer factorials
-    for i in 1..=25 {
+    for i in 1..=172 {
         assert_eq!(gamma(i as f64), factorial(i - 1));
     }
 }

From 70e89b9cc342d036a9a22782ef94df0437fa3128 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?=
 <44257381+JSorngard@users.noreply.github.com>
Date: Fri, 7 Nov 2025 13:17:23 +0100
Subject: [PATCH 27/35] Use fused multiply-add instructions in `polynomial` if
 they are available

---
 src/utils.rs | 17 ++++++++++++++++-
 1 file changed, 16 insertions(+), 1 deletion(-)

diff --git a/src/utils.rs b/src/utils.rs
index 58c3583..765313a 100644
--- a/src/utils.rs
+++ b/src/utils.rs
@@ -115,5 +115,20 @@ where
 {
     coefficients
         .into_iter()
-        .fold(0.0, |acc, coeff| acc * x + coeff)
+        .fold(0.0, |acc, coeff| mul_add(acc, x, coeff)
 }
+
+
+/// Multiplies `x` by `mul` and adds `add`.
+/// If the target CPU supports fused multiply-add instructions this function will use those.
+fn mul_add(x: f64, mul: f64, add: f64) -> f64 {
+    #[cfg(not(target_feature = "fma"))] 
+    {
+        x * mul + add
+    }
+
+    #[cfg(target_feature = "fma")]
+    {
+        x.mul_add(mul, add)
+    }
+}
\ No newline at end of file

From a4439138ff978ed4016695efb28cfab680c06969 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?=
 <44257381+JSorngard@users.noreply.github.com>
Date: Fri, 7 Nov 2025 13:19:35 +0100
Subject: [PATCH 28/35] Correct unclosed paren

---
 src/utils.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/utils.rs b/src/utils.rs
index 765313a..d7688e2 100644
--- a/src/utils.rs
+++ b/src/utils.rs
@@ -115,7 +115,7 @@ where
 {
     coefficients
         .into_iter()
-        .fold(0.0, |acc, coeff| mul_add(acc, x, coeff)
+        .fold(0.0, |acc, coeff| mul_add(acc, x, coeff))
 }
 
 

From e34459df29ded1295c610247d7214732bf34358a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?=
 <44257381+JSorngard@users.noreply.github.com>
Date: Fri, 7 Nov 2025 13:21:04 +0100
Subject: [PATCH 29/35] Fix some whitespace

---
 src/utils.rs | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/utils.rs b/src/utils.rs
index d7688e2..e1d6390 100644
--- a/src/utils.rs
+++ b/src/utils.rs
@@ -122,7 +122,7 @@ where
 /// Multiplies `x` by `mul` and adds `add`.
 /// If the target CPU supports fused multiply-add instructions this function will use those.
 fn mul_add(x: f64, mul: f64, add: f64) -> f64 {
-    #[cfg(not(target_feature = "fma"))] 
+    #[cfg(not(target_feature = "fma"))]
     {
         x * mul + add
     }
@@ -131,4 +131,4 @@ fn mul_add(x: f64, mul: f64, add: f64) -> f64 {
     {
         x.mul_add(mul, add)
     }
-}
\ No newline at end of file
+}

From e4e5304fb2253d30d2db8efb06fb3385ce2f553c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?=
 <44257381+JSorngard@users.noreply.github.com>
Date: Fri, 7 Nov 2025 13:22:05 +0100
Subject: [PATCH 30/35] Remove a newline

---
 src/utils.rs | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/utils.rs b/src/utils.rs
index e1d6390..cb7e19e 100644
--- a/src/utils.rs
+++ b/src/utils.rs
@@ -118,7 +118,6 @@ where
         .fold(0.0, |acc, coeff| mul_add(acc, x, coeff))
 }
 
-
 /// Multiplies `x` by `mul` and adds `add`.
 /// If the target CPU supports fused multiply-add instructions this function will use those.
 fn mul_add(x: f64, mul: f64, add: f64) -> f64 {

From 49bd5308ea38793d92bae589f2345d8410cf26c0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Wed, 12 Nov 2025 19:49:18 +0100
Subject: [PATCH 31/35] Add a doc-example to `polynomial`

---
 src/utils.rs | 17 +++++++++++++++--
 1 file changed, 15 insertions(+), 2 deletions(-)

diff --git a/src/utils.rs b/src/utils.rs
index cb7e19e..15ba71e 100644
--- a/src/utils.rs
+++ b/src/utils.rs
@@ -108,8 +108,21 @@ pub fn sign(a: f64, b: f64) -> f64 {
 ///
 /// The iterator is assumed to give the coefficients in decreasing order by the degree of the associated x-term.
 ///
-/// For example, to evaluate 4x^2 + 3x + 2, the iterator should give 4.0 then 3.0, and finally 2.0.
-pub fn polynomial<I>(x: f64, coefficients: I) -> f64
+/// # Example
+/// 
+/// Evaluate 4x^2 + 3x + 2 at x = 2:
+/// 
+/// ```
+/// use puruspe::utils::polynomial;
+/// 
+/// let x = 2.0;
+/// assert_eq!(
+///     polynomial(x, [4.0, 3.0, 2.0]),
+///     4.0*x*x + 3.0*x + 2.0,
+/// );
+/// ```
+// Should we change this `pub(crate)` to a `pub` and export it to the user?
+pub(crate) fn polynomial<I>(x: f64, coefficients: I) -> f64
 where
     I: IntoIterator<Item = f64>,
 {

From a6a9aeade08a75e2c4ccd57bf003ae689e3d9ebb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Wed, 12 Nov 2025 19:49:27 +0100
Subject: [PATCH 32/35] `cargo fmt`

---
 src/utils.rs | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/utils.rs b/src/utils.rs
index 15ba71e..e4392c4 100644
--- a/src/utils.rs
+++ b/src/utils.rs
@@ -109,12 +109,12 @@ pub fn sign(a: f64, b: f64) -> f64 {
 /// The iterator is assumed to give the coefficients in decreasing order by the degree of the associated x-term.
 ///
 /// # Example
-/// 
+///
 /// Evaluate 4x^2 + 3x + 2 at x = 2:
-/// 
+///
 /// ```
 /// use puruspe::utils::polynomial;
-/// 
+///
 /// let x = 2.0;
 /// assert_eq!(
 ///     polynomial(x, [4.0, 3.0, 2.0]),

From 6063b5bff59af9d6eb7ea59cd354435087060385 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Wed, 12 Nov 2025 19:50:23 +0100
Subject: [PATCH 33/35] Remove the doc-example since the function is not public

---
 src/utils.rs | 14 --------------
 1 file changed, 14 deletions(-)

diff --git a/src/utils.rs b/src/utils.rs
index e4392c4..5a43476 100644
--- a/src/utils.rs
+++ b/src/utils.rs
@@ -107,20 +107,6 @@ pub fn sign(a: f64, b: f64) -> f64 {
 /// Evaluate a polynomial at x using [Horner's method](https://en.wikipedia.org/wiki/Horner%27s_method).
 ///
 /// The iterator is assumed to give the coefficients in decreasing order by the degree of the associated x-term.
-///
-/// # Example
-///
-/// Evaluate 4x^2 + 3x + 2 at x = 2:
-///
-/// ```
-/// use puruspe::utils::polynomial;
-///
-/// let x = 2.0;
-/// assert_eq!(
-///     polynomial(x, [4.0, 3.0, 2.0]),
-///     4.0*x*x + 3.0*x + 2.0,
-/// );
-/// ```
 // Should we change this `pub(crate)` to a `pub` and export it to the user?
 pub(crate) fn polynomial<I>(x: f64, coefficients: I) -> f64
 where

From 7044241bf22c48181991c53765b59c22c104b87d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Wed, 12 Nov 2025 19:51:06 +0100
Subject: [PATCH 34/35] Re-add the explanatory comment

---
 src/utils.rs | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/utils.rs b/src/utils.rs
index 5a43476..27e9176 100644
--- a/src/utils.rs
+++ b/src/utils.rs
@@ -107,6 +107,8 @@ pub fn sign(a: f64, b: f64) -> f64 {
 /// Evaluate a polynomial at x using [Horner's method](https://en.wikipedia.org/wiki/Horner%27s_method).
 ///
 /// The iterator is assumed to give the coefficients in decreasing order by the degree of the associated x-term.
+///
+/// For example, to evaluate 4x^2 + 3x + 2, the iterator should give 4.0 then 3.0, and finally 2.0.
 // Should we change this `pub(crate)` to a `pub` and export it to the user?
 pub(crate) fn polynomial<I>(x: f64, coefficients: I) -> f64
 where

From deff23519caabc0c270db70a3de9d9471c806628 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Johanna=20S=C3=B6rng=C3=A5rd?= <jsorngard@gmail.com>
Date: Wed, 12 Nov 2025 19:52:12 +0100
Subject: [PATCH 35/35] Remove comment that was moved to github PR

---
 src/utils.rs | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/utils.rs b/src/utils.rs
index 27e9176..b8b90c7 100644
--- a/src/utils.rs
+++ b/src/utils.rs
@@ -109,7 +109,6 @@ pub fn sign(a: f64, b: f64) -> f64 {
 /// The iterator is assumed to give the coefficients in decreasing order by the degree of the associated x-term.
 ///
 /// For example, to evaluate 4x^2 + 3x + 2, the iterator should give 4.0 then 3.0, and finally 2.0.
-// Should we change this `pub(crate)` to a `pub` and export it to the user?
 pub(crate) fn polynomial<I>(x: f64, coefficients: I) -> f64
 where
     I: IntoIterator<Item = f64>,