home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / PASCAL / RMATH.ZIP / SPECFUN.PAS < prev   
Encoding:
Pascal/Delphi Source File  |  1990-11-01  |  6.1 KB  |  191 lines
  1. { unit SpecFun/SpecF87 }
  2. { Special functions for mathematical calculations }
  3. { Copyright 1990, by J. W. Rider. }
  4.  
  5. unit {$IFOPT N-} SpecFun {$ELSE} SpecF87 {$ENDIF} ;
  6.  
  7. interface
  8.  
  9. uses {$IFOPT N-} Math {$ELSE} Math87 {$ENDIF} ;
  10.  
  11. function factorial  (n:word):xfloat;
  12. function lnfactorial(n:word):xfloat;
  13.  
  14. function permutation(n,k:word):xfloat;
  15. function combination(n,k:word):xfloat;
  16.  
  17. function gamma  (x:xfloat):xfloat;
  18. function lngamma(x:xfloat):xfloat;
  19.  
  20. function beta  (x,y:xfloat):xfloat;
  21. function lnbeta(x,y:xfloat):xfloat;
  22.  
  23. function igamma (a,x:xfloat):xfloat;
  24. function igammac(a,x:xfloat):xfloat;
  25. function igammap(a,x:xfloat):xfloat;
  26. function igammaq(a,x:xfloat):xfloat;
  27.  
  28. function ibeta (a,b,x:xfloat):xfloat;
  29. function ibetap(a,b,x:xfloat):xfloat;
  30. function ibetaq(a,b,x:xfloat):xfloat;
  31.  
  32. function erf (x:xfloat):xfloat;
  33. function erfc(x:xfloat):xfloat;
  34.  
  35.  
  36. implementation
  37.  
  38. { IBETACF is a "continued fraction" evaluation of an incomplete
  39.   beta auxiliary function. (See HMF: eqn. 26.5.8; NR: BETACF ) }
  40.  
  41. function ibetacf(a,b,x:xfloat):xfloat;
  42. var am,ap,apb,app,azlast,az,bm,bp,bpp,bz,d,deca,inca:xfloat;
  43.     m,tm:longint;
  44. begin
  45.    am:=1; az:=1; bm:=1; m:=1; apb:=a+b; deca:=a-1; inca:=a+1;
  46.    bz:=1-apb*x/inca; tm:=2; d:=(b-1)*x/(inca*(a+2));
  47.    ap:=d+1; bp:=bz+d; d:=-inca*(apb+1)*x/((a+2)*(a+3));
  48.    app:=ap+d; bpp:=bp+d*bz; am:=ap/bpp; bm:=bp/bpp;
  49.    azlast:=az; az:=app/bpp;
  50.    while abs(az-azlast)>(abs(az)*FLOATEPS) do begin
  51.       inc(m); inc(tm,2); d:=m*(b-m)*x/((deca+tm)*(a+tm));
  52.       ap:=az+d*am; bp:=d*bm+1; d:=-(a+m)*(apb+m)*x/((a+tm)*(inca+tm));
  53.       app:=ap+d*az; bpp:=bp+d end;
  54.    ibetacf:=az end; { ibetacf }
  55.  
  56.  
  57. { IGAMMACF is a "continued fraction" evaluation of an incomplete
  58.   gamma auxiliary function.  ( See HMF: eqn. 6.5.31; NR: GCF ) }
  59.  
  60. function igammacf(a,x:xfloat):xfloat;
  61. var n:longint; a0,a1,b0,b1,f,fn,g,glast,na:xfloat;
  62. begin
  63.    a0:=1; a1:=x; b0:=0; b1:=1; f:=1; g:=0; glast:=1; n:=0;
  64.    repeat
  65.       inc(n); fn:=f*n; na:=n-a;
  66.       a0:=(a1+a0*na)*f; b0:=(b1+b0*na)*f;
  67.       a1:=x*a0+fn*a1;   b1:=x*b0+fn*b1;
  68.       if (a1<>0) then begin f:=1/a1; glast:=g; g:=b1*f end;
  69.       until abs(g-glast)<(abs(g)*FLOATEPS);
  70.    igammacf:=g; end; { igammacf }
  71.  
  72.  
  73. { IGAMMASER is a "power series" evaluation of the incomplete
  74.   gamma auxiliary function.  ( See HMF: eqn. 6.5.29; NR: GSER ) }
  75.  
  76. function igammaser(a,x:xfloat):xfloat; var d,s:xfloat;
  77. begin if x=0 then igammaser:=0
  78.       else begin
  79.          d:=1/a; s:=d;
  80.          repeat
  81.             a:=a+1; d:=d*x/a; s:=s+d;
  82.             until abs(d)<(abs(s)*FLOATEPS);
  83.          igammaser:=s end end; { igammaser }
  84.  
  85.  
  86. { *** interface'd functions in alphabetical order *** }
  87.  
  88. function beta; begin beta:=exp(lnbeta(x,y)) end;
  89.  
  90. function combination; begin
  91.    combination:=round(exp(lnfactorial(n)
  92.                 -lnfactorial(k)-lnfactorial(n-k))) end;
  93.  
  94. function erf; var y,z:xfloat;
  95. begin if abs(x)<FLOATEPS then erf:=x*R2_SQRTPI
  96.       else if abs(x)>10 then erf:=sgn(x)
  97.       else begin
  98.          y:=sqr(x); z:=x*exp(-y-SQRTPI);
  99.          if y<1.5 then erf:=z*igammaser(0.5,y)
  100.          else erf:=sgn(x)-z*igammacf(0.5,y) end end; { erf }
  101.  
  102. function erfc; begin erfc:=1-erf(x) end;
  103.  
  104. function factorial;
  105. const MAXATYPE = 32;
  106. type atype = array[0..MAXATYPE] of xfloat;
  107. const a : ^atype = NIL; atop:shortint = -1;
  108. var j:byte;
  109. begin
  110.    if a=NIL then begin
  111.       new(a); if a<>NIL then begin a^[0]:=1; atop:=0 end end;
  112.    if (a<>NIL) and (n>atop) and (n<=MAXATYPE) then begin
  113.       for j:=succ(atop) to n do a^[j]:=a^[pred(j)]*j;
  114.       atop:=n end;
  115.    if (n<=atop) then factorial:=a^[n]
  116.    else              factorial:=round(gamma(succ(n))) end; { factorial }
  117.  
  118. function gamma; var y:xfloat; begin
  119.    y:=exp(lngamma(x));
  120.    if x<0 then if not odd(trunc(x)) then y:=-y;
  121.    gamma:=y end;
  122.  
  123. function ibeta; var y:xfloat; begin
  124.    if      x<=0 then ibeta:=0
  125.    else if x>=1 then ibeta:=beta(a,b)
  126.    else begin
  127.       y:=exp(a*ln(x)+b*ln(1-x));
  128.       if (x*(a+b+2))<=(a+1) then ibeta:=y*ibetacf(a,b,x)/a
  129.       else ibeta:=beta(a,b)-y*ibetacf(b,a,1-x)/b end end; { ibeta }
  130.  
  131. function ibetap; var y:xfloat; begin
  132.    if      x<=0 then ibetap:=0
  133.    else if x>=1 then ibetap:=1
  134.    else begin
  135.       y:=exp(-lnbeta(a,b)+a*ln(x)+b*ln(1-x));
  136.       if x<=(a+1)/(a+b+2) then ibetap:=y*ibetacf(a,b,x)/a
  137.       else ibetap:=1-y*ibetacf(b,a,1-x)/b end end; { ibetap }
  138.  
  139. function ibetaq; begin ibetaq:=1-ibetap(b,a,1-x) end;
  140.  
  141. function igamma; begin
  142.    if      x<=0    then igamma:=0
  143.    else if x<(a+1) then igamma:=exp(-x+a*ln(x))*igammaser(a,x)
  144.    else                 igamma:=gamma(a)-igammac(a,x) end;
  145.  
  146. function igammac; begin
  147.    if x<(a+1) then igammac:=gamma(a)-igamma(a,x)
  148.    else            igammac:=exp(-x+a*ln(x))*igammacf(a,x) end;
  149.  
  150. function igammap; begin
  151.    if      x<=0    then igammap:=0
  152.    else if x<(a+1) then igammap:=exp(-x+a*ln(x)-lngamma(a))
  153.                                  *igammaser(a,x)
  154.    else                 igammap:=1-igammaq(a,x) end;
  155.  
  156. function igammaq; begin
  157.    if x<(a+1) then igammaq:=1-igammap(a,x)
  158.    else            igammaq:=exp(-x+a*ln(x)-lngamma(a))
  159.                             *igammacf(a,x) end;
  160.  
  161. function lnbeta; begin lnbeta:=lngamma(x)+lngamma(y)-lngamma(x+y) end;
  162.  
  163. function lnfactorial;
  164. const MAXATYPE = 99;
  165. type atype = array[0..MAXATYPE] of xfloat;
  166. const a:^atype = NIL;
  167. var j:byte;
  168. begin
  169.    if a=NIL then begin
  170.       new(a);
  171.       if a<>NIL then for j:=0 to MAXATYPE do a^[j]:=-1 end;
  172.    if (n<MAXATYPE) and (a<>NIL) then begin
  173.       if a^[n]<0 then a^[n]:=lngamma(succ(n));
  174.       lnfactorial:=a^[n] end
  175.    else lnfactorial:=lngamma(succ(n)) end; { lnfactorial }
  176.  
  177. function lngamma; var s,t:xfloat; begin
  178.    if      x<0 then lngamma:=LNPI-lngamma(1-x)-ln(sin(frac(x*0.5)*TWOPI))
  179.    else if x<1 then lngamma:=lngamma(x+1)-ln(x)
  180.    else begin
  181.       t:=x+4.5; t:=(x-0.5)*ln(t)-t;
  182.       s:=1 + 76.18009173/x     - 86.50532033  /(x+1) + 24.01409822/(x+2)
  183.            - 1.231739516/(x+3) + 1.20858003e-3/(x+4) - 5.3682E-6  /(x+5);
  184.       lngamma:=t+ln(SQRT2PI*s) end end; { lngamma }
  185.  
  186. function permutation; begin
  187.    permutation:=round(exp(lnfactorial(n)-lnfactorial(n-k))) end;
  188.  
  189. end. { unit specfun/specf87 }
  190.  
  191.