123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157 |
- #!/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);
|