home *** CD-ROM | disk | FTP | other *** search
/ Current Shareware 1994 January / SHAR194.ISO / engineer / ksprob21.zip / KSMISC.EXE / KSUTILS.PAS < prev    next >
Pascal/Delphi Source File  |  1993-05-19  |  10KB  |  378 lines
  1. unit ksutils;
  2.  
  3. {these are procedures and functions used by the small programs in
  4. the ksmisc archive. }
  5.  
  6. {copyright (C) 1993 Joseph C. Hudson 4903 Algonquin Blvd.
  7.  Clarkston MI 48348}
  8.  
  9. {$B-,D+,E+,F+,L+,N+,R-,V-}
  10.  
  11. interface
  12.  
  13. uses dos;
  14.  
  15. type
  16.   real = double;
  17.  
  18. const
  19.   ksFuzz : real = 1.0E-13; {used in comparing reals}
  20.  
  21.  
  22. var
  23.   rstack : array [0..123] of longint;
  24.   s2pi : real;
  25.   rptr1 : byte;
  26.   ksmissing : real;
  27.  
  28. function iseq(x,y : real) : boolean;
  29. function isne(x,y : real) : boolean; {is x <> y?}
  30. function isge(x,y : real) : boolean; {is x >= y?}
  31. function isgt(x,y : real) : boolean; {is x > y?}
  32. function isle(x,y : real) : boolean; {is x <= y?}
  33. function islt(x,y : real) : boolean; {is x < y?}
  34. function isfile(filn : pathstr) : boolean;
  35. procedure initrand1;
  36. function randnum1 : real;
  37. function pexp(x : real) : real;
  38. function zinv(cdf : real; sigD : byte) : real;
  39. function lngamma(x : real) : real;  {yields about 13 sig digits}
  40. function poispdf(x : integer; mu : real) : real;
  41. function poiscdf(x : integer; mu : real) : real;
  42. function poisinv(cdf, mu : real; var cdflo, cdfhi : real) : integer;
  43.  
  44. implementation
  45.  
  46. var
  47.   a : array[0..7] of byte;
  48.   b : array[0..5] of real;
  49.   i : byte;
  50.  
  51. function iseq(x,y : real) : boolean;
  52. var
  53.   max : real;
  54. begin
  55.   if (x = ksmissing) and (y = ksmissing) then begin iseq := true; exit end;
  56.   if ((x = ksmissing) and (y <> ksmissing)) or
  57.      ((x <> ksmissing) and (y = ksmissing)) then begin iseq := false; exit end;
  58.   if abs(x) < abs(y) then max := abs(y) else max := abs(x);
  59.   if abs(x - y) <= ksFuzz * max then iseq := true else iseq := false
  60. end; {iseq}
  61.  
  62. function isne(x,y : real) : boolean; {is x <> y?}
  63. begin
  64.   isne := not iseq(x,y)
  65. end; {iseq}
  66.  
  67. function isge(x,y : real) : boolean; {is x >= y?}
  68. begin
  69.   if x - y >= - ksFuzz * abs(0.5 * x + 0.5 * y)
  70.     then isge := true else isge := false
  71. end; {isge}
  72.  
  73. function isgt(x,y : real) : boolean; {is x > y?}
  74. begin
  75.   if x - y > ksFuzz * abs(0.5 * x + 0.5 * y)
  76.     then isgt := true else isgt := false
  77. end; {isgt}
  78.  
  79. function isle(x,y : real) : boolean; {is x <= y?}
  80. begin
  81.   if x - y <= ksFuzz * abs(0.5 * x + 0.5 * y)
  82.     then isle := true else isle := false
  83. end; {isle}
  84.  
  85. function islt(x,y : real) : boolean; {is x < y?}
  86. begin
  87.   if x - y < - ksFuzz * abs(0.5 * x + 0.5 * y)
  88.     then islt := true else islt := false
  89. end; {islt}
  90.  
  91. function isfile(filn : pathstr) : boolean;
  92. var
  93.   dirname,curdir : dirstr;
  94.   fname : namestr;
  95.   ext : extstr;
  96.   sname,result : pathstr;
  97. begin
  98.   fsplit(filn,dirname,fname,ext);
  99.   if (length(dirname) > 3) and (dirname[length(dirname)] = '\')
  100.     then dirname := copy(dirname,1,length(dirname) - 1);
  101.   getdir(0,curdir);
  102.   {$I-} chdir(dirname); {$I+}
  103.   if ioresult <> 0 then begin chdir(curdir); isfile := false; exit end;
  104.   sname := fname + ext;
  105.   result := fsearch(sname,dirname);
  106.   if result = '' then isfile := false
  107.                  else isfile := true;
  108.   chdir(curdir)
  109. end; {isfile}
  110.  
  111.  
  112. {the randnum1 generator is from an outline provided by Greg Hasshold. It is
  113.  fast and appears to have a very long cycle and nice randomness properties}
  114.  
  115. procedure initrand1;
  116. var
  117.   i,j : byte;
  118.   dword : array[1..4] of byte;
  119.   rnum : longint absolute dword;
  120. begin
  121.   randomize;
  122.   for i := 0 to 123 do
  123.   begin
  124.     for j := 1 to 4 do dword[j] := random(256);
  125.     rstack[i] := rnum
  126.   end;
  127.   rptr1 := 36
  128. end; {initrand1}
  129.  
  130. function randnum1 : real;
  131. var
  132.   rptr2 : byte;
  133. begin
  134.   inc(rptr1);
  135.   if rptr1 = 124 then rptr1 := 0;
  136.   if rptr1 > 36 then rptr2 := rptr1 - 37 else rptr2 := 87 + rptr1;
  137.   randnum1 := rstack[rptr2] / 4294967295.0 + 0.5;
  138.   rstack[rptr2] := rstack[rptr2] xor rstack[rptr1]
  139. end; {randnum1}
  140.  
  141. function pexp(x : real) : real; {exp protected from underflow and overflow}
  142. begin
  143.   if x < -1419 then pexp := 0.0
  144.     else if x > 709 then pexp := ksmissing
  145.       else pexp := exp(x)
  146. end; {pexp}
  147.  
  148. procedure setSigD(var sigD : byte);
  149. begin
  150.   if sigD = 0 then sigD := 12
  151.     else if sigD < 14 then sigD := sigD + 2
  152.       else sigD := 16
  153. end; {setSigD}
  154.  
  155. function zpdf(z : real) : real;
  156. begin
  157.   if z = ksmissing then begin zpdf := ksmissing; exit end;
  158.   zpdf := pexp(-z * z / 2.0) / s2pi
  159. end; {zpdf}
  160.  
  161. function zcdf(z : real; sigD : byte) : real;
  162. var
  163.   cdf, errTol : real;
  164.   zisNeg : boolean;
  165.  
  166.   procedure useSeries;
  167.   var
  168.     sum,term,zz : real;
  169.     kount : word;
  170.   begin
  171.     zz := z * z; term := z; sum := z; kount := 2;
  172.     while abs(term) > errTol * sum do
  173.     begin
  174.       term := - (term * zz / kount) * ((kount - 1) / (kount + 1));
  175.       sum := sum + term;
  176.       kount := kount + 2
  177.     end;
  178.     cdf := 0.5 + sum / s2pi;
  179.     if zisNeg then cdf := 1.0 - cdf;
  180.   end ; {useSeries}
  181.  
  182.   procedure useConFrac;
  183.   var
  184.     kount : word;
  185.     lnq,term,sum,zz,beta,lnerr,tmpr : real;
  186.   begin
  187.     zz := z * z; term := 1.0; sum := 1.0; kount := 1; beta := 1.0;
  188.     while abs(term) > errtol * sum do
  189.     begin
  190.       beta := 1.0 / (1.0 + kount * beta / zz);
  191.       term := (beta - 1.0) * term;
  192.       sum := sum + term;
  193.       inc(kount)
  194.     end;
  195.     lnq := ln(sum) - zz / 2.0 - ln(s2pi) - ln(z);
  196.     if zisNeg
  197.       then cdf := pexp(lnq)
  198.       else
  199.       begin
  200.         tmpr := pexp(lnq);
  201.         if tmpr = ksmissing then cdf := ksmissing
  202.                             else cdf := 1.0 - tmpr
  203.       end
  204.   end; {useConFrac}
  205.  
  206. begin {zcdf}
  207.   if z = ksmissing then begin zcdf := ksmissing; exit end;
  208.   setSigD(sigD);
  209.   errTol := 0.5 * pexp(-sigD * ln(10.0));
  210.   zisNeg := false; if z < 0.0 then begin zisNeg := true ; z := -z end;
  211.   if z = 0.0 then cdf := 0.5 else if z < 2.25 then useSeries else useConFrac;
  212.   zcdf := cdf
  213. end; {zcdf}
  214.  
  215. function zinvest(cdf : real) : real;
  216. var
  217.   reversed : boolean;
  218.   t, num, den : real;
  219. begin
  220.   if cdf > 0.5 then begin cdf := 1 - cdf; reversed := false end else reversed := true;
  221.   t := sqrt(-2.0 * ln(cdf));
  222.   num := 2.515517 + t * (0.802853 + t * 0.010328);
  223.   den := 1.0 + t * (1.432788 + t * (0.189269 + t * 0.001308));
  224.   t := t - num / den; {pr(Z > t) = cdf}
  225.   if reversed then zinvest := -t else zinvest := t
  226. end; {zinvest}
  227.  
  228. function zinv(cdf : real; sigD : byte) : real;
  229. var
  230.   errTol, z1, z2, p1, p2 : real;
  231.   zisNeg : boolean;
  232. begin
  233.   if isle(cdf,0.0) or isge(cdf,1.0) then begin zinv := ksMissing; exit end;
  234.   setSigD(sigD);
  235.   errTol := 0.5 * pexp(-sigD * ln(10.0));
  236.   zisNeg := false;
  237.   if cdf < 0.5 then begin zisNeg := true; cdf := 1.0 - cdf end;
  238.   if iseq(cdf,0.5) then begin zinv := 0.0; exit end
  239.   else
  240.   begin
  241.     z2 := zinvest(cdf); p2 := zcdf(z2,sigD);
  242.     repeat
  243.       z1 := z2; p1 := p2;
  244.       z2 := z1 - (p1 - cdf) / zpdf(z1);
  245.       p2 := zcdf(z2,sigD)
  246.     until abs(p2 - cdf) < errTol;
  247.     z1 := z2 + (z2 - z1) * (cdf - p2) / (p2 - p1);
  248.     if zisneg then zinv := -z1 else zinv := z1
  249.   end
  250. end; {zinv}
  251.  
  252. function lngamma(x : real) : real;
  253. var
  254.   i : word;
  255.   j,intx : longint;
  256.   y,lngam,y2,tmp1,tmp2 : real;
  257. begin
  258.   if isle(x,0.0) or isgt(x,65535.0) then begin lngamma := ksMissing; exit end;
  259.   if iseq(x,round(x)) then
  260.   begin
  261.     intx := round(x);
  262.     tmp1 := 0.0;
  263.     for i := 2 to intx - 1 do tmp1 := tmp1 + ln(i);
  264.     lngamma := tmp1;
  265.     exit
  266.   end;
  267.   intx := trunc(x);
  268.   if iseq((2.0 * x),2 * intx + 1) then
  269.   begin
  270.     tmp1 := 0.0;
  271.     for j := 2 to intx do tmp1 := tmp1 + ln(2 * j - 1);
  272.     lngamma := tmp1 + 0.5 * ln(pi) - intx * ln(2.0);
  273.     exit
  274.   end;
  275.   if intx < 12 then y := 12.0 + frac(x) else y := x;
  276.   y2 := 1.0 / sqr(y);
  277.   tmp1 := (y - 0.5) * ln(y) - y + 0.5 * ln(2.0 * pi);
  278.   tmp2 := b[5]; for i := 4 downto 0 do tmp2 := b[i] + y2 * tmp2;
  279.   lngam := tmp1 + tmp2 * (1.0 / (12.0 * y));
  280.   if intx < 12 then for i := 0 to 11 - intx do lngam := lngam  - ln(x + i);
  281.   lngamma := lngam
  282. end; {lngamma}
  283.  
  284. function poispdf(x : integer; mu : real) : real;
  285. begin
  286.   if isle(mu,0.0) or (mu = ksmissing) then begin poispdf := ksmissing; exit end;
  287.   poispdf := pexp( x * ln(mu) - mu - lngamma(x+1.0))
  288. end; {poispdf}
  289.  
  290. function poiscdf(x : integer; mu : real) : real;
  291. var
  292.   lnmu, lnpdf, pdf, cdf : real;
  293.   k : word;
  294. begin
  295.   if isle(mu,0.0) or (mu = ksmissing) then begin poiscdf := ksMissing; exit end;
  296.   lnmu := ln(mu);
  297.   lnpdf := -mu;
  298.   pdf := pexp(lnpdf);
  299.   cdf := pdf;
  300.   for k := 1 to x do
  301.   begin
  302.     lnpdf := lnpdf + lnmu - ln(k);
  303.     pdf := pexp(lnpdf);
  304.     cdf := cdf + pdf;
  305.   end;
  306.   poiscdf := cdf
  307. end; {poiscdf}
  308.  
  309. function poisinv(cdf, mu : real; var cdflo, cdfhi : real) : integer;
  310. var
  311.   templ : longint;
  312.   x : integer;
  313.   cdfold, cdfx, pdfx : real;
  314. begin
  315.   if isle(mu,0.0) or (mu = ksmissing) or islt(cdf,0.0) or isge(cdf,1.0) then
  316.   begin
  317.     cdflo := ksMissing;
  318.     cdfhi := ksMissing;
  319.     poisinv := 32767;
  320.     exit
  321.   end;
  322.   if iseq(cdf,0.0) then
  323.   begin
  324.     cdflo := 0.0;
  325.     cdfhi := pexp(-mu);
  326.     poisinv := -1;
  327.     exit
  328.   end;
  329.   templ := round(mu + sqrt(mu) * zinv(cdf,4));
  330.   if templ < 0 then x := 0 else if templ < 32767 then x := templ else
  331.     begin poisinv := 32767; cdflo := ksMissing; cdfhi := ksMissing; exit end;
  332.   cdfx := poiscdf(x,mu);
  333.   pdfx := poispdf(x,mu);
  334.   if iseq(cdf,cdfx) then
  335.   begin
  336.     poisinv := x;
  337.     cdflo := cdfx;
  338.     cdfhi := cdfx + poispdf(x + 1,mu);
  339.   end
  340.   else
  341.   if cdfx < cdf then
  342.   begin
  343.     repeat
  344.       inc(x);
  345.       cdfold := cdfx;
  346.       pdfx := pdfx * mu / x;
  347.       cdfx := cdfx + pdfx;
  348.     until cdfx > cdf;
  349.     poisinv := x - 1;
  350.     cdflo := cdfold;
  351.     cdfhi := cdfx
  352.   end
  353.   else
  354.   begin
  355.     repeat
  356.       dec(x);
  357.       cdfold := cdfx;
  358.       cdfx := cdfx - pdfx;
  359.       pdfx := pdfx * (x + 1.0) / mu
  360.     until cdfx <= cdf;
  361.     poisinv := x;
  362.     cdflo := cdfx;
  363.     cdfhi := cdfold
  364.   end
  365. end; {poisinv}
  366.  
  367. begin
  368.   s2pi := sqrt(2.0 * pi);
  369.   b[0] :=  1.0;
  370.   b[1] := -1.0   /    30.0;
  371.   b[2] :=  1.0   /   105.0;
  372.   b[3] := -1.0   /   140.0;
  373.   b[4] :=  1.0   /    99.0;
  374.   b[5] := -691.0 / 30300.0;
  375.   for i := 0 to 5 do a[i] := 0; a[6] := $f0; a[7] := $7f;
  376.   ksmissing := real(a)
  377. end.
  378.