diff --git a/CHANGELOG.md b/CHANGELOG.md index 5046eab..9229950 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 4545d55..cba31ce 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -1,4 +1,4 @@ -# Contributing to [Project] +# Contributing to `extra-math` Thank you for considering to contribute to our project! diff --git a/Cargo.lock b/Cargo.lock index be35886..2372726 100644 --- a/Cargo.lock +++ b/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", +] diff --git a/Cargo.nix b/Cargo.nix index b600c1c..713115c 100644 --- a/Cargo.nix +++ b/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; + }; }); } diff --git a/Cargo.toml b/Cargo.toml index f0a7390..5bd09aa 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,15 +1,23 @@ [package] -name = "rust-template" +name = "extra-math" version = "0.1.0" edition = "2021" authors = ["Charlotte 🦝 Delenk "] -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" diff --git a/README.md b/README.md index af24f77..3c36783 100644 --- a/README.md +++ b/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 diff --git a/flake.nix b/flake.nix index 7d3c0b2..708a530 100644 --- a/flake.nix +++ b/flake.nix @@ -1,5 +1,5 @@ { - description = "rust-template"; + description = "extra-math"; inputs = { nixpkgs.url = github:NixOS/nixpkgs; diff --git a/src/gamma.rs b/src/gamma.rs new file mode 100644 index 0000000..61eb638 --- /dev/null +++ b/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); + } +} diff --git a/src/gamma/f64.rs b/src/gamma/f64.rs new file mode 100644 index 0000000..073d0f7 --- /dev/null +++ b/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); + } +} diff --git a/src/lib.rs b/src/lib.rs index b93cf3f..3d961d2 100644 --- a/src/lib.rs +++ b/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;