Chapter 1: What is NumeriLib?

1.1 Overview:

Welcome to the official book of the NumeriLib Mathematics Library. This library is designed to provide a comprehensive set of mathematical functions for Rust programmers. Whether you are working on scientific computing, data analysis, or any other math-intensive project, this library aims to be your one-stop shop.

1.2 Work in Progress:

As a testament to my commitment to continuous improvement, NumeriLib is an ongoing project. I am constantly expanding and enhancing its capabilities, striving to provide a robust and comprehensive set of mathematical tools. Your feedback and contributions are invaluable in shaping the future of the library.

1.3 Getting Started:

1.3.1 Installation:

To install NumeriLib, simply add the following line to your Cargo.toml file:

...

[dependencies]
numerilib = "0.1.0"

...

Or, you can add it via the command line:

cargo add numerilib

Finally, if you want to use the latest version from the Git Repository, add the following line to your Cargo.toml file:

...

[dependencies]
numerilib = { git = "https://github.com/VLambda/NumeriLib.git", branch = "main" } 

...

Chapter 2: Basic Functions

Calculus/Basic Functions

  • Some Basic Functions in Mathematica

Derivatives:

Calculates the Derivative of a function at a given point: using the limit definition of a derivative. The Derivative of a function is the slope of the tangent line to the function at a given point.

Parameters:

  • func: A function that takes a single f64 argument and returns an f64. This is the function for which the derivative will be calculated.
  • x: The point at which the derivative will be calculated.

Returns:

The calculated derivative of the function at the given point.

Equation:

\[ f'\left( x \right)=\lim_{h \to 0} \frac{f\left( x+h \right)-f\left( x \right)}{h} \]

Example:

use mathematica::Functions;

fn main() {
    let function = |x: f64| x.powi(2);
    let x = 2;

    let derivative = Functions::derivative(function, x);

    println!("The Derivative of x^2 as x=2 is: {}", derivative);
}                     

Execute as:

>>> The Derivative of x^2 as x=2 is: 4

Summation:

Calculates Summations of a function. The Summation of a function is the sum of all the values of the function from a given start to a given limit.

Parameters:

  • start: The starting value of the summation.
  • limit: The ending value of the summation.
  • func: A function that takes a single f64 argument and returns an f64. This is the function for which the summation will be calculated.

Returns:

The summation of the function from the given start to the given limit.

Equation:

\[ \sum_{i=m}^{n}a_{i}=a_{m}+a_{m+1}+a_{m+2}+\cdot \cdot \cdot +a_{n-1}+a_{n} \]

Example #1: Constant

use ferrate::special::Functions;                             
                                                             
let start = 0;                                               
let limit = 9;                                               
let function = |x: f64| 3_f64;                               
                                                             
let summation = Functions::summation(start, limit, function);
                                                             
println!("The Summation of the constant, 3 from 0 to 9 is: {}", summation);

Execute as:

>>> The Summation of the constant, 3 from 0 to 9 is: 30

Example #2: Function

use ferrate::special::Functions;                             
                                                             
let start = 4.5;                                             
let limit = 100;                                             
let function = |x: f64| 1_f64 / x;                           
                                                             
let summation = Functions::summation(start, limit, function);
                                                             
println!("The Summation of 1/x from 4.5 to 100 is: {}", summation);

Execute as:

>>> The Summation of 1/x from 4.5 to 100 is: 3.104044184318839

Product:

Calculates the Product (a.k.a. Capital Pi Notation) of a function. The Product of a function is the product of all the values of the function from a given start to a given limit.

Parameters:

  • start: The starting value of the product.
  • limit: The ending value of the product.
  • func: A function that takes a single f64 argument and returns an f64. This is the function for which the product will be calculated.

Returns:

The product of the function from the given start to the given limit.

Equation:

\[ \prod_{i=m}^{n}a_{i}=a_{m}\times a_{m+1}\times a_{m+2}\times \cdot \cdot \cdot \times a_{n-1}\times a_{n} \]

Example #1: Constant

use ferrate::special::Functions;                             

fn main() {                                                             
    let start = 2_f64;                                           
    let limit = 7_f64;                                           
    let f = |x: f64| 3_f64;                                      
                                                                 
    let product_series = Functions::product(start, limit, f);    
    
    println!("The Product of the constant, 3 from 2 to 7 is: {}", product_series);
}

Execute as:

>>> The Product of the constant, 3 from 2 to 7 is: 2187

Example #2: Function

use ferrate::special::Functions;                             
     
fn main() {                                                        
    let start = 1_f64;                                           
    let limit = 6_f64;                                           
    let f = |x: f64| x.powi(2);                                  
                                                                 
    let product_series = Functions::product(start, limit, f);    
    
    println!("The Product of x^2 from 1 to 6 is: {}", product_series);
}

Execute as:

>>> The Product of x^2 from 1 to 6 is: 518400

Newton's Method:

Calculates the root of a function using Newton's Method. Newton's Method is an iterative method for finding the roots of a function. It is based on the idea that a continuous and differentiable function can be approximated by a straight line tangent to it.

Parameters:

  • x: The initial guess for the root of the function.
  • func: A function that takes a single f64 argument and returns an f64. This is the function for which the root will be calculated.

Returns:

The root of the function.

Equation:

\[ x_{n+1}=x_{n}-\frac{f\left( x_{n} \right)}{f'\left( x_{n} \right)} \]

Example:

use ferrate::special::Functions;             

fn main() {                                   
    let x = 1.5_f64;                             
    let function = |x: f64| x.powi(2) - 2_f64;   
    let newton = Functions::newmet(x, function); 
                                                 
    println!("The Root of x^2 - 2 is: {}", newton);
}                     

Execute as:

>>> The Root of x^2 - 2 is: 1.414213562373095

Definite Integration

  • Functions for Definite Integration in Mathematica

Right Endpoint Riemann Sum:

The Right Endpoint method to calculate a definite integral. The Riemann Sum is an approximation of the area under a curve using rectangles. The Right Endpoint method uses the right endpoint of each rectangle (interval) to calculate the height of the rectangle. The height of the rectangle is then multiplied by the width of the rectangle to calculate the area of the rectangle. The area of each rectangle is then summed to calculate the area under the curve.

Parameters:

  • function: A function that takes a single f64 argument and returns an f64. This is the function to be integrated.
  • lower_limit: The lower limit of integration.
  • upper_limit: The upper limit of integration.
  • intervals: The number of the intervals for the Riemann sum.

Returns:

The calculated definite integral using the Right Endpoint method.

Equation:

\[ \int_{a}^{b}f\left( x \right)dx\approx \Delta x\sum_{i=1}^{n}f\left( a + x_{i}\Delta x \right) \] \[ \Delta x=\frac{b-a}{\underbrace{n}_\text{Intervals}} \]

Example:

use mathematica::Functions;

fn main() {
    let lower_bound = 0_f64;
    let upper_bound = 6_f64;
    let intervals = 6_f64;
    let function = |x: f64| x.powi(2);

    let integral = Functions::right_riemann(function, lower_bound, upper_bound, intervals);

    println!("The Integral of x^2 at [0,6] using Right Endpoints is: {}", integral)
}

Execute as:

>>> The Integral of x^2 at [0,6] using Right Endpoints is: 91

Left Endpoint Riemann Sum:

The Left Endpoint method to calculate a definite integral. The Riemann Sum is an approximation of the area under a curve using rectangles. The Left Endpoint method uses the left endpoint of each rectangle (interval) to calculate the height of the rectangle. The height of the rectangle is then multiplied by the width of the rectangle to calculate the area of the rectangle. The area of each rectangle is then summed to calculate the area under the curve.

Parameters:

  • function: A function that takes a single f64 argument and returns an f64. This is the function to be integrated.
  • lower_limit: The lower limit of integration.
  • upper_limit: The upper limit of integration.
  • intervals: The number of the intervals for the Riemann sum.

Returns:

The calculated definite integral using the Left Endpoint method.

Equation:

\[ \int_{a}^{b}f\left( x \right)dx\approx \Delta x\sum_{i=0}^{n-1}f\left( a + x_{i}\Delta x \right) \] \[ \Delta x=\frac{b-a}{\underbrace{n}_\text{Intervals}} \]

Example:

use mathematica::Functions;

fn main() {
    let lower_bound = 0_f64;
    let upper_bound = 6_f64;
    let intervals = 6_f64;
    let function = |x: f64| x.powi(2);

    let integral = Functions::left_riemann(function, lower_bound, upper_bound, intervals);

    println!("The Integral of x^2 at [0,6] using Left Endpoints is: {}", integral)
}

Execute as:

>>> The Integral of x^2 at [0,6] using Left Endpoints is: 55

Midpoint Riemann Sum:

The Midpoint method to calculate a definite integral. The Riemann Sum is an approximation of the area under a curve using rectangles. The Midpoint method uses the midpoint of each rectangle (interval) to calculate the height of the rectangle. The height of the rectangle is then multiplied by the width of the rectangle to calculate the area of the rectangle. The area of each rectangle is then summed to calculate the area under the curve.

Parameters

  • function: A function that takes a single f64 argument and returns an f64. This is the function to be integrated.
  • lower_limit: The lower limit of integration.
  • upper_limit: The upper limit of integration.
  • intervals: The number of the intervals for the Riemann sum.

Returns

The calculated definite integral using the Midpoint method.

Equation:

\[ \int_{a}^{b}f\left( x \right)dx\approx \Delta x\sum_{i=0}^{n}f\left( a + x_{i + \frac{1}{2}}\Delta x^{2} \right) \] \[ \Delta x=\frac{b-a}{\underbrace{n}_\text{Intervals}} \]

Example

use mathematica::Functions;

fn main() {
    let lower_bound = 0_f64;
    let upper_bound = 6_f64;
    let intervals = 6_f64;
    let function = |x: f64| x.powi(2);

    let integral = Functions::midpoint_riemann(function, lower_bound, upper_bound, intervals);

    println!("The Integral of x^2 at [0,6] using Midpoints is: {}", integral)
}

Execute as:

>>> The Integral of x^2 at [0,6] using Midpoints is: 71.5

Trapezoid Rule:

The Trapezoid method to calculate a definite integral. The Trapezoid Rule is an approximation of the area under a curve using trapezoids. The Trapezoid Rule uses the left and right endpoints of each interval to calculate the height of the trapezoid. The height of the trapezoid is then multiplied by the width of the trapezoid to calculate the area of the trapezoid. The area of each trapezoid is then summed to calculate the area under the curve.

Parameters

  • function: A function that takes a single f64 argument and returns an f64. This is the function to be integrated.
  • lower_limit: The lower limit of integration.
  • upper_limit: The upper limit of integration.
  • intervals: The number of the intervals for the Riemann sum.

Returns

The calculated definite integral using the Trapezoidal rule.

Equation:

\[ \int_{a}^{b}f\left( x \right)dx \approx \frac{1}{2}\Delta x\sum_{i=0}^{n-1}f\left( a + x_{i} \Delta x \right) + f\left( a + x_{i+1} \Delta x \right) \] \[ \Delta x=\frac{b-a}{\underbrace{n}_\text{Intervals}} \]

Example

use mathematica::Functions;

fn main() {
    let lower_bound = 0_f64;
    let upper_bound = 6_f64;
    let intervals = 6_f64;
    let function = |x: f64| x.powi(2);

    let integral = Functions::trapezoid(function, lower_bound, upper_bound, intervals);

    println!("The Integral of x^2 at [0,6] using the Trapezoid Rule is: {}", integral)
}

Execute as:

>>> The Integral of x^2 at [0,6] using the Trapezoid Rule is: 73

Simpson's Rule:

Uses the Composite Simpson's 1/3rd Rule to calculate a definite integral. Simpson's Rule is an approximation of the area under a curve using parabolas. Simpson's Rule uses the left and right endpoints of each interval to calculate the height of the parabola. The height of the parabola is then multiplied by the width of the parabola to calculate the area of the parabola. The area of each parabola is then summed to calculate the area under the curve.

Parameters:

  • function: A function that takes a single f64 argument and returns an f64. This is the function to be integrated.
  • lower_limit: The lower limit of integration.
  • upper_limit: The upper limit of integration.
  • intervals: The number of the intervals for the Riemann sum.

Returns:

The calculated definite integral using Simpson's Rule.

Equation:

\[ \int_{a}^{b}f\left( x \right)dx \approx \frac{1}{3} \Delta x \left[ f\left( a \right) + 4 \sum_{i=0}^{\frac{n}{2}} f\left(a + x_{2i-1} \frac{\Delta x}{2} \right) + 2 \sum_{i=1}^{\frac{n}{2} - 1} f\left( a + x_{2i} \frac{\Delta x}{2} \right) + f\left( b \right) \right] \] \[ \Delta x=\frac{b-a}{\underbrace{n}_\text{Intervals}} \]

Example:

use mathematica::Functions;

fn main() {
    let lower_bound = 0_f64;
    let upper_bound = 6_f64;
    let intervals = 6_f64;
    let function = |x: f64| x.powi(2);

    let integral = Functions::simpson(function, lower_bound, upper_bound, intervals);

    println!("The Integral of x^2 at [0,6] using Simpson's Rule is: {}", integral)
}

Execute as:

>>> The Integral of x^2 at [0,6] using Simpson's Rule is: 72

Boole's Rule:

Uses Boole's Rule to calculate a definite integral. Boole's Rule is an approximation of the area under a curve using quartic polynomials. Boole's Rule uses the left and right endpoints of each interval to calculate the height of the quartic polynomial. The height of the quartic polynomial is then multiplied by the width of the quartic polynomial to calculate the area of the quartic polynomial. The area of each quartic polynomial is then summed to calculate the area under the curve.

Parameters:

  • function: A function that takes a single f64 argument and returns an f64. This is the function to be integrated.
  • lower_limit: The lower limit of integration.
  • upper_limit: The upper limit of integration.

Returns:

The calculated definite integral using Boole's Rule.

Equation:

\[ \int_{a}^{b}f\left( x \right)dx \to \int_{x_{0}}^{x_{4}}f\left( x \right)dx \]

\[ \Delta x = \frac{b-a}{4} \]

Intervals
\( x_{0} = a \)
\( x_{1} = a + {\Delta x} \)
\( x_{2} = a + {2 \Delta x} \)
\( x_{3} = a + {3 \Delta x} \)
\( x_{4} = b \)

\[ \int_{x_{0}}^{x_{4}}f\left( x \right)dx\approx \frac{2\Delta x}{45}\left[ 7f\left( x_{0} \right) + 32f\left( x_{1} \right) + 12f\left( x_2 \right) + 32f\left( x_{3} \right) + 7f\left( x_{4} \right) \right] \]

Example:

use numerilib::Functions;

fn main() {
    let lower_bound = 0_f64;
    let upper_bound = 6_f64;
    let function = |x: f64| x.powi(2);
    
    let integral = Functions::boole_rule(function, lower_bound, upper_bound);
    
    println!("The Integral of x^2 at [0,6] using Boole's Rule is: {}", integral)
}

Execute as:

>>> The Integral of x^2 at [0,6] using Boole's Rule is: 72

Adaptive Quadrature:

Uses Adaptive Quadrature to calculate a definite integral. Adaptive Quadrature is an approximation of the area under a curve using Simpson's Rule. Adaptive Quadrature uses Simpson's Rule to calculate the area under the curve. If the error is greater than the tolerance, the interval is split in half and the area under the curve is calculated again. This process is repeated until the error is less than the tolerance.

Parameters:

  • function: A function that takes a single f64 argument and returns an f64. This is the function to be integrated.
  • lower_limit: The lower limit of integration.
  • upper_limit: The upper limit of integration.
  • tolerance: The level of precision (ie: 1e-6) that is passed.

Returns:

The calculated definite integral using Adaptive Quadrature.

Procedure:

Adaptive Quadrature is Described as the following:

1. Integrate the function using Simpson's Rule.
   I = ∫[a,b]f(x)dx
   c = (a + b) / 2
   Q = ∫[a,c]f(x)dx 
   R = ∫[c,b]f(x)dx
   I ≈ Q + R

2. If the error is greater than the tolerance
   split the interval in half and integrate again.
   
   ε = |Q - R|
   If ε > tolerance, then:
   
    a. Split the intervals in half.
       d = (a + c) / 2
       e = (c + b) / 2
   
    b. Integrate the function using Simpson's Rule.
         M1 = ∫[a,d]f(x)dx + ∫[d,b]f(x)dx
         M2 = ∫[a,e]f(x)dx + ∫[e,b]f(x)dx
         I ≈ M1 + M2
   
3. Repeat until the error is less than the tolerance.
   In which, return the value.

Example:

use numerilib::Functions;

fn main() {
    let lower_bound = 0_f64;
    let upper_bound = 6_f64;
    let tolerance = 1e-12;
    let function = |x: f64| x.powi(2);
    
    let integral = Functions::adaptive_quadrature(function, lower_bound, upper_bound, tolerance);
    
    println!("The Integral of x^2 at [0,6] using Adaptive Quadrature is: {}", integral)
}

Execute as:

>>> The Integral of x^2 at [0,6] using Adaptive Quadrature is: 72

Complex Numbers & Functions

  • Supports Complex Numbers and Functions relating to Complex Numbers

Complex Numbers:

Complex Numbers are numbers that have both a real and imaginary part i. They are represented as a + bi where a is the real part and b is the imaginary part. i is the imaginary unit, where i^2 = -1.

Complex Numbers are supported by the struct Complex in the numerlib::Complex module. To create a new Complex Number, use the new() function.

Parameters:

  • real: The real part of the Complex Number, represented by a f64 value.
  • imag: The imaginary part of the Complex Number, represented by a f64 value.

Returns:

A new Complex Number, with the complex type.

Equation:

\[ a + bi \]

Example:

use numerilib::Complex;

fn main() {
    let z = Complex::new(1.0, 2.0);

    println!("Complex Number z: {}", z);
}

Executes as:

>>> Complex Number z: 1+2i

Returning Parts of a Complex Number:

Inorder to return the real or imaginary part of a Complex Number, use the .real_part() or imag_part() functions.

Parameters:

  • &self: The Complex Number to return the real or imaginary part of.

Returns:

The real or imaginary part of the Complex Number, represented by a f64 value.

Returning Real Value Example:

use numerilib::Complex;

fn main() {
    let z = Complex::new(1.0, 2.0);

    let real = z.real_part();
    
    println!("Real Part of z: {}", real);
}

Executes as:

>>> Real Part of z: 1

Returning Imaginary Value Example:

use numerilib::Complex;

fn main() {
    let z = Complex::new(1.0, 2.0);

    let imag = z.imag_part();
    
    println!("Imaginary Part of z: {}", imag);
}

Executes as:

>>> Imaginary Part of z: 2

The Magnitude of a Complex Number:

The Magnitude of a Complex Number is the distance from the origin to the Complex Number. It is represented by the absolute value of the Complex Number.

Parameters:

  • &self: The Complex Number to return the magnitude of.

Returns:

The magnitude of the Complex Number, represented by a f64 value.

Equation:

\[ \left| z \right| = \sqrt{a^2 + b^2} \]

Example:

use numerilib::Complex;

fn main() {
    let z = Complex::new(1.0, 2.0);

    let magnitude = z.magnitude();
    
    println!("Magnitude of z: {}", magnitude);
}

Executes as:

>>> Magnitude of z: 2.23606797749979 // √5

Complex Arithmetic:

Complex Arithmetic is supported by the Complex struct in the numerilib::Complex module. One Quick Notice: Complex Arithmetic can only be done by using operations in the Complex type.

Addition Example:

use numerilib::Complex;

fn main() {
    let a = Complex::new(3.0, 2.0);
    let b = Complex::new(1.0, 4.0);

    let sum = a + b;
    
    println!("Sum of 3+2i + 1+4i = {}", sum);
}

Executes as:

>>> Sum of 3+2i + 1+4i = 4+6i

Subtraction Example:

use numerilib::Complex;

fn main() {
    let a = Complex::new(3.0, 2.0);
    let b = Complex::new(1.0, 4.0);

    let difference = a - b;
    
    println!("Difference of 3+2i - 1+4i = {}", difference);
}

Executes as:

>>> Difference of 3+2i - 1+4i = 2-2i

Multiplication Example:

use numerilib::Complex;

fn main() {
    let a = Complex::new(3.0, 2.0);
    let b = Complex::new(1.0, 4.0);

    let product = a * b;
    
    println!("Product of 3+2i * 1+4i = {}", product);
}

Executes as:

>>> Product of 3+2i * 1+4i = -5+14i

Division Example:

use numerilib::Complex;

fn main() {
    let a = Complex::new(3.0, 2.0);
    let b = Complex::new(1.0, 4.0);

    let quotient = a / b;
    
    println!("Quotient of 3+2i / 1+4i = {}", quotient);
}

Executes as:

>>> Quotient of 3+2i / 1+4i = 0.5-0.5i

Power Raised to a i32:

use numerilib::Complex;

fn main() {
    let a = Complex::new(3.0, 2.0);

    let power = a.powi(2);
    
    println!("Power of 3+2i ^ 2 = {}", power);
}

Executes as:

>>> Power of 3+2i ^ 2 = 5+11.999999999999998i

Power Raised to a f64:

use numerilib::Complex;

fn main() {
    let a = Complex::new(3.0, 2.0);

    let power = a.powf(2_f64);
    
    println!("Power of 3+2i ^ 2 = {}", power);
}

Executes as:

>>> Power of 3+2i ^ 2 = 5+11.999999999999998i

Chapter 3: Special Functions

Error Functions

  • Error Function and its relating functions

The Error Function/Erf:

Uses the Maclaurin Series Expansion of the Error Function. The Error Function is the integral of the Gaussian Function. The Error Function is used in statistics to calculate the probability of a given value of a normal distribution.

Parameters:

  • z: The number for which the Error Function will be calculated.

Returns:

The Error Function of the given number.

Equation:

\[ \text{erf}\left( z \right)=\frac{2}{\sqrt{\pi}}\int_{0}^{z}e^{-t^{2}}dt \] \[ \text{erf}\left( z \right)=\frac{2}{\sqrt{\pi}}\sum_{n=0}^{\infty }\frac{\left( -1 \right)^{n}z^{2x+1}}{n!\left( 2n+1 \right)} \]

Example:

use ferrate::special::Error;        
                         
fn main() {                
    let z = 2_f64;        
                          
    let erf = Error::erf(z);
                          
    println!("The Error Function of {} is: {}", z, erf);
}

Executes as:

>>> The Error Function of 2 is: 0.9999999845946841

The Complementary Error Function/Erfc:

Uses the definition of the Complementary Error Function to compute. The Complementary Error Function is the integral of the Gaussian Function.

Parameters:

  • z: The number for which the Complementary Error Function will be calculated.

Returns:

The Complementary Error Function of the given number.

Equation:

\[ \text{erfc}\left( z \right)=1-\text{erf}\left( z \right) \]

Example:

use ferrate::special::Error;                     

fn main() {                                        
    let z = 2_f64;                                 
                                                     
    let erfc = Error::erfc(z);                      
                                                     
    println!("The Complementary Error Function of {} is: {}", z, erfc);
}

Executes as:

>>> The Complementary Error Function of 2 is: 0.000000015405315911820594

The Inverse Error Function

Approximates the Inverse of the Error Function.

Parameters:

  • z: The number for which the Inverse Error Function will be calculated.

Returns:

The Inverse Error Function of the given number.

Equation:

\[ z==\text{erf}^{-1}\left( w \right) \] \[ w==\text{erf}\left( z \right) \]

Example:

use ferrate::special::Error;           

fn main() {                              
    let z = 0.975;                       
                                         
    let inverf = Error::inverf(z);       
                                         
    println!("The Inverse Error Function of {} is: {}", z, inverf);
}

Executes as:

>>> The Inverse Error Function of 0.975 is: 1.5849110680594818

Gamma Functions

  • All the Gamma Functions

Stirling's Approximation:

Uses the Definition of Stirling's Approximation. Stirling's Approximation is an approximation for factorials.

Parameters:

  • n: The number for which the factorial will be approximated via Stirling's Approximation.

Returns:

The approximation of the factorial of the given number.

Equation:

\[ n!\sim \sqrt{2\pi n}\left( \frac{n}{e} \right)^{n} \]

Example:

use ferrate::special::Gamma;                 

fn main() {                               
    let n = 2_f64;                               
    let stirling = Gamma::stirling(n);           
                                             
    println!("The Stirling Approximation of {}! is: {}", n, stirling);
}

Executes as:

>>> The Stirling Approximation of 2! is: 1.9190043514880554

Log Gamma Lanczos:

The Log Gamma Function is calculated via the Lanczos Approximation. The Lanczos Approximation is an algorithm uses precomputed coefficients to perform calculations for the Gamma and Log-Gamma Function with fixed precision.

Parameters:

  • n: The number for which the Log Gamma Function will be calculated.

Returns:

The Log Gamma Function of the given number.

Procedure:

Lanczos Approximation for the Log Gamma Function goes as follows:

1. The Constants are defined:

const G: f64 = 5f64;

const P: [f64; 7] = [
    1.000000000189712,
    76.18009172948503,
    -86.50532032927205,
    24.01409824118972,
    -1.2317395783752254,
    0.0012086577526594748,
    -0.00000539702438713199,
];

2. A loop calculates for the Final Constant:

let A(z) = P[0] + ∑[k=1; P.len()] P[k] / (z + k)];

3. The Log Gamma Function is calculated:

lnΓ(z) = ln(√(2π)) + ln(A(z)) - (z + G - 0.5) * ln(z + G - 0.5) * (z + 0.5);

Example:

use ferrate::special::Gamma;              

fn main() {                                          
    let n = 6_f64;                            
    let lanczos_ln = Gamma::lanczosln(n);     
                                          
    println!("The Log Gamma Function of {} is: {}", n, lanczos_ln);
}

Executes as:

>>> The Log Gamma Function of 6 is: 4.787491742782046

Gamma Function Lanczos:

The Gamma function is a generalization of the factorial function to non-integer and complex numbers. It is defined for all non-positive integers. The Gamma Function is calculated via the Lanczos Approximation. The Lanczos Approximation is an algorithm uses precomputed coefficients to perform calculations for the Gamma and Log-Gamma Function with fixed precision.

Parameters:

  • n: The number for which the Gamma Function will be calculated.

Returns:

The Gamma Function of the given number.

Equation:

\[ \Gamma\left( x \right)=\left( x-1 \right)!\] \[ \Gamma\left( x \right)=\int_{0}^{\infty }t^{x-1}e^{-t}dt\] \[ \Gamma\left( x \right)=e^{ln\Gamma\left( x \right)} \]

Example:

use ferrate::special::Gamma;    

fn main() {                               
    let n = 6_f64;                               
    let lanczos = Gamma::lanczos(n);           
                                             
    println!("The Gamma Function of {} is: {}", n, lanczos);
}

Executes as:

>>> The Gamma Function of 6 is: 120

Incomplete Gamma Function (Lower)

Uses the Series Definition of the Lower Incomplete Gamma Function

Parameters:

  • a: The value for which the Lower Incomplete Gamma Function will be calculated.
  • x: The integral bound at which the Lower Incomplete Gamma Function will be calculated.

Returns:

The Lower Incomplete Gamma Function of the given number.

Equation:

\[ \gamma\left( s,x \right)=\int_{0}^{x }t^{s-1}e^{-t}dt \] \[ \gamma\left( s,x \right)=x^{s}\sum_{k=0}^{\infty }\frac{\left( -1 \right)^{k}x^{k}}{\left( s+k \right)k!} \]

Example:

use ferrate::special::Gamma;               
                
fn main() {                                  
    let a = 3_f64;                            
    let x = 1_f64;                            
                                              
    let gamma = Gamma::incgamma(a, x);        
                                              
    println!("The Lower Incomplete Gamma Function of {} and {} is: {}", a, x, gamma);
}

Executes as:

>>> The Lower Incomplete Gamma Function of 3 and 1 is: 0.16060279414278839

Incomplete Gamma Function (Upper)

Uses the Definition of the Upper Incomplete Gamma Function.

Parameters:

  • a: The value for which the Upper Incomplete Gamma Function will be calculated.
  • x: The integral bound at which the Upper Incomplete Gamma Function will be calculated.

Returns:

The Upper Incomplete Gamma Function of the given number.

Equation:

\[ \Gamma\left( s,x \right)=\int_{x}^{\infty }t^{s-1}e^{-t}dt \] \[ \Gamma\left( s,x \right)+\gamma\left( s,x \right)=\Gamma\left( x \right) \] \[ \Gamma\left( s,x \right)=\Gamma\left( x \right)-\gamma\left( s,x \right) \]

Example:

use ferrate::special::Gamma;              
             
fn main() {                                  
    let a = 3_f64;                            
    let x = 1_f64;                            
                                              
    let gamma = Gamma::incgammac(a, x);       
                                              
    println!("The Upper Incomplete Gamma Function of {} and {} is: {}", a, x, gamma);
}

Executes as:

>>> The Upper Incomplete Gamma Function of 3 and 1 is: 1.8393972058572117

The Regularized Incomplete Gamma Function

Uses the Definition of the Regularized Incomplete Gamma Function.

Parameters:

  • a: The value for which the Regularized Incomplete Gamma Function will be calculated.
  • x: The integral bound at which the Regularized Incomplete Gamma Function will be calculated.

Returns:

The Regularized Incomplete Gamma Function of the given number.

Equation:

\[ \text{P}\left( s,x \right)=\frac{\gamma\left( s,x \right)}{\Gamma\left( x \right)} \]

Example:

use ferrate::special::Gamma;                   
                  
fn main() {                                     
    let a = 3_f64;                               
    let x = 1_f64;                               
                                                 
    let gamma = Gamma::reggamma(a, x);           
                                                 
    println!("The Regularized Incomplete Gamma Function of {} and {} is: {}", a, x, gamma);
}

Executes as:

>>> The Regularized Incomplete Gamma Function of 3 and 1 is: 0.052653017343711174

CDF for Poisson Random Variables

Uses the Definition of the CDF for Poisson Random Variables.

Parameters:

  • a: The value for which the CDF for Poisson Random Variables will be calculated.
  • x: The integral bound at which the CDF for Poisson Random Variables will be calculated.

Returns:

The CDF for Poisson Random Variables of the given number.

Equation:

\[ \text{Q}\left( s,x \right)=\frac{\Gamma\left( s,x \right)}{\Gamma\left( x \right)} \] \[ \text{Q}\left( s,x \right)=1-\text{P}\left( s,x \right) \]

Example:

use ferrate::special::Gamma;                  
                           
fn main() {                                      
    let a = 3_f64;                               
    let x = 1_f64;                               
                                                 
    let gamma = Gamma::reggammac(a, x);          
                                                 
    println!("The CDF for Poisson Random Variables of {} and {} is: {}", a, x, gamma);
}

Executes as:

>>> The CDF for Poisson Random Variables of 3 and 1 is: 0.9473469826562888

Polygamma

  • The polygamma functions, denoted as (\( \psi^{n}\left( z \right) \)), are mathematical functions that are extensions of the digamma function (\( \psi\left( z \right) \)), also known as the psi or di-gamma function. The digamma function is the logarithmic derivative of the gamma function (\( \Gamma\left( z \right) \)), which is a widely used function in mathematics and has applications in various fields.

Digamma:

Calculates the Digamma function with its series representation. The Digamma function is the logarithmic derivative of the gamma function (\( \Gamma\left( z \right) \)).

Parameters:

  • z: The number for which the Digamma function will be calculated.

Returns:

The Digamma function of the given number.

Equation:

\[ \psi\left( z \right)=\frac{d }{dz}\Gamma\left( z \right)=\frac{\Gamma'\left( z \right)}{\Gamma\left( z \right)} \]

Example:

use ferrate::special::Polygamma;         

fn main() {                               
    let z = 2_f64;                        
                                          
    let digamma = Polygamma::digamma(z);  
                                          
    println!("The Digamma Function of {} is: {}", z, digamma);
}

Executes as:

>>> The Digamma Function of 2 is: 0.42278438084235914

Polygamma:

Calculates the Polygamma function with its series representation. The Polygamma function is the (\( m \))-th derivative of the Digamma function (\( \psi\left( z \right) \)).

Parameters:

  • degree: The degree of the Polygamma function. This is the (\( m \)) in (\( \psi^{m}\left( z \right) \)).
  • z: The number for which the Polygamma function will be calculated.

Returns:

The Polygamma function of the given number.

Equation:

\[ \psi^{\left( m \right)}\left( z \right):=\frac{d^{m}}{dz^{m}}\psi\left( z \right)=\left( -1 \right)^{\left( m+1 \right)}m!\sum_{k=0}^{\infty }\frac{1}{\left( z+k \right)^{m+1}}\qquad m\gt 1 \]

\[ \psi^{-1}\left( z \right)=\int\psi\left( z \right)=\text{ln}\Gamma\left( z \right) \]

Example:

use ferrate::special::Polygamma;                
               
fn main() {                                  
    let degree = 0;                                 
    let z = 5_f64;                                  
                                                
    let polygamma = Polygamma::polygamma(degree, z);

    println!("The Polygamma Function of {} at degree {} is: {}", z, degree, polygamma);                                                
}

Executes as:

>>> The Polygamma Function of 5 at degree 0 is: 1.5061177964312848

Beta Functions

  • The beta function is closely related to the gamma function (\( \Gamma\left( z \right) \)) and has various applications in mathematics, statistics, and other fields.

Log Beta:

Uses the definition of Logarithms and the Beta Function to calculate the Logarithm of the Beta Function.

Parameters:

  • z1: The first parameter of the Beta Function.
  • z2: The second parameter of the Beta Function.

Returns:

The Logarithm of the Beta Function of the given numbers.

Equation:

\[ \text{ln}B\left( z_{1},z_{2} \right)=\text{ln}\Gamma\left( z_{1} \right)+\text{ln}\Gamma\left( z_{2} \right)-\text{ln}\Gamma\left( z_{1}+z_{2} \right) \]

Example:

use ferrate::special::Beta;             

fn main() {                            
    let z1 = 1_f64;                    
    let z2 = 2_f64;                    
                                       
    let lnbeta = Beta::lnbeta(z1, z2); 
                                       
    println!("The Log Beta Function of {} and {} is: {}", z1, z2, lnbeta);
}

Executes as:

>>> The Log Beta Function of 1 and 2 is: -0.6931471805616405

Beta Function:

Calculates the Beta Function. The Beta Function is a special function that is closely related to the Gamma Function (\( \Gamma\left( z \right) \)). The Beta Function is used in various fields of mathematics and statistics. For example, the Beta Function is used in the definition of the Student's t-distribution and the F-distribution.

Parameters:

  • z1: The first parameter of the Beta Function.
  • z2: The second parameter of the Beta Function.

Returns:

The Beta Function of the given numbers.

Equation:

\[ B\left( z_{1},z_{2} \right)=\frac{\Gamma\left( z_{1} \right)\Gamma\left( z_{2} \right)}{\Gamma\left( z_{1}+z_{2} \right)} \]

Example:

use ferrate::special::Beta;             

fn main() {                            
    let z1 = 1_f64;                    
    let z2 = 2_f64;                    
                                       
    let beta = Beta::beta(z1, z2);     
                                       
    println!("The Beta Function of {} and {} is: {}", z1, z2, beta);
}

Executes as:

>>> The Beta Function of 1 and 2 is: 0.5

Incomplete Beta Function:

Uses the Definition of an Incomplete Beta Function to calculate. The Incomplete Beta Function is a special function that is closely related to the Beta Function (\( B\left( z_{1},z_{2} \right) \)). The Incomplete Beta Function is used in various fields of mathematics and statistics. For example, the Incomplete Beta Function is used in the definition of the Student's t-distribution and the F-distribution.

Parameters:

  • x: The integral bound at which the Incomplete Beta Function will be calculated.
  • z1: The first parameter of the Incomplete Beta Function.
  • z2: The second parameter of the Incomplete Beta Function.

Returns:

The Incomplete Beta Function of the given numbers.

Equation:

\[ I_{x}\left( z_{1},z_{2} \right)=\frac{B_{x}\left( z_{1},z_{2} \right)}{B\left( z_{1},z_{2} \right)} \]

\[ B_{x}\left( z_{1},z_{2} \right)={B\left( z_{1},z_{2} \right)}\cdot I_{x}\left( z_{1},z_{2} \right) \]

Example:

use ferrate::special::Beta;                

fn main() {                               
    let x = 1_f64 / 7_f64;                 
    let z1 = 0.5_f64;                      
    let z2 = 3_f64;                        
                                           
    let incbeta = Beta::incbeta(x, z1, z2);
                                           
    println!("The Incomplete Beta Function of {}, {} and {} is: {}", x, z1, z2, incbeta);
}

Executes as:

>>> The Incomplete Beta Function of 0.14285714285714285, 0.5 and 3 is: 0.6870211373344728

The Regularized Incomplete Beta Function:

Uses a Series Expansion of the Incomplete Beta Function.

Parameters:

  • x: The integral bound at which the Regularized Incomplete Beta Function will be calculated.
  • z1: The first parameter of the Regularized Incomplete Beta Function.
  • z2: The second parameter of the Regularized Incomplete Beta Function.

Returns:

The Regularized Incomplete Beta Function of the given numbers.

Equation:

\[ I_{x}\left( z_{1},z_{2} \right)=\frac{B_{x}\left( z_{1},z_{2} \right)}{B\left( z_{1},z_{2} \right)} \] \[ I_{x}\left( z_{1},z_{2} \right)=\frac{x^{z_{1}}\left( 1-x \right)^{z_{2}}}{z_{1}B\left( z_{1},z_{2} \right)}\left[ 1+\sum_{n=0}^{\infty }\frac{B\left( z_{1}+1,n+1 \right)}{B\left( z_{1}+z_{2},n+1 \right)}x^{n+1} \right] \]

Example:

use ferrate::special::Beta;                  

fn main() {                                 
    let x = 1_f64 / 7_f64;                   
    let z1 = 1_f64;                          
    let z2 = 2_f64;                          
                                             
    let regincbeta = Beta::regincbeta(z1, z2, x);
                                             
    println!("The Regularized Incomplete Beta Function of {}, {} and {} is: {}", x, z1, z2, regincbeta);
}

Executes as:

>>> The Regularized Incomplete Beta Function of 0.14285714285714285, 1 and 2 is: 0.6440823162530317

The Inverse of the Regularized Incomplete Beta Function:

  • Uses Newton's Method to Calculate the Inverse of the Regularized Incomplete Beta Function.

Parameters:

  • x: The integral bound at which the Regularized Incomplete Beta Function will be calculated.
  • z1: The first parameter of the Regularized Incomplete Beta Function.
  • z2: The second parameter of the Regularized Incomplete Beta Function.

Returns:

The Inverse of the Regularized Incomplete Beta Function of the given numbers.

Equation:

\[ I_{x}^{-1}\left( I_{x}\left( z_{1},z_{2} \right) \right)=1 \]

Example:

use ferrate::special::Beta;                  

fn main() {                                 
    let z1 = 1_f64;                          
    let z2 = 2_f64;                          
    let x = 0.590401_f64;                    
                                             
    let inverse = Beta::invregincbeta(z1, z2, x);
                                             
    println!("The Inverse of the Regularized Incomplete Beta Function of {}, {} and {} is: {}", x, z1, z2, inverse);
}

Executes as:

>>> The Inverse of the Regularized Incomplete Beta Function of 0.590401, 1 and 2 is: 0.3600007812492397