home *** CD-ROM | disk | FTP | other *** search
/ CP/M / CPM_CDROM.iso / simtel / sigm / vols000 / vol082 / exp.pas < prev    next >
Encoding:
Pascal/Delphi Source File  |  1984-04-29  |  1.9 KB  |  97 lines
  1. extern
  2. function exp ( x : real ): real;
  3.  
  4. const
  5. factor2 = 0.50000000000000e+00;
  6. factor3 = 0.16666666666667e+00;
  7. factor4 = 0.41666666666667e-01;
  8. factor5 = 0.83333333333333e-02;
  9. factor6 = 0.13888888888889e-02;
  10. factor7 = 0.19841269841270e-03;
  11. factor8 = 0.24801587301587e-04;
  12. factor9 = 0.27557319223986e-05;
  13. factor10 = 0.27557319223986e-06;
  14. factor11 = 0.25052108385442e-07;
  15. factor12 = 0.20876756987868e-08;
  16. factor13 = 0.16059043836822e-09;
  17. factor14 = 0.11470745597730e-10;
  18. factor15 = 0.76471637318198e-12;
  19. factor16 = 0.47794773323874e-13;
  20. factor17 = 0.28114572543455e-14;
  21. factor18 = 0.15619206968586e-15;
  22. factor19 = 0.82206352466245e-17;
  23. factor20 = 0.41103176233122e-18;
  24. factor21 = 0.19572941063392e-19;
  25.  
  26. type
  27. str = array [1..22] of char;
  28.  
  29. var
  30. result,xa : real;
  31. p : array [1..21] of real;
  32. count,i,j : integer;
  33. sign : char;
  34.  
  35. begin (* exp *)
  36.  
  37. if x >145.0 then exp:=0.99999999999999e+63
  38.   else if x < -145.0 then exp:=0.0
  39. else
  40. begin
  41.  
  42. if x < 0.0 then
  43.   begin
  44.   sign:='-';
  45.   x:=-x;
  46.   end
  47. else sign:='+';
  48.  
  49. (* compute exp(x) for 0<x<1.5 *)
  50. count:=0;
  51. if x > 1.5 then
  52.   if x < 6.0 then begin count:=1; x:=x/4 end
  53.   else if x < 24.0 then begin count:=2; x:=x/16 end
  54.        else if x < 96.0 then
  55.         begin count:=3; x:=x/64 end
  56.             else begin count:=4; x:=x/256.0 end;
  57.  
  58. if x < 0.007 then
  59.   result:=x+1.0
  60. else
  61. begin
  62. xa:=x;
  63. for i:=1 to 21 do
  64.   begin
  65.   p[i]:=xa;
  66.   xa:=xa*x;
  67.   end;
  68.  
  69. result:=
  70.     1+x
  71.     +(p[2]*factor2)
  72.     +(p[3]*factor3)
  73.     +(p[4]*factor4)
  74.     +(p[5]*factor5)
  75.     +(p[6]*factor6)
  76.     +(p[7]*factor7)
  77.     +(p[8]*factor8)
  78.     +(p[9]*factor9)
  79.     +(p[10]*factor10)
  80.     +(p[11]*factor11)
  81.     +(p[12]*factor12)
  82.     +(p[13]*factor13)
  83.     +(p[14]*factor14)
  84.     +(p[15]*factor15)
  85.     +(p[16]*factor16)
  86.     +(p[17]*factor17)
  87.     +(p[18]*factor18);
  88.  
  89. end; (* else *)
  90. for i:=1 to count do result:=sqr( sqr(result) );
  91.  
  92. if sign = '-' then exp:=1.0/result else exp:=result;
  93.  
  94. end; (* else *)
  95.  
  96. end;. (* exp *)
  97.