Anonymous user
Gamma function: Difference between revisions
→{{header|REXX}}: added a 2nd version (using Sponge's approximation).
No edit summary |
(→{{header|REXX}}: added a 2nd version (using Sponge's approximation).) |
||
Line 2,955:
=={{header|REXX}}==
This version uses a Taylor series with 80-digits coefficients with much more accuracy.
<br>As a result, the gamma value for <big> ½ </big> is now 25 decimal digits more accurate than the previous version
<br>(which only used 20 digit coefficients).
Line 3,043 ⟶ 3,044:
</span>
<br><br>to 110 digits past the decimal point, the vinculum (overbar) marks the ''difference digit'' from the computed value (by this REXX program) of gamma(<b><big>½</big></b>). <br><br>
===Spouge's approximation, using 87 digit coefficients===
{{trans|Phix}}
{{trans|C}}
This REXX version is a translation of '''Phix''' but with more (decimal digits) precision and more ''steps''.
Many of the "normal" high-level mathematical functions aren't available in REXX, so some of them (RYO) are included here.
<lang rexx>/*REXX program calculates the gamma function using Spouge's approximation with 87 digits*/
e=2.71828182845904523536028747135266249775724709369995957496696762772407663035354759457138
numeric digits length(e) - length(.) /*use the number of decimal digits in E*/
c.= 0
# = 40 /*#: the number of steps in GAMMA func*/
call sq gamma(-3/2), 3/4
call sq gamma(-1/2), -1/2
call sq gamma( 1/2), 1
call si gamma( 1 )
call sq gamma( 3/2), 2
call si gamma( 2 )
call sq gamma( 5/2), 4/3
call si gamma( 3 )
call sq gamma( 7/2), 8/15
call si gamma( 4 )
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
gamma: procedure expose c. e #; parse arg z; #p= # + 1
accm = c.1
if accm==0 then do; accm= sqrt( 2*pi() )
c.1 = accm
kfact = 1
do k=2 to #
c.k= exp(#p-k) * pow(#p-k, k-1.5) / kfact
kfact = kfact * -(k-1)
end /*k*/
end
do j=2 to #; accm = accm + c.j / (z+j-1)
end /*k*/
return (accm * exp(-(z+#)) * pow(z+#, z+0.5) ) / z
/*──────────────────────────────────────────────────────────────────────────────────────*/
pi: return 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348
fmt: parse arg n,p,a; _= format(n,p,a); L= length(_); return left( strip0(_), L)
isInt: return datatype(arg(1), 'W') /*is the argument an integer? */
sq: procedure expose #; parse arg x,mu; say fmt(x,9,#) fmt((x*mu)**2,9,#); return
si: procedure expose #; parse arg x; say fmt(x,9,#); return
strip0: procedure; arg _; if pos(., _)\==0 then _= strip(strip(_,'T',0),'T',.); return _
/*──────────────────────────────────────────────────────────────────────────────────────*/
exp: procedure expose e; arg x; ix= x%1; if abs(x-ix)>.5 then ix=ix+sign(x); x= x-ix; z=1
_=1; w=1; do j=1; _= _*x/j; z= (z+_)/1; if z==w then leave; w=z
end /*j*/; if z\==0 then z= e**ix * z; return z
/*──────────────────────────────────────────────────────────────────────────────────────*/
ln: procedure; parse arg x; call e; ig= x>1.5; is= 1-2*(ig\==1); ii= 0; xx= x
do while ig & xx>1.5 | \ig & xx<.5; _=e
do k=-1; iz=xx*_**-is; if k>=0&(ig&iz<1|\ig&iz>.5) then leave; _=_*_; izz=iz; end
xx= izz; ii= ii+is*2**k; end /*while*/; x= x*e**-ii-1; z=0; _= -1; p=z
do k=1; _=-_*x; z=z+_/k; if z=p then leave; p=z; end; /*k*/; return z+ii
/*──────────────────────────────────────────────────────────────────────────────────────*/
pow: procedure; parse arg x,y; if y=0 then return 1; if x=0 then return 0
if isInt(y) then return x**y; if isInt(1/y) then return root(x, 1/y)
if abs(y//1)=.5 then return sqrt(x)**sign(y)*x**(y%1); return exp( y*ln(x) )
/*──────────────────────────────────────────────────────────────────────────────────────*/
root: procedure; parse arg x 1 ox,y 1 oy; if x=0 | y=1 then return x/1
if \isInt(y) then return $pow(x, 1/y)
if y==2 then return sqrt(x); if y==-2 then return 1/sqrt(x); return rooti(x,y)/1
/*──────────────────────────────────────────────────────────────────────────────────────*/
rooti: x=abs(x); y=abs(y); a= digits() + 5; m= y-1; d= 5
parse value format(x,2,1,,0) 'E0' with ? 'E' _ .; g= (?/y'E'_ % y) + (x>1)
do until d==a; d=min(d+d, a); numeric digits d; o=0
do until o=g; o=g; g= format((m*g**y+x)/y/g**m,,d-2); end; end
_= g*sign(ox); if oy<0 then _= 1/_; return _
/*──────────────────────────────────────────────────────────────────────────────────────*/
sqrt: procedure; parse arg x; if x=0 then return 0; d=digits(); numeric digits; h=d+6
numeric form; m.=9; parse value format(x,2,1,,0) 'E0' with g "E" _ .; g=g *.5'e'_ %2
do j=0 while h>9; m.j=h; h=h%2+1; end /*j*/
do k=j+5 to 0 by -1; numeric digits m.k; g=(g+x/g)*.5; end /*k*/
numeric digits d; return g/1</lang>
{{out|output|text= when using the default input:}}
<pre>
2.3632718012073547030642233111215269103967 3.1415926535897932384626433832795028841972
-3.5449077018110320545963349666822903655951 3.1415926535897932384626433832795028841972
1.7724538509055160272981674833411451827975 3.1415926535897932384626433832795028841972
1
0.8862269254527580136490837416705725913988 3.1415926535897932384626433832795028841972
1
1.3293403881791370204736256125058588870982 3.1415926535897932384626433832795028841972
2
3.3233509704478425511840640312646472177454 3.1415926535897932384626433832795028841972
6
</pre>
=={{header|Ring}}==
| |||