add incomplete gamma function

CHANGELOG.md

@@ -7,4 +7,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+### Added
+- Added implementations for incomplete gamma functions, converted from Math.Net
+
 [Unreleased]: https://git.chir.rs/ProcyOS/rust-template

CONTRIBUTING.md

@@ -1,4 +1,4 @@
-# Contributing to [Project]
+# Contributing to `extra-math`
 
 Thank you for considering to contribute to our project!
 

Cargo.lock

@@ -3,5 +3,31 @@
 version = 3
 
 [[package]]
-name = "rust-template"
+name = "autocfg"
+version = "1.3.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "0c4b4d0bd25bd0b74681c0ad21497610ce1b7c91b1022cd21c80c6fbdd9476b0"
+
+[[package]]
+name = "extra-math"
 version = "0.1.0"
+dependencies = [
+ "libm",
+ "num-traits",
+]
+
+[[package]]
+name = "libm"
+version = "0.2.8"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "4ec2a862134d2a7d32d7983ddcdd1c4923530833c9f2ea1a44fc5fa473989058"
+
+[[package]]
+name = "num-traits"
+version = "0.2.19"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841"
+dependencies = [
+ "autocfg",
+ "libm",
+]

Cargo.nix

@@ -3,7 +3,7 @@
 args @ {
   release ? true,
   rootFeatures ? [
-    "rust-template/default"
+    "extra-math/default"
   ],
   rustPackages,
   buildRustPackages,
@@ -22,7 +22,7 @@ args @ {
   workspaceSrc,
   ignoreLockHash,
 }: let
-  nixifiedLockHash = "4598ab22b158163d169b6aea289131643ca58806898bb56966fba495f4660c02";
+  nixifiedLockHash = "9a1ee9305c225b74d8dc742d0cdc3b6fa8f016ddd5427a73b6ec00ee90d7f0f7";
   workspaceSrc =
     if args.workspaceSrc == null
     then ./.
@@ -59,12 +59,63 @@ in
   in {
     cargo2nixVersion = "0.11.0";
     workspace = {
-      rust-template = rustPackages.unknown.rust-template."0.1.0";
+      extra-math = rustPackages.unknown.extra-math."0.1.0";
     };
-    "unknown".rust-template."0.1.0" = overridableMkRustCrate (profileName: rec {
-      name = "rust-template";
+    "registry+https://github.com/rust-lang/crates.io-index".autocfg."1.3.0" = overridableMkRustCrate (profileName: rec {
+      name = "autocfg";
+      version = "1.3.0";
+      registry = "registry+https://github.com/rust-lang/crates.io-index";
+      src = fetchCratesIo {
+        inherit name version;
+        sha256 = "0c4b4d0bd25bd0b74681c0ad21497610ce1b7c91b1022cd21c80c6fbdd9476b0";
+      };
+    });
+
+    "unknown".extra-math."0.1.0" = overridableMkRustCrate (profileName: rec {
+      name = "extra-math";
       version = "0.1.0";
       registry = "unknown";
       src = fetchCrateLocal workspaceSrc;
+      features = builtins.concatLists [
+        (lib.optional (rootFeatures' ? "extra-math/default") "default")
+        (lib.optional (rootFeatures' ? "extra-math/std") "std")
+      ];
+      dependencies = {
+        libm = (rustPackages."registry+https://github.com/rust-lang/crates.io-index".libm."0.2.8" {inherit profileName;}).out;
+        num_traits = (rustPackages."registry+https://github.com/rust-lang/crates.io-index".num-traits."0.2.19" {inherit profileName;}).out;
+      };
+    });
+
+    "registry+https://github.com/rust-lang/crates.io-index".libm."0.2.8" = overridableMkRustCrate (profileName: rec {
+      name = "libm";
+      version = "0.2.8";
+      registry = "registry+https://github.com/rust-lang/crates.io-index";
+      src = fetchCratesIo {
+        inherit name version;
+        sha256 = "4ec2a862134d2a7d32d7983ddcdd1c4923530833c9f2ea1a44fc5fa473989058";
+      };
+      features = builtins.concatLists [
+        ["default"]
+      ];
+    });
+
+    "registry+https://github.com/rust-lang/crates.io-index".num-traits."0.2.19" = overridableMkRustCrate (profileName: rec {
+      name = "num-traits";
+      version = "0.2.19";
+      registry = "registry+https://github.com/rust-lang/crates.io-index";
+      src = fetchCratesIo {
+        inherit name version;
+        sha256 = "071dfc062690e90b734c0b2273ce72ad0ffa95f0c74596bc250dcfd960262841";
+      };
+      features = builtins.concatLists [
+        ["libm"]
+        (lib.optional (rootFeatures' ? "extra-math/std") "std")
+      ];
+      dependencies = {
+        libm = (rustPackages."registry+https://github.com/rust-lang/crates.io-index".libm."0.2.8" {inherit profileName;}).out;
+      };
+      buildDependencies = {
+        autocfg = (buildRustPackages."registry+https://github.com/rust-lang/crates.io-index".autocfg."1.3.0" {profileName = "__noProfile";}).out;
+      };
     });
   }

Cargo.toml

@@ -1,15 +1,23 @@
 [package]
-name = "rust-template"
+name = "extra-math"
 version = "0.1.0"
 edition = "2021"
 authors = ["Charlotte 🦝 Delenk <lotte@chir.rs>"]
-description = "Template for rust libraries"
+description = "Additional math functions not found in the standard library or libm"
 license = "MIT OR Apache-2.0"
 repository = "https://git.chir.rs/ProcyOS/rust-template"
-keywords = ["template"]
-categories = ["template"]
+keywords = ["mathematics"]
+categories = ["mathematics", "no-std", "no-std::no-alloc"]
 
 [dependencies]
+libm = "0.2.8"
+num-traits = { version = "0.2.19", default-features = false, features = [
+    "libm",
+] }
+
+[features]
+default = []
+std = ["num-traits/std"]
 
 [lints.rust]
 deprecated-safe = "forbid"
@@ -33,8 +41,7 @@ trivial-numeric-casts = "warn"
 unit-bindings = "deny"
 unnameable-types = "warn"
 unreachable-pub = "warn"
-unsafe-code = "deny"
-unsafe-op-in-unsafe-fn = "warn"
+unsafe-code = "forbid"
 unused-crate-dependencies = "warn"
 unused-extern-crates = "warn"
 unused-import-braces = "warn"
@@ -62,5 +69,4 @@ rc-buffer = "warn"
 rc-mutex = "deny"
 std-instead-of-alloc = "warn"
 std-instead-of-core = "warn"
-undocumented-unsafe-blocks = "deny"
 unwrap-used = "forbid"

README.md

@@ -1,64 +1,21 @@
-# Project Title
+# extra_math
 
-One Paragraph of the project description
+Mathematical functions that are used in procyos’s code that aren’t found in either `libm` or the standard library.
 
-Initially appeared on
-[gist](https://gist.github.com/PurpleBooth/109311bb0361f32d87a2). But the page cannot open anymore so that is why I have moved it here.
-
-## Getting Started
-
-These instructions will give you a copy of the project up and running on
-your local machine for development and testing purposes. See deployment
-for notes on deploying the project on a live system.
-
-### Prerequisites
-
-Requirements for the software and other tools to build, test and push 
-- [Example 1](https://www.example.com)
-- [Example 2](https://www.example.com)
-
-### Installing
-
-A step by step series of examples that tell you how to get a development
-environment running
-
-Say what the step will be
-
-    Give the example
-
-And repeat
-
-    until finished
-
-End with an example of getting some data out of the system or using it
-for a little demo
+It is not intended to be accurate, just useful in a no-std environment.
 
 ## Running the tests
 
-Explain how to run the automated tests for this system
-
-### Sample Tests
-
-Explain what these tests test and why
-
-    Give an example
-
-### Style test
-
-Checks if the best practices and the right coding style has been used.
-
-    Give an example
-
-## Deployment
-
-Add additional notes to deploy this on a live system
+```sh
+cargo test
+nix flake check
+```
 
 ## Built With
 
   - [Contributor Covenant](https://www.contributor-covenant.org/) - Used
     for the Code of Conduct
-  - [Creative Commons](https://creativecommons.org/) - Used to choose
-    the license
+  - [Math.net](https://github.com/mathnet/mathnet-numerics/blob/master/src/Numerics/SpecialFunctions/Gamma.cs) for some of the code
 
 ## Contributing
 
@@ -90,9 +47,3 @@ Licensed under either of
    [LICENSE-MIT](LICENSE-MIT)
 
 at your option.
-
-## Acknowledgments
-
-  - Hat tip to anyone whose code is used
-  - Inspiration
-  - etc

flake.nix

@@ -1,5 +1,5 @@
 {
-  description = "rust-template";
+  description = "extra-math";
 
   inputs = {
     nixpkgs.url = github:NixOS/nixpkgs;

src/gamma.rs

@@ -0,0 +1,67 @@
+//! Gamma function approximation
+//!
+//! This module calculates various gamma functions
+
+use num_traits::Float;
+
+mod f64;
+
+/// Trait implementing the gamma functions
+pub trait Gamma: Float {
+    /// Calculates the logarithm of the Gamma function.
+    fn ln_gamma(self) -> Self {
+        self.gamma().ln()
+    }
+    /// Calculates the Gamma function.
+    fn gamma(self) -> Self;
+    /// Returns the upper incomplete regularized gamma function Q(a, x)
+    fn upper_gamma_regularized(self, x: Self) -> Self {
+        Self::one() - self.lower_gamma_regularized(x)
+    }
+    /// Returns the upper incomplete gamma function Q(a, x)
+    fn upper_gamma_incomplete(self, x: Self) -> Self {
+        self.upper_gamma_regularized(x) * self.gamma()
+    }
+    /// Returns the lower incomplete regularized gamma function P(a, x)
+    fn lower_gamma_regularized(self, x: Self) -> Self;
+    /// Returns the lower incomplete gamma function P(a, x)
+    fn lower_gamma_incomplete(self, x: Self) -> Self {
+        self.lower_gamma_regularized(x) * self.gamma()
+    }
+}
+
+impl Gamma for f32 {
+    fn gamma(self) -> Self {
+        Gamma::gamma(self as f64) as f32
+    }
+
+    fn lower_gamma_regularized(self, x: Self) -> Self {
+        Gamma::lower_gamma_regularized(self as f64, x as f64) as f32
+    }
+
+    fn ln_gamma(self) -> Self {
+        Gamma::ln_gamma(self as f64) as f32
+    }
+
+    fn upper_gamma_regularized(self, x: Self) -> Self {
+        Gamma::upper_gamma_regularized(self as f64, x as f64) as f32
+    }
+
+    fn upper_gamma_incomplete(self, x: Self) -> Self {
+        Gamma::upper_gamma_incomplete(self as f64, x as f64) as f32
+    }
+
+    fn lower_gamma_incomplete(self, x: Self) -> Self {
+        Gamma::lower_gamma_incomplete(self as f64, x as f64) as f32
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use crate::gamma::Gamma;
+
+    #[test]
+    fn upper_gamma_accuracy() {
+        assert!((Gamma::upper_gamma_regularized(1.5, 0.5) - 0.801252f64).abs() < 0.001);
+    }
+}

src/gamma/f64.rs

@@ -0,0 +1,239 @@
+//! 64 bit float implementation for gamma
+//!
+//! This implementation is based on the Math.NET library.
+
+// Copyright (c) 2009-2010 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+
+use core::f64::NAN;
+#[allow(unused_imports)]
+use num_traits::Float;
+
+use super::Gamma;
+
+impl Gamma for f64 {
+    fn ln_gamma(self) -> Self {
+        libm::lgamma(self)
+    }
+    fn gamma(self) -> Self {
+        libm::tgamma(self)
+    }
+    fn lower_gamma_regularized(self, x: Self) -> Self {
+        if self.is_nan() {
+            return NAN;
+        }
+        if x.is_nan() {
+            return NAN;
+        }
+        const EPSILON: f64 = 0.000000000000001;
+        const BIG: f64 = 4503599627370496.0;
+        const BIG_INV: f64 = 2.22044604925031308085e-16;
+        if self < 0.0 {
+            return NAN;
+        }
+        if x < 0.0 {
+            return NAN;
+        }
+        if self < 0.000000000000001 {
+            return 1.0;
+        }
+        if x < 0.000000000000001 {
+            return 0.0;
+        }
+        let ax = self * x.ln() - x - Gamma::ln_gamma(self);
+        if ax < -709.78271289338399 {
+            return if self < x { 1.0 } else { 0.0 };
+        }
+        if x <= 1.0 || x <= self {
+            let mut r2 = self;
+            let mut c2 = 1.0;
+            let mut ans2 = 1.0;
+            loop {
+                r2 += 1.0;
+                c2 *= x / r2;
+                ans2 += c2;
+                if (c2 / ans2) <= EPSILON {
+                    break;
+                }
+            }
+            return ax.exp() * ans2 / self;
+        }
+        let mut c = 0;
+        let mut y = 1.0 - self;
+        let mut z = x + y + 1.0;
+        let mut p3 = 1.0;
+        let mut q3 = x;
+        let mut p2 = x + 1.0;
+        let mut q2 = z * x;
+        let mut ans = p2 / q2;
+        let mut error;
+
+        loop {
+            c += 1;
+            y += 1.0;
+            z += 2.0;
+            let yc = y * (c as f64);
+            let p = p2 * z - p3 * yc;
+            let q = q2 * z - q3 * yc;
+            if q != 0.0 {
+                let nextans = p / q;
+                error = ((ans - nextans) / nextans).abs();
+                ans = nextans;
+            } else {
+                error = 1.0;
+            }
+            p3 = p2;
+            p2 = p;
+            q3 = q2;
+            q2 = q;
+            if p.abs() > BIG {
+                p3 *= BIG_INV;
+                p2 *= BIG_INV;
+                q3 *= BIG_INV;
+                q2 *= BIG_INV;
+            }
+            if error <= EPSILON {
+                break;
+            }
+        }
+        1.0 - ax.exp() * ans
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use crate::gamma::Gamma;
+    macro_rules! assert_almost_eq {
+        ($a:expr, $b:expr, $prec:expr) => {
+            #[allow(trivial_numeric_casts)]
+            if (($a - $b) as f64).abs() > $prec {
+                panic!(
+                    "assertion failed: `abs(left - right) < {:e}`, (left: `{}`, right: `{}`)",
+                    $prec, $a, $b
+                );
+            }
+        };
+    }
+
+    #[test]
+    fn test_gamma_lr() {
+        assert!(Gamma::lower_gamma_regularized(f64::NAN, f64::NAN).is_nan());
+        assert_almost_eq!(
+            Gamma::lower_gamma_regularized(0.1, 1.0),
+            0.97587265627367222115949155252812057714751052498477013,
+            1e-14
+        );
+        assert_eq!(
+            Gamma::lower_gamma_regularized(0.1, 2.0),
+            0.99432617602018847196075251078067514034772764693462125
+        );
+        assert_eq!(
+            Gamma::lower_gamma_regularized(0.1, 8.0),
+            0.99999507519205198048686442150578226823401842046310854
+        );
+        assert_almost_eq!(
+            Gamma::lower_gamma_regularized(1.5, 1.0),
+            0.42759329552912016600095238564127189392715996802703368,
+            1e-13
+        );
+        assert_almost_eq!(
+            Gamma::lower_gamma_regularized(1.5, 2.0),
+            0.73853587005088937779717792402407879809718939080920993,
+            1e-15
+        );
+        assert_eq!(
+            Gamma::lower_gamma_regularized(1.5, 8.0),
+            0.99886601571021467734329986257903021041757398191304284
+        );
+        assert_almost_eq!(
+            Gamma::lower_gamma_regularized(2.5, 1.0),
+            0.15085496391539036377410688601371365034788861473418704,
+            1e-13
+        );
+        assert_almost_eq!(
+            Gamma::lower_gamma_regularized(2.5, 2.0),
+            0.45058404864721976739416885516693969548484517509263197,
+            1e-14
+        );
+        assert_almost_eq!(
+            Gamma::lower_gamma_regularized(2.5, 8.0),
+            0.99315592607757956900093935107222761316136944145439676,
+            1e-15
+        );
+        assert_almost_eq!(
+            Gamma::lower_gamma_regularized(5.5, 1.0),
+            0.0015041182825838038421585211353488839717739161316985392,
+            1e-15
+        );
+        assert_almost_eq!(
+            Gamma::lower_gamma_regularized(5.5, 2.0),
+            0.030082976121226050615171484772387355162056796585883967,
+            1e-14
+        );
+        assert_almost_eq!(
+            Gamma::lower_gamma_regularized(5.5, 8.0),
+            0.85886911973294184646060071855669224657735916933487681,
+            1e-14
+        );
+        assert_almost_eq!(Gamma::lower_gamma_regularized(100.0, 0.5), 0.0, 1e-188);
+        assert_almost_eq!(Gamma::lower_gamma_regularized(100.0, 1.5), 0.0, 1e-141);
+        assert_almost_eq!(
+            Gamma::lower_gamma_regularized(100.0, 90.0),
+            0.1582209891864301681049696996709105316998233457433473,
+            1e-13
+        );
+        assert_almost_eq!(
+            Gamma::lower_gamma_regularized(100.0, 100.0),
+            0.5132987982791486648573142565640291634709251499279450,
+            1e-13
+        );
+        assert_almost_eq!(
+            Gamma::lower_gamma_regularized(100.0, 110.0),
+            0.8417213299399129061982996209829688531933500308658222,
+            1e-13
+        );
+        assert_almost_eq!(Gamma::lower_gamma_regularized(100.0, 200.0), 1.0, 1e-14);
+        assert_eq!(Gamma::lower_gamma_regularized(500.0, 0.5), 0.0);
+        assert_eq!(Gamma::lower_gamma_regularized(500.0, 1.5), 0.0);
+        assert_almost_eq!(Gamma::lower_gamma_regularized(500.0, 200.0), 0.0, 1e-70);
+        assert_almost_eq!(
+            Gamma::lower_gamma_regularized(500.0, 450.0),
+            0.0107172380912897415573958770655204965434869949241480,
+            1e-14
+        );
+        assert_almost_eq!(
+            Gamma::lower_gamma_regularized(500.0, 500.0),
+            0.5059471461707603580470479574412058032802735425634263,
+            1e-13
+        );
+        assert_almost_eq!(
+            Gamma::lower_gamma_regularized(500.0, 550.0),
+            0.9853855918737048059548470006900844665580616318702748,
+            1e-14
+        );
+        assert_almost_eq!(Gamma::lower_gamma_regularized(500.0, 700.0), 1.0, 1e-15);
+        assert_eq!(Gamma::lower_gamma_regularized(1000.0, 10000.0), 1.0);
+        assert_eq!(Gamma::lower_gamma_regularized(1e+50, 1e+48), 0.0);
+        assert_eq!(Gamma::lower_gamma_regularized(1e+50, 1e+52), 1.0);
+    }
+}

src/lib.rs

@@ -1,14 +1,4 @@
-pub fn add(left: u64, right: u64) -> u64 {
-    left + right
-}
+#![doc = include_str!(concat!(env!("CARGO_MANIFEST_DIR"), "/README.md"))]
+#![no_std]
 
-#[cfg(test)]
-mod tests {
-    use super::*;
-
-    #[test]
-    fn it_works() {
-        let result = add(2, 2);
-        assert_eq!(result, 4);
-    }
-}
+pub mod gamma;