home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / PASCAL / RMATH.ZIP / CMATH.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1990-08-17  |  4.9 KB  |  142 lines
  1. { unit CMATH/CMATH87; last modfied 17AUG90 }
  2.  
  3. { Complex mathematics for Turbo Pascal }
  4. { Copyright 1990, by J. W. Rider }
  5.  
  6. unit {$ifopt N-} cmath {$else} cmath87 {$endif} ;
  7.  
  8. interface
  9.  
  10. uses {$ifopt N-} math {$else} math87 {$endif} ;
  11.  
  12. const maxcxray = {$ifopt N-} 5459 {$else} 3275 {$endif} ;
  13.  
  14. type complex = record rp, ip: float; end;
  15.      cxray = array [0..maxcxray] of complex;
  16.  
  17. function cabs (x: complex): float;
  18. function phase(x: complex): float;
  19.  
  20. procedure conj  (var x: complex); { complex conjugate }
  21. procedure cneg  (var x: complex); { x:=-x      }
  22. procedure cjx   (var x: complex); { x:= j*x, j = sqrt(-1) }
  23. procedure cnjx  (var x: complex); { x:=-j*x, j = sqrt(-1) }
  24. procedure crecip(var x: complex); { x:=1/x     }
  25. procedure csqr  (var x: complex); { x:=x*x     }
  26. procedure csqrt (var x: complex); { x:=sqrt(x) }
  27. procedure cln   (var x: complex); { x:=ln(x)   }
  28. procedure cexp  (var x: complex); { x:=exp(x)  }
  29. procedure csin  (var x: complex); { x:=sin(x)  }
  30. procedure ccos  (var x: complex); { x:=cos(x)  }
  31. procedure ctan  (var x: complex); { x:=tan(x)  }
  32. procedure csinh (var x: complex); { x:=sinh(x) }
  33. procedure ccosh (var x: complex); { x:=cosh(x) }
  34. procedure ctanh (var x: complex); { x:=tanh(x) }
  35. procedure casin (var x: complex); { x:=asin(x) }
  36. procedure cacos (var x: complex); { x:=acos(x) }
  37. procedure catan (var x: complex); { x:=atan(x) }
  38. procedure casinh(var x: complex); { x:=asinh(x)}
  39. procedure cacosh(var x: complex); { x:=acosh(x)}
  40. procedure catanh(var x: complex); { x:=atanh(x)}
  41.  
  42. procedure cadd (var x,y: complex); { x:=x+y     }
  43. procedure csub (var x,y: complex); { x:=x-y     }
  44. procedure cmult(var x,y: complex); { x:=x*y     }
  45. procedure cdiv (var x,y: complex); { x:=x/y     }
  46. procedure cpow (var x,y: complex); { x:=pow(x,y)}
  47.  
  48. procedure cpoly (var x:complex; degree: integer; var coeffs);
  49.  
  50.  
  51. implementation
  52.  
  53. { "alphbeta" is used by the "cacos" and "casin" routines.  I
  54.   derived it from Abramowitz and Stegun, Handbook of Mathematical
  55.   Functions. (as well as most of the routines) }
  56.  
  57. procedure alphbeta(var x:complex); var z:xfloat;
  58.    begin with x do begin
  59.       z:=system.sqrt(sqr(rp+1)+sqr(ip));
  60.       rp:=system.sqrt(sqr(rp-1)+sqr(ip));
  61.       ip:=(z+rp)/2; rp:=ip-rp; ip:=ln(ip+sqrt(sqr(ip)-1)) end end;
  62.  
  63. function cabs; begin cabs:=hypot(x.rp,x.ip) end;
  64. procedure cacos; begin alphbeta(x); x.rp:=acos(x.rp) end;
  65. procedure cacosh; begin cacos(x); cjx(x) end;
  66. procedure cadd; begin x.rp:=x.rp+y.rp; x.ip:=x.ip+y.ip end;
  67. procedure casin; begin alphbeta(x); x.rp:=asin(x.rp) end;
  68. procedure casinh; begin cjx(x); casin(x); cnjx(x) end;
  69.  
  70. procedure catan; var rp2,ip2,z:xfloat; begin
  71.    rp2:=sqr(x.rp); ip2:=sqr(x.ip); z:=atan(2*x.rp/(1-rp2-rp2))/2;
  72.    x.ip:=ln((rp2+ip2+2*x.ip+1)/(rp2+ip2-2*x.ip+1))/4;
  73.    x.rp:=z end;
  74.  
  75. procedure catanh; begin cjx(x); catan(x); cnjx(x) end;
  76.  
  77. procedure ccos; var z:xfloat; begin
  78.    z:=cos(x.rp)*cosh(x.ip);
  79.    x.ip:=-sin(x.rp)*sinh(x.ip); x.rp:=z end;
  80.  
  81. procedure ccosh; var z:xfloat; begin
  82.    z:=cosh(x.rp)*cos(x.ip);
  83.    x.ip:=sinh(x.rp)*sin(x.ip); x.rp:=z end;
  84.  
  85. procedure cdiv; var w,z:xfloat; begin
  86.    w:=sqr(y.rp)+sqr(y.ip); z:=(x.rp*y.rp+x.ip*y.ip)/w;
  87.    x.ip:=(x.ip*y.rp-x.rp*y.ip)/w; x.rp:=z end;
  88.  
  89. procedure cexp; var z:xfloat; begin
  90.    z:=exp(x.rp); x.rp:=z*cos(x.ip); x.ip:=z*sin(x.ip) end;
  91.  
  92. procedure cjx; var z:xfloat; begin z:=-x.ip; x.ip:=x.rp; x.rp:=z end;
  93.  
  94. procedure cln; var z:xfloat; begin
  95.    z:=ln(cabs(x)); x.ip:=phase(x); x.rp:=z end;
  96.  
  97. procedure cmult; var z:xfloat; begin
  98.    z:=x.rp*y.rp-x.ip*y.ip; x.ip:=x.rp*y.ip+x.ip*y.rp; x.rp:=z end;
  99.  
  100. procedure cneg; begin x.rp:=-x.rp; x.ip:=-x.ip end;
  101. procedure cnjx; var z:xfloat; begin z:=x.ip; x.ip:=-x.rp; x.rp:=z end;
  102. procedure conj; begin x.ip:=-x.ip end;
  103.  
  104. procedure cpoly;
  105.    { coeffs is assumed to be an array [0..degree] of "complex" }
  106.    var y:complex;
  107.    begin y.rp:=0; y.ip:=0; while degree>=0 do begin
  108.       cmult(y,x); cadd(y,cxray(coeffs)[degree]); dec(degree) end;
  109.    x:=y end;
  110.  
  111. procedure cpow; begin cln(x); cmult(x,y); cexp(x) end;
  112.  
  113. procedure crecip; var w:xfloat; begin
  114.    w:=sqr(x.rp)+sqr(x.ip); x.rp:=x.rp/w; x.ip:=-x.ip/w end;
  115.  
  116. procedure csin; var z:xfloat; begin
  117.    z:=sin(x.rp)*cosh(x.ip); x.ip:=cos(x.rp)*sinh(x.ip); x.rp:=z end;
  118.  
  119. procedure csinh; var z:xfloat; begin
  120.    z:=sinh(x.rp)*cos(x.ip); x.ip:=cosh(x.rp)*sin(x.ip); x.rp:=z end;
  121.  
  122. procedure csqr; var z:xfloat; begin
  123.    z:=sqr(x.rp)-sqr(x.ip); x.ip:=2*x.rp*x.ip; x.rp:=z end;
  124.  
  125. procedure csqrt; var absx,z:xfloat; begin
  126.    absx:=cabs(x); z:=sqrt((absx+x.rp)/2);
  127.    x.ip:=sgn(x.ip)*sqrt((absx-x.rp)/2); x.rp:=z end;
  128.  
  129. procedure csub; begin x.rp:=x.rp-y.rp; x.ip:=x.ip-y.ip end;
  130.  
  131. procedure ctan; var z:xfloat; begin
  132.    z:=cos(2*x.rp)+cosh(2*x.ip);
  133.    x.rp:=sin(2*x.rp)/z; x.ip:=sinh(2*x.rp)/z end;
  134.  
  135. procedure ctanh; var z:xfloat; begin with x do begin
  136.    z:=cosh(2*rp)+cos(2*ip); rp:=sinh(2*rp)/z; ip:=sin(2*ip)/z
  137.    end end;
  138.  
  139. function phase; var z:xfloat; begin z:=atan2(x.rp,x.ip) end;
  140.  
  141. end.
  142.