Shammas Polynomials, Pade-Shammas Polynomials, Fourier-Shammas Series

[Albert Chan]

Hi, Namir

Just a suggestion for your reference function, to estimate Ψ(x)

On page 10, your code use: (BTW, your formula had typos, with 2 excess ')')

Ψ(x) = Γ'(x) / Γ(x) ≈ (Γ(x+h) - Γ(x-h))/(2h) / Γ(x)

The problem is Γ(x) have growth rate even faster than e^x
Much more weight is placed to Γ(x+h), and less to Γ(x-h), thus over-estimated Γ'(x).

Instead, we can flatten the curve, and use an equivalent formula (*)

Ψ(x) = (lgamma(x))' ≈ (lgamma(x+h) - lgamma(x-h)) / (2h)

XCas> h := 0.001
XCas> digamma_1(x) := (Gamma(x+h) - Gamma(x-h)) / (2h) / Gamma(x)
XCas> digamma_2(x) := (lgamma(x+h) - lgamma(x-h)) / (2h)
XCas> rel_err(x) := 1 .- [digamma_1(x), digamma_2(x)] ./ Psi(x)

XCas> for(k:=10; k<100; k+=10) {print(k, rel_err(k));}

10, [-8.96829302599e-07, 8.19068590729e-10]
20, [-1.49615930911e-06, 1.45581324773e-10]
30, [-1.92596196058e-06, 5.6931570569e-11]
40, [-2.26519091839e-06, 3.20117266028e-11]
50, [-2.54765957997e-06, 3.02801117513e-11]
60, [-2.79093496314e-06, 2.03200789528e-11]
70, [-3.00534823494e-06, 2.88802315396e-12]
80, [-3.19750830902e-06, 1.85684800869e-12]
90, [-3.37196553724e-06, -4.30255830963e-12]

(*) Googled "matlab lgamma", the equivalent function is named gammaln()