generated from ProcyOS/rust-template
WIP: rewrite #7
14 changed files with 3490 additions and 443 deletions
3
.gitignore
vendored
3
.gitignore
vendored
|
|
@ -1,4 +1,5 @@
|
||||||
/target
|
/target
|
||||||
/.direnv
|
/.direnv
|
||||||
/result
|
/result
|
||||||
/resul-bin
|
/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).
|
and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
|
||||||
|
|
||||||
## [Unreleased]
|
## [Unreleased]
|
||||||
|
### Removed
|
||||||
|
Rewrite from scratch
|
||||||
|
|
||||||
|
### Added
|
||||||
|
- a 256 bit integer type, `U256`.
|
||||||
|
|
||||||
## [0.1.1] - 2024-09-09
|
## [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.
|
# It is not intended for manual editing.
|
||||||
version = 3
|
version = 3
|
||||||
|
|
||||||
[[package]]
|
|
||||||
name = "autocfg"
|
|
||||||
version = "1.3.0"
|
|
||||||
source = "registry+https://github.com/rust-lang/crates.io-index"
|
|
||||||
checksum = "0c4b4d0bd25bd0b74681c0ad21497610ce1b7c91b1022cd21c80c6fbdd9476b0"
|
|
||||||
|
|
||||||
[[package]]
|
[[package]]
|
||||||
name = "extra-math"
|
name = "extra-math"
|
||||||
version = "0.1.1"
|
version = "0.2.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",
|
|
||||||
]
|
|
||||||
|
|
|
||||||
58
Cargo.nix
58
Cargo.nix
|
|
@ -22,7 +22,7 @@ args @ {
|
||||||
workspaceSrc,
|
workspaceSrc,
|
||||||
ignoreLockHash,
|
ignoreLockHash,
|
||||||
}: let
|
}: let
|
||||||
nixifiedLockHash = "c2eb36fad35932f80d6eec4c083c16822de79c56cace0639af21e9933c2ca6f9";
|
nixifiedLockHash = "74f9e44cc52b70b46b1a8e880aaaa90232d5db51c515bd0708c4b6732faf788b";
|
||||||
workspaceSrc =
|
workspaceSrc =
|
||||||
if args.workspaceSrc == null
|
if args.workspaceSrc == null
|
||||||
then ./.
|
then ./.
|
||||||
|
|
@ -59,63 +59,17 @@ in
|
||||||
in {
|
in {
|
||||||
cargo2nixVersion = "0.11.0";
|
cargo2nixVersion = "0.11.0";
|
||||||
workspace = {
|
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 {
|
"unknown".extra-math."0.2.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 {
|
|
||||||
name = "extra-math";
|
name = "extra-math";
|
||||||
version = "0.1.1";
|
version = "0.2.0";
|
||||||
registry = "unknown";
|
registry = "unknown";
|
||||||
src = fetchCrateLocal workspaceSrc;
|
src = fetchCrateLocal workspaceSrc;
|
||||||
features = builtins.concatLists [
|
features = builtins.concatLists [
|
||||||
(lib.optional (rootFeatures' ? "extra-math/default") "default")
|
(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]
|
[package]
|
||||||
name = "extra-math"
|
name = "extra-math"
|
||||||
version = "0.1.1"
|
version = "0.2.0"
|
||||||
edition = "2021"
|
edition = "2021"
|
||||||
authors = ["Charlotte 🦝 Delenk <lotte@chir.rs>"]
|
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"
|
license = "MIT OR Apache-2.0"
|
||||||
repository = "https://git.chir.rs/ProcyOS/rust-template"
|
repository = "https://git.chir.rs/ProcyOS/rust-template"
|
||||||
keywords = ["mathematics"]
|
keywords = ["mathematics"]
|
||||||
categories = ["mathematics", "no-std", "no-std::no-alloc"]
|
categories = ["mathematics", "no-std", "no-std::no-alloc"]
|
||||||
|
|
||||||
[dependencies]
|
[dependencies]
|
||||||
libm = "0.2.8"
|
|
||||||
num-traits = { version = "0.2.19", default-features = false, features = [
|
|
||||||
"libm",
|
|
||||||
] }
|
|
||||||
|
|
||||||
[features]
|
[features]
|
||||||
default = []
|
default = []
|
||||||
std = ["num-traits/std"]
|
unstable = []
|
||||||
|
u256 = []
|
||||||
|
|
||||||
[lints.rust]
|
[lints.rust]
|
||||||
deprecated-safe = "forbid"
|
deprecated-safe = "forbid"
|
||||||
|
|
|
||||||
|
|
@ -1,8 +1,12 @@
|
||||||
# extra_math
|
# 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
|
## 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": {
|
"locked": {
|
||||||
"lastModified": 1713199118,
|
"lastModified": 1726552619,
|
||||||
"narHash": "sha256-MlLdAvk+zXCFUy280sY6LqtykqWXIkKVXo72J7a6HlU=",
|
"narHash": "sha256-ytTBILVMnRZYvjiLYz+J6IFf/TOXdGuP6RDesMx9qgA=",
|
||||||
"owner": "cargo2nix",
|
"owner": "DarkKirb",
|
||||||
"repo": "cargo2nix",
|
"repo": "cargo2nix",
|
||||||
"rev": "1efb03f2f794ad5eed17e807e858c4da001dbc3e",
|
"rev": "baa12124e2de09e1cbbdac320f14809fa55af1a2",
|
||||||
"type": "github"
|
"type": "github"
|
||||||
},
|
},
|
||||||
"original": {
|
"original": {
|
||||||
"owner": "cargo2nix",
|
"owner": "DarkKirb",
|
||||||
"ref": "main",
|
|
||||||
"repo": "cargo2nix",
|
"repo": "cargo2nix",
|
||||||
"type": "github"
|
"type": "github"
|
||||||
}
|
}
|
||||||
|
|
|
||||||
|
|
@ -2,17 +2,17 @@
|
||||||
description = "extra-math";
|
description = "extra-math";
|
||||||
|
|
||||||
inputs = {
|
inputs = {
|
||||||
nixpkgs.url = github:NixOS/nixpkgs;
|
nixpkgs.url = "github:NixOS/nixpkgs";
|
||||||
flake-utils.url = github:numtide/flake-utils;
|
flake-utils.url = "github:numtide/flake-utils";
|
||||||
|
|
||||||
rust-overlay = {
|
rust-overlay = {
|
||||||
url = github:oxalica/rust-overlay;
|
url = "github:oxalica/rust-overlay";
|
||||||
inputs.nixpkgs.follows = "nixpkgs";
|
inputs.nixpkgs.follows = "nixpkgs";
|
||||||
inputs.flake-utils.follows = "flake-utils";
|
inputs.flake-utils.follows = "flake-utils";
|
||||||
};
|
};
|
||||||
|
|
||||||
cargo2nix = {
|
cargo2nix = {
|
||||||
url = github:cargo2nix/cargo2nix/main;
|
url = "github:DarkKirb/cargo2nix";
|
||||||
inputs.nixpkgs.follows = "nixpkgs";
|
inputs.nixpkgs.follows = "nixpkgs";
|
||||||
inputs.flake-utils.follows = "flake-utils";
|
inputs.flake-utils.follows = "flake-utils";
|
||||||
inputs.rust-overlay.follows = "rust-overlay";
|
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"))]
|
#![doc = include_str!(concat!(env!("CARGO_MANIFEST_DIR"), "/README.md"))]
|
||||||
#![no_std]
|
#![no_std]
|
||||||
|
|
||||||
pub mod gamma;
|
#[cfg(feature = "u256")]
|
||||||
pub mod n_choose_k;
|
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