Gamma function: Difference between revisions

Line 2,813:

=={{header|Phix}}==

with javascript_semantics

sequence c = repeat(0,12)

function spouge_gamma(atom z)

function gamma(atom z)

atom accm = c[1]

if accm=0 then

Line 2,834:

Line 2,832:

accm *= exp(-(z+12))*power(z+12,z+0.5) -- Gamma(z+1)

return accm/z

end function

function taylor_gamma(atom x)

-- (good for values between 0 and 1, apparently)

constant t = { 1.00000_00000_00000_00000,

0.57721_56649_01532_86061,

-0.65587_80715_20253_88108,

-0.04200_26350_34095_23553,

0.16653_86113_82291_48950,

-0.04219_77345_55544_33675,

-0.00962_19715_27876_97356,

0.00721_89432_46663_09954,

-0.00116_51675_91859_06511,

-0.00021_52416_74114_95097,

0.00012_80502_82388_11619,

-0.00002_01348_54780_78824,

-0.00000_12504_93482_14267,

0.00000_11330_27231_98170,

-0.00000_02056_33841_69776,

0.00000_00061_16095_10448,

0.00000_00050_02007_64447,

-0.00000_00011_81274_57049,

0.00000_00001_04342_67117,

0.00000_00000_07782_26344,

-0.00000_00000_03696_80562,

0.00000_00000_00510_03703,

-0.00000_00000_00020_58326,

-0.00000_00000_00005_34812,

0.00000_00000_00001_22678,

-0.00000_00000_00000_11813,

0.00000_00000_00000_00119,

0.00000_00000_00000_00141,

-0.00000_00000_00000_00023,

0.00000_00000_00000_00002 }

atom y = x-1,

s = t[$]

for n=length(t)-1 to 1 by -1 do

s = s*y + t[n]

end for

return 1/s

end function

function lanczos_gamma(atom z)

if z<0.5 then

return PI / (sin(PI*z)*lanczos_gamma(1-z))

end if

-- use a lanczos approximation:

atom x = 0.99999999999980993,

t = z + 6.5;

sequence p = { 676.5203681218851,

-1259.1392167224028,

771.32342877765313,

-176.61502916214059,

12.507343278686905,

-0.13857109526572012,

9.9843695780195716e-6,

1.5056327351493116e-7 }

z -= 1

for i=1 to length(p) do

x += p[i] / (z + i)

end for

return sqrt(2*PI) * power(t,z+0.5) * exp(-t) * x

end function

~~procedure~~ sq~~(atom~~ x, atom ~~mul~~)

constant sqPI = sqrt(PI)

atom p = x*mul

~~printf~~(1,~~"%18.16g,%18.15g\n"~~,{x~~<span~~ ~~style="color: #0000FF;">,~~p*p})

procedure sq(sequence zm, string fmt="%19.16f")

atom {z, mul} = zm

atom e = sqPI/mul

sequence s = {spouge_gamma(z),

taylor_gamma(z),

lanczos_gamma(z)},

error = sq_abs(sq_sub(s,e))

string t = join(s,", ",fmt:=fmt)&", "

integer bdx = smallest(error,return_index:=true)

atom best = s[bdx],

p = s[bdx]*mul

for i=1 to length(s) do -- (potentially mark >1)

if s[i]=best then

t[i*22-2..i*22-1] = "*,"

end if

end for

string es = sprintf(fmt,e)

printf(1,"%5g: %s %s, %19.16f\n",{z,t,es,p*p})

end procedure

~~procedure~~ si(~~atom~~ x)

printf(1," z ------ spouge ----- ----- taylor ------ ----- lanczos ----- ---- expected ----- %19.16f\n",PI)

papply({{-3/2,3/4},{-1/2,-1/2},{1/2,1},{1,sqPI},{3/2,2},{2,sqPI},{5/2,4/3},{3,sqPI/2},{7/2,8/15},{4,sqPI/6}},sq)

printf(1,"%18.15f\n",{x})

sq({0.001,sqPI/999.4237725},"%19.15f")

end procedure

sq({0.01,sqPI/99.43258512},"%19.16f")

sq(~~gamma~~(-3/2),~~3/4~~)

sq({0.1,sqPI/9.513507699},"%19.16f")

sq(~~gamma~~(-1/2),-~~1/2~~)

sq({10,sqPI/362880},"%19.12f")

sq(~~gamma~~(1/2),1)

sq({100,sqPI/9.332621544e155},"%19.13g")

si~~<span~~ ~~style="color: #0000FF;">(~~~~gamma~~(1))

if machine_bits()=64 then

sq(~~gamma~~(3/2),2)

sq({150,sqPI/3.808922638e260},"%19.13g") -- (fatal power overflow error on 32 bits)

si~~<span~~ ~~style="color: #0000FF;">(gamma(2~~))

end if

sq(gamma(5/2),4/3)

si(gamma(3))

sq(gamma(7/2),8/15)

si(gamma(4))

The closest to the expected result for each z (row) is marked with a trailing asterisk.

(slightly different results under pwa/p2js...)

The final column is the value of PI (to 16dp) we would get from that best/starred result.

<pre>

z ------ spouge ----- ----- taylor ------ ----- lanczos ----- ---- expected ----- 3.1415926535897932

2.363271801207354, 3.14159265358979

-1.5: 2.3632718012073547*, 2.3632718095606211, 2.3632718012073532, 2.3632718012073547, 3.1415926535897932

-3.544907701811032, 3.14159265358979

-0.5: -3.5449077018110320*, -3.5449077018110306, -3.5449077018110308, -3.5449077018110321, 3.1415926535897932

1.772453850905515, 3.14159265358979

0.5: 1.7724538509055158, 1.7724538509055160*, 1.7724538509055166, 1.7724538509055160, 3.1415926535897932

1.000000000000001

1: 0.9999999999999998, 1.0000000000000000*, 1.0000000000000002, 1.0000000000000000, 3.1415926535897932

0.8862269254527643, 3.14159265358984

1.5: 0.8862269254527577, 0.8862269254527580*, 0.8862269254527583, 0.8862269254527580, 3.1415926535897932

1.000000000000010

2: 0.9999999999999994, 1.0000000000000000*, 1.0000000000000005, 1.0000000000000000, 3.1415926535897932

1.329340388179146, 3.14159265358984

2.5: 1.3293403881791359, 1.3293403881791365*, 1.3293403881791379, 1.3293403881791370, 3.1415926535897906

2.000000000000024

3: 1.9999999999999978, 1.9999999999939679, 2.0000000000000016*, 2.0000000000000000, 3.1415926535897981

3.323350970447942, 3.14159265358998

3.5: 3.3233509704478376, 3.3233509583896768, 3.3233509704478456*, 3.3233509704478426, 3.1415926535897990

6.000000000000175

4: 5.9999999999999884, 5.9999914100724727, 6.0000000000000063*, 6.0000000000000000, 3.1415926535897998

0.001: 999.423772484595421*, 999.423772484595404, 999.423772484595254, 999.423772500000000, 3.1415926534929476

0.01: 99.4325851191505990, 99.4325851191506035*, 99.4325851191505828, 99.4325851200000000, 3.1415926535361195

0.1: 9.5135076986687313, 9.5135076986687318*, 9.5135076986687300, 9.5135076990000000, 3.1415926533710076

10: 362879.999999996094, 0.000000029163, 362880.000000000725*, 362880.000000000000, 3.1415926535898058

100: 9.332621544394e+155, 7.510232292979e-39, 9.332621544394e+155*, 9.332621544e+155, 3.1415926538549222

150: 3.80892263763e+260*, 5.128102530869e-44, 3.80892263763e+260, 3.808922638e+260, 3.1415926529801383

</pre>

=== mpfr version ===

Above translated to mpfr, ~~with~~ ~~higher~~ ~~accuracy~~ ~~and~~ ~~more~~ ~~iterations~~ as ~~per~~ ~~REXX~~, and compared against the builtin.

Above translated to mpfr, but spouge only since there's not much point transferring inherent inaccuracies in taylor/lanczos constants, and compared against the builtin.

without javascript_semantics -- (no mpfr_exp(), mpfr_gamma() in pwa/p2js)

~~requires~~(~~<span~~ ~~style="color: #008000;">"1.0.0"~~) ~~-- (mpfr_set_default_prec[ision] has been renamed)~~

constant dp = 30

string fmt = "%5s: %33s, %33s, %32s\n"

requires("1.0.2") -- (mpfr_get_fixed(maxlen), mpfr_gamma)

include mpfr.e

mpfr_set_default_precision(-87) -- 87 decimal places.

sequence mc = mpfr_inits(40)

~~sequence~~ c = ~~mpfr_inits~~~~<span~~ ~~style="color: #0000FF;">(~~40)

function mpfr_spouge_gamma(mpfr z)

mpfr accm = mc[1]

~~function~~ ~~gamma~~(~~atom~~ z)

if mpfr_cmp_si(accm,0)=0 then

-- mc[1] := sqrt(2*PI)

mpfr accm = c[1]

if mpfr_cmp_si(accm,0)=0 then

-- c[1] := sqrt(2*PI)

mpfr_const_pi(accm)

mpfr_mul_si(accm,accm,2)

Line 2,892:

Line 2,974:

tmk = mpfr_init(),

p = mpfr_init()

for k=2 to length(c) do

for k=2 to length(mc) do

-- c[k] = exp(13-k)*power(13-k,k-1.5)/k1_factrl

-- mc[k] = exp(13-k)*power(13-k,k-1.5)/k1_factrl

mpfr_set_si(tmk,length(c)+1-k)

mpfr_set_si(tmk,length(mc)+1-k)

mpfr_exp(c[k],tmk)

mpfr_exp(mc[k],tmk)

mpfr_set_d(p,k-1.5)

mpfr_pow(p,tmk,p)

mpfr_div(p,p,k1_factrl)

mpfr_mul(c[k],c[k],p)

mpfr_mul(mc[k],mc[k],p)

-- k1_factrl *= -(k-1)

mpfr_mul_si(k1_factrl,k1_factrl,-(k-1))

end for

end if

accm = mpfr_init_set(accm)

accm = mpfr_init_set(accm)

for k=2 to length(c) do

for k=2 to length(mc) do

-- accm += c[k]/(z+k-1)

-- accm += mc[k]/(z+k-1)

mpfr ck = mpfr_init_set(c[k]),

mpfr ck = mpfr_init_set(mc[k]),

zk = ~~mpfr_init~~(z~~+k-1~~)

zk = mpfr_init_set(z)

mpfr_add_si(zk,zk,k-1)

mpfr_div(ck,ck,zk)

mpfr_add(accm,accm,ck)

end for

-- atom zc = z+length(mc)

atom zc = z+length(c)

~~~~-- accm *= exp(-zc)*power(zc,z+0.5) -- Gamma(z+1)

-- accm *= exp(-zc)*power(zc,z+0.5) -- Gamma(z+1)

mpfr ez = ~~mpfr_init~~(-zc),

mpfr p = mpfr_init_set(z),

p = mpfr_init(~~zc~~),

ez = mpfr_init(),

zh = mpfr_init(~~z+~~0.5)

zh = mpfr_init(0.5)

mpfr_add_si(p,p,length(mc))

mpfr_neg(ez,p)

mpfr_exp(ez,ez)

mpfr_add(zh,zh,z)

mpfr_pow(p,p,zh)

mpfr_mul(accm,accm,ez)

mpfr_mul(accm,accm,p)

-- return accm/z

~~mpfr_set_d~~(ez,z)

mpfr_div(accm,accm,z)

mpfr_div(accm,accm,ez)

return accm

end function

~~function~~ ~~gamma2~~~~(atom~~ z)

constant mPI = mpfr_init(),

~~<span~~ ~~style="color:~~ ~~#004080;">mpfr~~ r = mpfr_init(~~z~~)

mqPI = mpfr_init()

~~mpfr_gamma~~(~~r,r~~)

mpfr_const_pi(mPI)

string pistr = mpfr_get_fixed(mPI,dp)

return r

mpfr_sqrt(mqPI,mPI)

function makempfr(object x)

mpfr res

if string(x) then

x = split(x,'/')

res = mpfr_init(x[1])

mpfr_div_si(res,res,to_integer(x[2]))

elsif sequence(x) then

{mpfr x1, object x2} = x

res = mpfr_init_set(x1)

if string(x2) then

mpfr d = mpfr_init(x2)

mpfr_div(res,res,d)

else

mpfr_div_d(res,res,x2)

end if

else

res = mpfr_init(x)

end if

return res

end function

procedure sq(~~mpfr~~ x, integer~~ n,~~ d=1)

procedure mq(sequence zm, integer d=dp)

mpfr p = ~~mpfr_init_set~~(x)

mpfr {z, mul} = apply(zm,makempfr)

~~mpfr_mul_si~~~~<span~~ ~~style="color: #0000FF;">(~~p,p,n)

mpfr s = mpfr_spouge_gamma(z)

~~mpfr_div_si~~(p,p,d)

string t = mpfr_get_fixed(s,d,10,maxlen:=dp+2)

mpfr e = mpfr_init()

mpfr_gamma(e,z)

mpfr p = mpfr_init_set(s)

mpfr_mul(p,p,mul)

mpfr_mul(p,p,p)

string xs = mpfr_get_fixed(x,40),

string zs = mpfr_get_fixed(z,3),

ps = mpfr_get_fixed(p,40)

es = mpfr_get_fixed(e,d,10,maxlen:=dp+2),

~~printf~~(1,~~"%43s,%43s\n"~~,{xs,ps})

ps = mpfr_get_fixed(p,dp,10)

printf(1,fmt,{zs,t,es,ps})

end procedure

printf(1," z %s %s %s\n",{pad(" spouge ",dp+2,"BOTH",'-'),pad(" expected ",dp+2,"BOTH",'-'),pistr})

procedure si(mpfr x)

papply({{-1.5,0.75},{-0.5,-0.5},{0.5,1},{1,{mqPI,1}},{1.5,2},{2,{mqPI,1}},{2.5,"4/3"},{3,{mqPI,2}},{3.5,"8/15"},{4,{mqPI,6}}},mq)

string xs = mpfr_get_fixed(x,40)

~~printf~~(1,~~"%s\n"~~,~~trim_tail(xs~~,~~".0"~~))

mq({"1/1000",{mqPI,"999.4237724845954661149822012996"}},28)

mq({"1/100",{mqPI,"99.43258511915060371353298887051"}},29)

end procedure

mq({"1/10",{mqPI,"9.513507698668731836292487177265"}},30)

sq(~~gamma~~~~(-3/~~2),3,4)

mq({10,{mqPI,362880}},25)

sq(~~gamma~~(-1/~~2),-1~~,2)

mq({100,{mqPI,"9.332621544394415268169923885627e155"}},0)

sq(~~gamma~~(1/2),1)

mq({150,{mqPI,"3.80892263763056972698595524350735e260"}},0)

si(gamma(1))

sq(gamma(3/2),2)

si(gamma(2))

sq(gamma(5/2),4,3)

si(gamma(3))

sq(gamma(7/2),8,15)

si(gamma(4))

puts(1,"mpfr_gamma():\n")

sq(gamma2(-3/2),3,4)

sq(gamma2(-1/2),-1,2)

sq(gamma2(1/2),1)

si(gamma2(1))

sq(gamma2(3/2),2)

si(gamma2(2))

sq(gamma2(5/2),4,3)

si(gamma2(3))

sq(gamma2(7/2),8,15)

si(gamma2(4))

<pre>

z ------------ spouge ------------ ----------- expected ----------- 3.141592653589793238462643383279

2.3632718012073547030642233111215269103967, 3.1415926535897932384626433832795028841972

-1.5: 2.363271801207354703064223311121, 2.363271801207354703064223311121, 3.141592653589793238462643383279

-3.5449077018110320545963349666822903655951, 3.1415926535897932384626433832795028841972

-0.5: -3.544907701811032054596334966e0, -3.544907701811032054596334966e0, 3.141592653589793238462643383279

1.7724538509055160272981674833411451827975, 3.1415926535897932384626433832795028841972

0.5: 1.772453850905516027298167483341, 1.772453850905516027298167483341, 3.141592653589793238462643383279

1

1: 1.000000000000000000000000000000, 1, 3.141592653589793238462643383279

0.8862269254527580136490837416705725913988, 3.1415926535897932384626433832795028841972

1.5: 0.886226925452758013649083741670, 0.886226925452758013649083741670, 3.141592653589793238462643383279

1

2: 1.000000000000000000000000000000, 1, 3.141592653589793238462643383279

1.3293403881791370204736256125058588870982, 3.1415926535897932384626433832795028841972

2.5: 1.329340388179137020473625612505, 1.329340388179137020473625612505, 3.141592653589793238462643383279

2

3: 2.000000000000000000000000000000, 2, 3.141592653589793238462643383279

3.3233509704478425511840640312646472177454, 3.1415926535897932384626433832795028841972

3.5: 3.323350970447842551184064031264, 3.323350970447842551184064031264, 3.141592653589793238462643383279

6

4: 6.000000000000000000000000000000, 6, 3.141592653589793238462643383279

mpfr_gamma():

0.001: 999.4237724845954661149822012996, 999.4237724845954661149822012996, 3.141592653589793238462643383279

2.3632718012073547030642233111215269103967, 3.1415926535897932384626433832795028841972

0.009: 99.43258511915060371353298887051, 99.43258511915060371353298887051, 3.141592653589793238462643383279

-3.5449077018110320545963349666822903655951, 3.1415926535897932384626433832795028841972

0.1: 9.513507698668731836292487177265, 9.513507698668731836292487177265, 3.141592653589793238462643383279

1.7724538509055160272981674833411451827975, 3.1415926535897932384626433832795028841972

10: 362880.0000000000000000000000000, 362880, 3.141592653589793238462643383279

1

100: 9.33262154439441526816992388e155, 9.33262154439441526816992388e155, 3.141592653589793238462643383279

0.8862269254527580136490837416705725913988, 3.1415926535897932384626433832795028841972

150: 3.80892263763056972698595524e260, 3.80892263763056972698595524e260, 3.141592653589793238462643383279

1

1.3293403881791370204736256125058588870982, 3.1415926535897932384626433832795028841972

2

3.3233509704478425511840640312646472177454, 3.1415926535897932384626433832795028841972

6

</pre>