home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 June / SIMTEL_0692.cdr / msdos / turbopas / pas_sci.arc / ERFD3.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1985-09-28  |  1.3 KB  |  81 lines
  1. program erfd3;        { -> 330 }
  2. { evaluation of the gaussian error function }
  3.  
  4. var    x,er,ec        : real;
  5.     done        : boolean;
  6.  
  7. function erf(x: real): real;
  8. { infinite series expansion of the Gaussian error function }
  9.  
  10. const    sqrtpi        = 1.7724538;
  11.     tol        = 1.0E-4;
  12.  
  13. var    x2,sum,sum1,term: real;
  14.     i        : integer;
  15.  
  16. begin
  17.       x2:=x*x;
  18.       sum:=x;
  19.       term:=x;
  20.       i:=0;
  21.       repeat
  22.     i:=i+1;
  23.     sum1:=sum;
  24.     term:=2.0*term*x2/(1.0+2.0*i);
  25.     sum:=term+sum1
  26.       until term<tol*sum;
  27.       erf:=2.0*sum*exp(-x2)/sqrtpi
  28. end;    { erf }
  29.  
  30. function erfc(x: real): real;
  31. { complement of error function }
  32. const    sqrtpi        = 1.7724538;
  33.     terms        = 12;
  34.  
  35. var    x2,u,v,sum    : real;
  36.     i        : integer;
  37. begin
  38.   x2:=x*x;
  39.   v:=1.0/(2.0*x2);
  40.   u:=1.0+v*(terms+1.0);
  41.   for i:=terms downto 1 do
  42.     begin
  43.       sum:=1.0+i*v/u;
  44.       u:=sum
  45.     end;
  46.   erfc:=exp(-x2)/(x*sum*sqrtpi)
  47. end;        { ercf }
  48.  
  49. begin        { main }
  50.   ClrScr;
  51.   done:=false;
  52.   writeln;
  53.   repeat
  54.     write('Arg? ');
  55.     readln(x);
  56.     if x<0.0 then done:=true
  57.     else
  58.       begin
  59.     if x=0.0 then
  60.       begin
  61.         er:=0.0;
  62.         ec:=1.0
  63.       end
  64.     else
  65.       begin
  66.         if x<1.5 then
  67.           begin
  68.         er:=erf(x);
  69.         ec:=1.0-er
  70.           end
  71.         else
  72.           begin
  73.         ec:=erfc(x);
  74.         er:=1.0-ec
  75.         end    { if }
  76.     end;
  77.     writeln('X= ',x:8:4,' Erf= ',er:12:8,', Erfc= ',ec:12)
  78.       end    { if }
  79.     until done
  80. end.
  81.