sidef-scripts/Math/factorial_approximations.sf

#!/usr/bin/ruby

define ℯ = Num.e
define π = Num.pi
define τ = Num.tau
define γ = Num.EulerGamma
define ϕ = Num.phi

# Unknown
# sqrt(2*pi*e) * e^(-((3 - sqrt(3)) / 6)) * (((3 - sqrt(3)) / 6 + n) / e)^(n + 1/2)
func f1(n) {
    define w = ((3 - √3) / 6)
    define d = (√(τ * ℯ) * exp(-w))
    d * pow((n+w) / ℯ, n + 1/2)
}

# S. Ramanujan's formula (simplified)
# sqrt(pi) * (n/e)^n * (n*(4*n*(2*n + 1) + 1))^(1/6)
func f2(n) {
    define d = √π
    d * exp(-n) * n**n * root(n*(4*n*(2*n + 1) + 1), 6)
}

# Stirling's formula + Trizen's approximation
# sqrt(2*pi*n) * (n/e)^n * (digamma(n+1) + gamma) / (log(n + 1/(phi * sqrt(n))) + gamma) / zeta(n)
func f3(n) {
    (√(τ*n) * exp(-n) * n**n * ((digamma(n+1) + γ) / (log(n + 1/(ϕ * √n)) + γ)) / zeta(n))
}

# Stirling's formula
# sqrt(2*pi*n) * (n/e)^n
func f4(n) {
    √(τ * n) * exp(-n) * n**n
}

# Stirling's formula + Laplace's method (1)
# sqrt(2*pi*n) * (n/e)^n * (1 + 1/(12*n))
func f5(n) {
    √(τ * n) * exp(-n) * n**n * (1 + 1/(12*n))
}

# S. Ramanujan (2)
# sqrt(pi) * n^n * e^(-n) * ((8 * n^3) + (4 * n^2) + n + 1/30)^(1/6)
func f6(n) {
    define w = √π
    w * pow(n, n) * exp(-n) * pow(((8 * n**3) + (4 * n**2) + n + 1/30), 1/6)
}

# N. Batir
# sqrt(2*pi) * n^n * e^(-n) * sqrt(n + 1/2) * e^(-(1 / (6 * (n + 3/8))))
func f7(n) {
    define w = √τ
    w * pow(n, n) * exp(-n) * √(n + 1/2) * exp(-(1 / (6 * (n + 3/8))))
}

# N. Batir(2)
# sqrt(2*pi) * n^n * e^(-n) * sqrt(n + 1/6 + 1/(72*n) - 31/(6480 * n^2) - 139/(155520 * n^3) + 9871/(6531840 * n^4))
func f8(n) {
    define w = √τ
    w * pow(n, n) * exp(-n) * √(n + 1/6 + 1/(72*n) - 31/(6480 * n**2) - 139/(155520 * n**3) + 9871/(6531840 * n**4))
}

# N. Batir(2) derivation by Wolfram|Alpha
# 1/216 * sqrt(pi/70) * exp(-n) * n^(n-2) * sqrt(42*n*(24*n*(90*n*(12*n*(6*n + 1) + 1) - 31) - 139) + 9871)
func f9(n) {
    1/216 * √(π/70) * exp(-n) * n**(n-2) * √(42*n*(24*n*(90*n*(12*n*(6*n + 1) + 1) - 31) - 139) + 9871);
}

# Gamma approximation
# n! ~ (1-i) * (-1)^(1/4) * sqrt(pi) * 12^(-n) * exp(-2*n) * sqrt(n) * ((12*e*n^2 + e)/n)^n
# n! ~ sqrt(pi) * 2^(1/2 - 2*n) * exp(-n) * n^(1/2-n) * (4*n^2 + 1/3)^n
func f10(n) {
    sqrt(π) * 2**(1/2 - 2*n) * exp(-n) * n**(1/2 - n) * (4 * n**2 + 1/3)**n
}

# Gamma approximation (2)
# n! ~ limit_{k -> Infinity} sqrt(2*pi) * ((12*n^2 + 24*k*n + 12*k^2 + 1)/(12*e*n + 12*k*e))^(n+k) / ((n+k-1)!/(n!)) / sqrt(n+k)
func f11(n) {
    define k = 50       # higher value of k => higher precision
    sqrt(τ) * ((12*n**2 + (24*k)*n + 12*k**2 + 1)/(12*ℯ*n + 12*k*ℯ))**(n+k) / prod(1..^k, {|j| n + j }) / sqrt(n+k)
}

# Gamma approximation (3)
# n! ~ (sqrt(2*pi) * ((12*n^2 + 24*n + 12 + 1)/(12*e*n + 12*e))^(n + 1))/sqrt(n + 1)
func f12(n) {
    sqrt(τ) * ((12*n**2 + 24*n + 12 + 1)/(12*ℯ*n + 12*ℯ))**(n+1) / sqrt(n+1)
}

const FORMULAS = [
    [Unknown            => f1],
    [Ramanujan          => f2],
    ["Stirling+Trizen"  => f3],
    [Stirling           => f4],
    ["Stirling+Laplace" => f5],
    [Ramanujan2         => f6],
    [Batir              => f7],
    [Batir2             => f8],
    [Batir2WA           => f9],
    [Gamma              => f10],
    [Gamma2             => f11],
    [Gamma3             => f12],
]

#
## TESTS
#

var report = Hash()
for n in (15..30) {

    for name,f in (FORMULAS) {
        report{name} := [] << f(n)
    }

    report{:EXACT} := [] << n!

    say "\n=> #{n}!"
    for k,v in (report.sort_by {|_,v| v }) {
        "%20s: %s\n".printf(k, v[-1]);
    }
}

func A(a) { a.sum / a.len }
func G(a) { a.prod.root(a.len) }
func H(a) { a.len / a.map{1/_}.sum }

var r1 = Hash()
var r2 = Hash()
var r3 = Hash()

for k,v in (report) {
    k == :EXACT && next;
    var arr = v.range.map {|i| abs(report{:EXACT}[i] - v[i])}
    r1{k} = A(arr);
    r2{k} = G(arr);
    r3{k} = H(arr);
}

func print_report(h) {
    for k,v in (h.sort_by {|_,v| v }) {
        "%20s: %s\n".printf(k, v);
    }
}

say ("\n", '-'*80);

say "\n=>> A <<=";
print_report(r1);

say "\n=>> G <<=";
print_report(r2)

say "\n=>> H <<=";
print_report(r3);

say ("\n", '-'*80);