generated from ProcyOS/rust-template
WIP: rewrite #7
14 changed files with 3490 additions and 443 deletions
1
.gitignore
vendored
1
.gitignore
vendored
|
|
@ -2,3 +2,4 @@
|
|||
/.direnv
|
||||
/result
|
||||
/resul-bin
|
||||
.vscode/
|
||||
|
|
@ -6,6 +6,11 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/),
|
|||
and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
|
||||
|
||||
## [Unreleased]
|
||||
### Removed
|
||||
Rewrite from scratch
|
||||
|
||||
### Added
|
||||
- a 256 bit integer type, `U256`.
|
||||
|
||||
## [0.1.1] - 2024-09-09
|
||||
|
||||
|
|
|
|||
28
Cargo.lock
generated
28
Cargo.lock
generated
|
|
@ -2,32 +2,6 @@
|
|||
# It is not intended for manual editing.
|
||||
version = 3
|
||||
|
||||
[[package]]
|
||||
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.1"
|
||||
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",
|
||||
]
|
||||
version = "0.2.0"
|
||||
|
|
|
|||
58
Cargo.nix
58
Cargo.nix
|
|
@ -22,7 +22,7 @@ args @ {
|
|||
workspaceSrc,
|
||||
ignoreLockHash,
|
||||
}: let
|
||||
nixifiedLockHash = "c2eb36fad35932f80d6eec4c083c16822de79c56cace0639af21e9933c2ca6f9";
|
||||
nixifiedLockHash = "74f9e44cc52b70b46b1a8e880aaaa90232d5db51c515bd0708c4b6732faf788b";
|
||||
workspaceSrc =
|
||||
if args.workspaceSrc == null
|
||||
then ./.
|
||||
|
|
@ -59,63 +59,17 @@ in
|
|||
in {
|
||||
cargo2nixVersion = "0.11.0";
|
||||
workspace = {
|
||||
extra-math = rustPackages.unknown.extra-math."0.1.1";
|
||||
extra-math = rustPackages.unknown.extra-math."0.2.0";
|
||||
};
|
||||
"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.1" = overridableMkRustCrate (profileName: rec {
|
||||
"unknown".extra-math."0.2.0" = overridableMkRustCrate (profileName: rec {
|
||||
name = "extra-math";
|
||||
version = "0.1.1";
|
||||
version = "0.2.0";
|
||||
registry = "unknown";
|
||||
src = fetchCrateLocal workspaceSrc;
|
||||
features = builtins.concatLists [
|
||||
(lib.optional (rootFeatures' ? "extra-math/default") "default")
|
||||
(lib.optional (rootFeatures' ? "extra-math/std") "std")
|
||||
(lib.optional (rootFeatures' ? "extra-math/u256") "u256")
|
||||
(lib.optional (rootFeatures' ? "extra-math/unstable") "unstable")
|
||||
];
|
||||
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;
|
||||
};
|
||||
});
|
||||
}
|
||||
|
|
|
|||
11
Cargo.toml
11
Cargo.toml
|
|
@ -1,23 +1,20 @@
|
|||
[package]
|
||||
name = "extra-math"
|
||||
version = "0.1.1"
|
||||
version = "0.2.0"
|
||||
edition = "2021"
|
||||
authors = ["Charlotte 🦝 Delenk <lotte@chir.rs>"]
|
||||
description = "Additional math functions not found in the standard library or libm"
|
||||
description = "Support for extra maths functions and floating point types"
|
||||
license = "MIT OR Apache-2.0"
|
||||
repository = "https://git.chir.rs/ProcyOS/rust-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"]
|
||||
unstable = []
|
||||
u256 = []
|
||||
|
||||
[lints.rust]
|
||||
deprecated-safe = "forbid"
|
||||
|
|
|
|||
|
|
@ -1,8 +1,12 @@
|
|||
# extra_math
|
||||
|
||||
Mathematical functions that are used in procyos’s code that aren’t found in either `libm` or the standard library.
|
||||
Maths support code for various floating point types
|
||||
|
||||
It is not intended to be accurate, just useful in a no-std environment.
|
||||
## Features
|
||||
|
||||
Each floating point type is feature gated. Some other functionality is enabled in some cases.
|
||||
|
||||
- `u256`: Enable a 256 bit unsigned integer type. This is *NOT* a constant time implementation.
|
||||
|
||||
## Running the tests
|
||||
|
||||
|
|
|
|||
138
build.rs
Normal file
138
build.rs
Normal file
|
|
@ -0,0 +1,138 @@
|
|||
//! Build script for compiler feature detection
|
||||
use std::env;
|
||||
|
||||
fn main() {
|
||||
#[cfg(all(feature = "unstable"))]
|
||||
{
|
||||
// Feature detection whether native f16/f128 support works
|
||||
// This is taken from https://github.com/rust-lang/rust/blob/master/library/std/build.rs
|
||||
println!("cargo:rerun-if-changed=build.rs");
|
||||
let target_arch =
|
||||
env::var("CARGO_CFG_TARGET_ARCH").expect("CARGO_CFG_TARGET_ARCH was not set");
|
||||
let target_os = env::var("CARGO_CFG_TARGET_OS").expect("CARGO_CFG_TARGET_OS was not set");
|
||||
let target_vendor =
|
||||
env::var("CARGO_CFG_TARGET_VENDOR").expect("CARGO_CFG_TARGET_VENDOR was not set");
|
||||
let target_pointer_width: u32 = env::var("CARGO_CFG_TARGET_POINTER_WIDTH")
|
||||
.expect("CARGO_CFG_TARGET_POINTER_WIDTH was not set")
|
||||
.parse()
|
||||
.unwrap();
|
||||
let is_miri = env::var_os("CARGO_CFG_MIRI").is_some();
|
||||
// Emit these on platforms that have no known ABI bugs, LLVM selection bugs, lowering bugs,
|
||||
// missing symbols, or other problems, to determine when tests get run.
|
||||
// If more broken platforms are found, please update the tracking issue at
|
||||
// <https://github.com/rust-lang/rust/issues/116909>
|
||||
//
|
||||
// Some of these match arms are redundant; the goal is to separate reasons that the type is
|
||||
// unreliable, even when multiple reasons might fail the same platform.
|
||||
println!("cargo:rustc-check-cfg=cfg(reliable_f16)");
|
||||
println!("cargo:rustc-check-cfg=cfg(reliable_f128)");
|
||||
|
||||
// This is a step beyond only having the types and basic functions available. Math functions
|
||||
// aren't consistently available or correct.
|
||||
println!("cargo:rustc-check-cfg=cfg(reliable_f16_math)");
|
||||
println!("cargo:rustc-check-cfg=cfg(reliable_f128_math)");
|
||||
|
||||
let has_reliable_f16 = match (target_arch.as_str(), target_os.as_str()) {
|
||||
// We can always enable these in Miri as that is not affected by codegen bugs.
|
||||
_ if is_miri => true,
|
||||
// Selection failure until recent LLVM <https://github.com/llvm/llvm-project/issues/93894>
|
||||
// FIXME(llvm19): can probably be removed at the version bump
|
||||
("loongarch64", _) => false,
|
||||
// Selection failure <https://github.com/llvm/llvm-project/issues/50374>
|
||||
("s390x", _) => false,
|
||||
// Unsupported <https://github.com/llvm/llvm-project/issues/94434>
|
||||
("arm64ec", _) => false,
|
||||
// MinGW ABI bugs <https://gcc.gnu.org/bugzilla/show_bug.cgi?id=115054>
|
||||
("x86_64", "windows") => false,
|
||||
// Apple has a special ABI for `f16` that we do not yet support
|
||||
// FIXME(builtins): fixed by <https://github.com/rust-lang/compiler-builtins/pull/675>
|
||||
("x86" | "x86_64", _) if target_vendor == "apple" => false,
|
||||
// Missing `__gnu_h2f_ieee` and `__gnu_f2h_ieee`
|
||||
("powerpc" | "powerpc64", _) => false,
|
||||
// Missing `__gnu_h2f_ieee` and `__gnu_f2h_ieee`
|
||||
("mips" | "mips32r6" | "mips64" | "mips64r6", _) => false,
|
||||
// Missing `__extendhfsf` and `__truncsfhf`
|
||||
("riscv32" | "riscv64", _) => false,
|
||||
// Most OSs are missing `__extendhfsf` and `__truncsfhf`
|
||||
(_, "linux" | "macos") => true,
|
||||
// Almost all OSs besides Linux and MacOS are missing symbols until compiler-builtins can
|
||||
// be updated. <https://github.com/rust-lang/rust/pull/125016> will get some of these, the
|
||||
// next CB update should get the rest.
|
||||
_ => false,
|
||||
};
|
||||
|
||||
let has_reliable_f128 = match (target_arch.as_str(), target_os.as_str()) {
|
||||
// We can always enable these in Miri as that is not affected by codegen bugs.
|
||||
_ if is_miri => true,
|
||||
// Unsupported <https://github.com/llvm/llvm-project/issues/94434>
|
||||
("arm64ec", _) => false,
|
||||
// ABI and precision bugs <https://github.com/rust-lang/rust/issues/125109>
|
||||
// <https://github.com/rust-lang/rust/issues/125102>
|
||||
("powerpc" | "powerpc64", _) => false,
|
||||
// Selection bug <https://github.com/llvm/llvm-project/issues/95471>
|
||||
("nvptx64", _) => false,
|
||||
// ABI unsupported <https://github.com/llvm/llvm-project/issues/41838>
|
||||
("sparc", _) => false,
|
||||
// MinGW ABI bugs <https://gcc.gnu.org/bugzilla/show_bug.cgi?id=115054>
|
||||
("x86_64", "windows") => false,
|
||||
// 64-bit Linux is about the only platform to have f128 symbols by default
|
||||
(_, "linux") if target_pointer_width == 64 => true,
|
||||
// Same as for f16, except MacOS is also missing f128 symbols.
|
||||
_ => false,
|
||||
};
|
||||
|
||||
// Configure platforms that have reliable basics but may have unreliable math.
|
||||
|
||||
// LLVM is currently adding missing routines, <https://github.com/llvm/llvm-project/issues/93566>
|
||||
let has_reliable_f16_math = has_reliable_f16
|
||||
&& match (target_arch.as_str(), target_os.as_str()) {
|
||||
// FIXME: Disabled on Miri as the intrinsics are not implemented yet.
|
||||
_ if is_miri => false,
|
||||
// x86 has a crash for `powi`: <https://github.com/llvm/llvm-project/issues/105747>
|
||||
("x86" | "x86_64", _) => false,
|
||||
// Assume that working `f16` means working `f16` math for most platforms, since
|
||||
// operations just go through `f32`.
|
||||
_ => true,
|
||||
};
|
||||
|
||||
let has_reliable_f128_math = has_reliable_f128
|
||||
&& match (target_arch.as_str(), target_os.as_str()) {
|
||||
// FIXME: Disabled on Miri as the intrinsics are not implemented yet.
|
||||
_ if is_miri => false,
|
||||
// LLVM lowers `fp128` math to `long double` symbols even on platforms where
|
||||
// `long double` is not IEEE binary128. See
|
||||
// <https://github.com/llvm/llvm-project/issues/44744>.
|
||||
//
|
||||
// This rules out anything that doesn't have `long double` = `binary128`; <= 32 bits
|
||||
// (ld is `f64`), anything other than Linux (Windows and MacOS use `f64`), and `x86`
|
||||
// (ld is 80-bit extended precision).
|
||||
("x86_64", _) => false,
|
||||
(_, "linux") if target_pointer_width == 64 => true,
|
||||
_ => false,
|
||||
};
|
||||
|
||||
if has_reliable_f16 {
|
||||
println!("cargo:rustc-cfg=reliable_f16");
|
||||
}
|
||||
if has_reliable_f128 {
|
||||
println!("cargo:rustc-cfg=reliable_f128");
|
||||
}
|
||||
if has_reliable_f16_math {
|
||||
println!("cargo:rustc-cfg=reliable_f16_math");
|
||||
}
|
||||
if has_reliable_f128_math {
|
||||
println!("cargo:rustc-cfg=reliable_f128_math");
|
||||
}
|
||||
}
|
||||
|
||||
println!("cargo:rustc-check-cfg=cfg(overflow_checks_stable)");
|
||||
|
||||
if std::panic::catch_unwind(|| {
|
||||
#[allow(arithmetic_overflow)]
|
||||
let _ = 255_u8 + 1;
|
||||
})
|
||||
.is_err()
|
||||
{
|
||||
println!("cargo:rustc-cfg=overflow_checks_stable");
|
||||
}
|
||||
}
|
||||
11
flake.lock
11
flake.lock
|
|
@ -14,16 +14,15 @@
|
|||
]
|
||||
},
|
||||
"locked": {
|
||||
"lastModified": 1713199118,
|
||||
"narHash": "sha256-MlLdAvk+zXCFUy280sY6LqtykqWXIkKVXo72J7a6HlU=",
|
||||
"owner": "cargo2nix",
|
||||
"lastModified": 1726552619,
|
||||
"narHash": "sha256-ytTBILVMnRZYvjiLYz+J6IFf/TOXdGuP6RDesMx9qgA=",
|
||||
"owner": "DarkKirb",
|
||||
"repo": "cargo2nix",
|
||||
"rev": "1efb03f2f794ad5eed17e807e858c4da001dbc3e",
|
||||
"rev": "baa12124e2de09e1cbbdac320f14809fa55af1a2",
|
||||
"type": "github"
|
||||
},
|
||||
"original": {
|
||||
"owner": "cargo2nix",
|
||||
"ref": "main",
|
||||
"owner": "DarkKirb",
|
||||
"repo": "cargo2nix",
|
||||
"type": "github"
|
||||
}
|
||||
|
|
|
|||
|
|
@ -2,17 +2,17 @@
|
|||
description = "extra-math";
|
||||
|
||||
inputs = {
|
||||
nixpkgs.url = github:NixOS/nixpkgs;
|
||||
flake-utils.url = github:numtide/flake-utils;
|
||||
nixpkgs.url = "github:NixOS/nixpkgs";
|
||||
flake-utils.url = "github:numtide/flake-utils";
|
||||
|
||||
rust-overlay = {
|
||||
url = github:oxalica/rust-overlay;
|
||||
url = "github:oxalica/rust-overlay";
|
||||
inputs.nixpkgs.follows = "nixpkgs";
|
||||
inputs.flake-utils.follows = "flake-utils";
|
||||
};
|
||||
|
||||
cargo2nix = {
|
||||
url = github:cargo2nix/cargo2nix/main;
|
||||
url = "github:DarkKirb/cargo2nix";
|
||||
inputs.nixpkgs.follows = "nixpkgs";
|
||||
inputs.flake-utils.follows = "flake-utils";
|
||||
inputs.rust-overlay.follows = "rust-overlay";
|
||||
|
|
|
|||
67
src/gamma.rs
67
src/gamma.rs
|
|
@ -1,67 +0,0 @@
|
|||
//! 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);
|
||||
}
|
||||
}
|
||||
239
src/gamma/f64.rs
239
src/gamma/f64.rs
|
|
@ -1,239 +0,0 @@
|
|||
//! 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);
|
||||
}
|
||||
}
|
||||
|
|
@ -1,5 +1,5 @@
|
|||
#![doc = include_str!(concat!(env!("CARGO_MANIFEST_DIR"), "/README.md"))]
|
||||
#![no_std]
|
||||
|
||||
pub mod gamma;
|
||||
pub mod n_choose_k;
|
||||
#[cfg(feature = "u256")]
|
||||
pub mod u256;
|
||||
|
|
|
|||
|
|
@ -1,36 +0,0 @@
|
|||
//! Calculates the binomial coefficient N-choose-K.
|
||||
|
||||
use num_traits::{Num, NumCast, PrimInt};
|
||||
|
||||
/// Calculates the binomial coefficient
|
||||
pub trait BinomialCoefficient: Num + NumCast {
|
||||
/// Calculates the binary coefficient N-choose-K.
|
||||
///
|
||||
/// # Return
|
||||
/// Returns None if k cannot be represented in the current number type.
|
||||
fn choose<I: PrimInt>(self, k: I) -> Option<Self>;
|
||||
}
|
||||
|
||||
impl<N: Num + NumCast + Copy> BinomialCoefficient for N {
|
||||
fn choose<I: PrimInt>(self, k: I) -> Option<Self> {
|
||||
let mut num = Self::one();
|
||||
let mut denom = Self::one();
|
||||
let mut i = I::one();
|
||||
while i <= k {
|
||||
num = num * (self + Self::one() - Self::from(i)?);
|
||||
denom = denom * Self::from(i)?;
|
||||
i = i + I::one();
|
||||
}
|
||||
Some(num / denom)
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod tests {
|
||||
use crate::n_choose_k::BinomialCoefficient;
|
||||
|
||||
#[test]
|
||||
fn test_binomial() {
|
||||
assert_eq!(BinomialCoefficient::choose(4, 2), Some(6));
|
||||
}
|
||||
}
|
||||
3317
src/u256/mod.rs
Normal file
3317
src/u256/mod.rs
Normal file
File diff suppressed because it is too large
Load diff
Loading…
Reference in a new issue