home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / PASCAL / RMATH.ZIP / MATH.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1990-11-08  |  10.3 KB  |  281 lines
  1. { unit math/math87; last modified: 30OCT90 }
  2.  
  3. { elementary analytical and transcendental math routines }
  4. { Copyright 1990, by J. W. Rider }
  5.  
  6. unit {$ifopt N-} math {$else} math87 {$endif} ;
  7.  
  8. interface
  9.  
  10. {$IFOPT N+} const maxfloatray = 6551;
  11.             type  float=double; xfloat=extended;
  12. {$else}     const maxfloatray = 10919;
  13.             type  float=real;   xfloat=real;
  14. {$endif}
  15.  
  16. type floatray = array [0..maxfloatray] of float;
  17.  
  18. const FLOATEPS = {$ifopt N-} 1.5e-6 {$else} 3.4e-9 {$endif} ;
  19.  
  20. const {strictly overkill on precision}
  21.  
  22. CBRT2     = 1.2599210498948731647672;{ cube root of 2          }
  23. CBRT3     = 1.4422495703074083823216;{ cube root of 3          }
  24. D2R       = 0.017453292519943295769237;{ # radians in a degree }
  25. E         = 2.7182818284590452353603;{ base of natural logs    }
  26. EULERC    = 0.5772156649015328606065;{ Euler's constant        }
  27. LN2       = 0.6931471805599453094172;{ natural logarithm of 2  }
  28. LN3       = 1.0986122886681096913952;{ natural logarithm of 3  }
  29. LN10      = 2.3025850929940456840280;{ natural logarithm of 10 }
  30. LNPI      = 1.1447298858494001741434;{ natural logarithm of PI }
  31. LOG2E     = 1.4426950408889634073599;{ 1/ln(2)                 }
  32. LOG10E    = 0.4342944819032518276511;{ 1/ln(10)                }
  33. M2R       = 0.000290888208665721596154;{  # rads in a minute   }
  34. PI        = 3.1415926535897932384626;{ plain ol' round pi      }
  35. PI_2      = 1.5707963267948966192313;{ pi/2                    }
  36. PI_4      = 0.7853981633974483096156;{ pi/4                    }
  37. R1_E      = 0.3678794411714423215955;{ 1/e                     }
  38. R1_PI     = 0.3183098861837906715378;{ 1/pi                    }
  39. R1_SQRTPI = 0.5641895835477562869481;{ 1/sqrt(pi)              }
  40. R2_PI     = 0.6366197723675813430755;{ 2/pi                    }
  41. R2_SQRTPI = 1.12837916709551257390;  { 2/sqrt(pi)              }
  42. R2D       = 57.295779513082320876798;{ # degrees in a radian   }
  43. S2R       = 0.000004848136811095359936;{  # rads in a second   }
  44. SQRPI     = 9.8696044010893586188345;{ sqr(pi) (not round!)    }
  45. SQRT_2    = 0.707106781186547524401; { sqrt(1/2)               }
  46. SQRT2     = 1.4142135623730950488017;{ sqrt(2)                 }
  47. SQRT2PI   = 2.5066282746310005024158;{ sqrt(2*pi)              }
  48. SQRT3     = 1.7320508075688772935;   { sqrt(3)                 }
  49. SQRT5     = 2.2360679774997896964;   { sqrt(5)                 }
  50. SQRTPI    = 1.7724538509055160272982;{ sqrt(pi) = gamma(1/2)   }
  51. TWOPI     = 6.2831853071795864769253;{ 2*pi                    }
  52.  
  53. function acos   (x: xfloat): xfloat ;
  54. function acosd  (x: xfloat): xfloat ;
  55. function acosh  (x: xfloat): xfloat ;
  56. function asin   (x: xfloat): xfloat ;
  57. function asind  (x: xfloat): xfloat ;
  58. function asinh  (x: xfloat): xfloat ;
  59. function atan   (x: xfloat): xfloat ;
  60. function atand  (x: xfloat): xfloat ;
  61. function atanh  (x: xfloat): xfloat ;
  62. function ceil   (x: xfloat): xfloat ;
  63. function cosd   (x: xfloat): xfloat ;
  64. function cosh   (x: xfloat): xfloat ;
  65. function deg2rad(x: xfloat): xfloat ;
  66. function exp    (x: xfloat): xfloat ;
  67. function exp2   (x: xfloat): xfloat ;
  68. function exp10  (x: xfloat): xfloat ;
  69. function floor  (x: xfloat): xfloat ;
  70. function ln     (x: xfloat): xfloat ;
  71. function log    (x: xfloat): xfloat ;
  72. function log2   (x: xfloat): xfloat ;
  73. function log10  (x: xfloat): xfloat ;
  74. function pow10  (x: xfloat): xfloat ;
  75. function rad2deg(x: xfloat): xfloat ;
  76. function sgn    (x: xfloat): xfloat ;
  77. function sind   (x: xfloat): xfloat ;
  78. function sinh   (x: xfloat): xfloat ;
  79. function sqrt   (x: xfloat): xfloat ;
  80. function tan    (x: xfloat): xfloat ;
  81. function tand   (x: xfloat): xfloat ;
  82. function tanh   (x: xfloat): xfloat ;
  83.  
  84. function atan2  (x,y: xfloat): xfloat ;
  85. function atan2d (x,y: xfloat): xfloat ;
  86. function fdiv   (x,y: xfloat): xfloat ;
  87. function fmod   (x,y: xfloat): xfloat ;
  88. function frem   (x,y: xfloat): xfloat ;
  89. function hypot  (x,y: xfloat): xfloat ;
  90. function ipow   (x,y: xfloat): xfloat ;
  91. function iroot  (x,y: xfloat): xfloat ;
  92. function max    (x,y: xfloat): xfloat ;
  93. function min    (x,y: xfloat): xfloat ;
  94. function pow    (x,y: xfloat): xfloat ;
  95. function root   (x,y: xfloat): xfloat ;
  96.  
  97. function iif    (p: boolean; t,f: xfloat): xfloat ;
  98. function powi   (x: xfloat; n: integer): xfloat ;
  99. function round2 (x: xfloat; n: shortint): xfloat ;
  100.  
  101. function frexp  (x: xfloat; var epart: longint): xfloat ;
  102. function ldexp  (x: xfloat;     epart: longint): xfloat ;
  103.  
  104. function modf   (x: xfloat; var ipart: float): xfloat ;
  105. function remf   (x: xfloat; var ipart: float): xfloat ;
  106.  
  107. function dms2deg (d,m,s: xfloat): xfloat;
  108. function dms2rad (d,m,s: xfloat): xfloat;
  109.  
  110. procedure deg2dms (deg: xfloat; var d,m,s: float);
  111. procedure rad2dms (  r: xfloat; var d,m,s: float);
  112.  
  113. function poly   (x: xfloat; degree: integer; var coeffs ): xfloat;
  114.  
  115.  
  116. implementation
  117.  
  118. function acos; var y:xfloat; begin
  119.    y:=sqrt(1-sqr(x));
  120.    if abs(x)<=y then acos:=PI_2-arctan(x/y)
  121.    else              acos:=arctan(y/x)+iif(x<0,PI,0) end;
  122.  
  123. function acosd; begin acosd:=rad2deg(acos(x)) end;
  124. function acosh; begin acosh:=ln(sqrt(sqr(x)-1)+x) end;
  125.  
  126. function asin; var y:xfloat; begin
  127.    y:=sqrt(1-sqr(x));
  128.    if abs(x)<=y then asin:=arctan(x/y)
  129.    else              asin:=sgn(x)*PI_2-arctan(y/x) end;
  130.  
  131. function asind; begin asind:=rad2deg(asin(x)) end;
  132. function asinh; begin asinh:=system.ln(system.sqrt(sqr(x)+1)+x) end;
  133.  
  134. function atan; begin atan:=arctan(x) end;
  135. function atand; begin atand:=rad2deg(arctan(x)) end;
  136.  
  137. function atan2; begin
  138.       if abs(x)>abs(y) then atan2:=arctan(y/x)+iif(x>0,0,sgn(y)*PI)
  139.  else if abs(x)<abs(y) then atan2:=(PI_2-arctan(x/y))*sgn(y)
  140.  else                       atan2:=sgn(y)*iif(x<0,3,1)*PI_4 end;
  141.  
  142. function atan2d; begin atan2d:=rad2deg(atan2(x,y)) end;
  143. function atanh; begin atanh:=ln((1+x)/(1-x))*0.5 end;
  144.  
  145. function ceil; var intx:xfloat; begin
  146.    intx:=int(x); if intx>=x then ceil:=intx else ceil:=intx+1 end;
  147.  
  148. function cosd; begin cosd:=cos(deg2rad(x)) end;
  149.  
  150. function cosh; var expx:xfloat; begin
  151.    expx:=system.exp(x); cosh:=(sqr(expx)+1)/(expx*2) end;
  152.  
  153. procedure deg2dms; begin s:=remf(remf(deg,d)*60,m)*60 end;
  154. function deg2rad; begin deg2rad:=x*D2R end;
  155.  
  156. function dms2deg; const R1_60 = 1.0/60.0; begin
  157.    dms2deg:=(s*R1_60+m)*R1_60+d end;
  158.  
  159. function dms2rad; begin dms2rad:=d*D2R+m*M2R+s*S2R end;
  160.  
  161. function exp;
  162.    const MAXEXPD = {$IFOPT N-} 88.0296919311
  163.                    {$ELSE } 11356.52340629414394 {$ENDIF};
  164.    begin if x < -MAXEXPD then exp:=0 else exp:=system.exp(x) end;
  165.  
  166. function exp2; begin exp2:=exp(x*LN2) end;
  167. function exp10; begin exp10:=exp(x*LN10) end;
  168. function fdiv; begin fdiv:=int(x/y) end;
  169.  
  170. function floor; var intx:xfloat; begin
  171.    intx:=int(x); if intx<=x then floor:=intx else floor:=intx-1 end;
  172.  
  173. function fmod; begin fmod:=x-floor(x/y)*y end;
  174. function frem; begin frem:=x-int(x/y)*y end;
  175.  
  176. { WARNING: "frexp" accesses the representation of floating point
  177.   numbers at the bit level.  This makes the routine very
  178.   non-portable. }
  179. function frexp;
  180.    type byteray = array [0..pred(sizeof(xfloat))] of byte;
  181.         wordray = array [0..pred(sizeof(xfloat) div 2)] of word;
  182.    begin if x=0 then begin frexp:=0; epart:=0 end
  183.    else begin
  184.       {$IFOPT N-}
  185.          epart:=integer(byteray(x)[0])-128;
  186.          byteray(x)[0]:=128;
  187.       {$ELSE}
  188.          epart:=longint(wordray(x)[4] and $7fff)-16382;
  189.          wordray(x)[4]:=(wordray(x)[4] and $8000) or 16382;
  190.       {$ENDIF}
  191.       frexp:=x end end; { frexp }
  192.  
  193. function hypot; begin hypot:=system.sqrt(sqr(x)+sqr(y)) end;
  194. function iif; begin if p then iif:=t else iif:=f end;
  195.  
  196. function ipow; begin
  197.    if (x>=0) or (y=0) then ipow:= 0
  198.    else ipow:= sin(frac(y*0.5)*TWOPI)*pow(abs(x),y) end;
  199.  
  200. function iroot; begin iroot:=ipow(x,1/y) end;
  201.  
  202. { WARNING: "ldexp" accesses the representation of floating point
  203.   numbers at the bit level.  This makes the routine very
  204.   non-portable. }
  205. function ldexp;
  206.    type byteray = array [0..pred(sizeof(xfloat))] of byte;
  207.         wordray = array [0..pred(sizeof(xfloat) div 2)] of word;
  208.    var y:word;
  209.    begin if x=0 then ldexp:=0
  210.    else begin
  211.       {$IFOPT N-}
  212.          inc(byteray(x)[0],epart);
  213.       {$ELSE}
  214.          y:=wordray(x)[4];
  215.          wordray(x)[4]:=(y and $8000) or word((y and $7fff)+epart);
  216.       {$ENDIF}
  217.       ldexp:=x end end; { ldexp }
  218.  
  219. function ln; begin
  220.    if abs(x-1)<FLOATEPS then ln:=x-1
  221.    else                      ln:=system.ln(abs(x)) end;
  222.  
  223. function log; begin log:=ln(x) end;
  224. function log2; begin log2:=ln(x)*LOG2E end;
  225. function log10; begin log10:=ln(x)*LOG10E end;
  226. function max; begin max:=iif(x>=y,x,y) end;
  227. function min; begin min:=iif(x<=y,x,y) end;
  228. function modf; begin ipart:=floor(x); modf:=x-ipart end;
  229.  
  230. function poly; var y:xfloat; begin
  231.    { coeffs is assumed to be an array [0..degree] of "float" }
  232.    y:=0;
  233.    while degree>=0 do begin
  234.       y:= y*x + floatray(coeffs)[degree]; dec(degree) end;
  235.    poly:=y end;
  236.  
  237. function pow; begin
  238.    if      x>0 then pow:= exp(ln(x)*y)
  239.    else if x=0 then pow:= 0 { pow(0,0)=0 instead of "indeterminate" }
  240.    else if y=0 then pow:= 1
  241.    else             pow:= cos(frac(y*0.5)*TWOPI)*pow(abs(x),y) end;
  242.  
  243. function pow10; begin pow10:=exp10(x) end;
  244.  
  245. function powi; var z:xfloat; begin
  246.    if n>0 then begin
  247.       z:=sqr(powi(x,n shr 1));
  248.       if odd(n) then    powi:=z*x
  249.       else              powi:=z end
  250.    else if n=0 then     powi:=iif(x=0,0,1)
  251.    else { if n<0 then } powi:=1/powi(x,-n) end;
  252.  
  253. function rad2deg; begin rad2deg:=x*R2D end;
  254. procedure rad2dms; begin deg2dms(rad2deg(r),d,m,s) end;
  255. function remf; begin ipart:=int(x); remf:=frac(x) end;
  256. function root; begin root:=pow(x,1/y) end;
  257.  
  258. function round2; var d:xfloat; begin
  259.    d:=powi(10,abs(n));
  260.    if n<0 then round2:=round(x/d)*d
  261.    else        round2:=round(frac(x)*d)/d+int(x) end;
  262.  
  263. function sgn; begin if x=0 then sgn:=0 else sgn:=iif(x>0,1,-1) end;
  264. function sind; begin sind:=sin(deg2rad(x)) end;
  265.  
  266. function sinh; var expx:xfloat; begin
  267.    expx:=system.exp(x); sinh:=(sqr(expx)-1)/(expx*2) end;
  268.  
  269. function sqrt; begin sqrt:=system.sqrt(abs(x))*sgn(x) end;
  270. function tan; begin tan:=sin(x)/cos(x) end;
  271. function tand; begin tand:=tan(deg2rad(x)) end;
  272.  
  273. function tanh;
  274.    const MAXTANHD = {$IFOPT N-} 44.0148459655
  275.                     {$ELSE }  5678.26170314707197 {$ENDIF};
  276.    var exp2x:xfloat;
  277.    begin if abs(x)>=MAXTANHD then tanh:=sgn(x)
  278.    else begin exp2x:=exp(2*x); tanh:=(exp2x-1)/(exp2x+1) end end;
  279.  
  280. end. { Unit MATH/MATH87 }
  281.