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