home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / zpp / complex.cpp next >
Encoding:
C/C++ Source or Header  |  1992-04-03  |  8.5 KB  |  434 lines
  1. /* -------------------------------------------------------------------- */
  2. /* complex.cpp                                */
  3. /*                                    */
  4. /* Z++ Version 1.0                Last revised 04/03/92    */
  5. /*                                    */
  6. /* Complex number class for Turbo C++/Borland C++.            */
  7. /* Copyright 1992 by Carl W. Moreland                    */
  8. /* This source code may be freely distributed as long as the copyright    */
  9. /* notice remains intact.                        */
  10. /* -------------------------------------------------------------------- */
  11.  
  12. #include "complex.h"
  13.  
  14. unsigned char complex::zArgMode   = Z_RADIANS;
  15. unsigned char complex::zPrintMode = Z_COMMA;
  16. unsigned char complex::zLetter    = 'i';
  17.  
  18. /* ----- Constructors ------------------------------------------------- */
  19.  
  20. complex::complex(void)
  21. {
  22.   re = 0;
  23.   im = 0;
  24. }
  25.  
  26. complex::complex(const double real, const double imag)
  27. {
  28.   re = real;
  29.   im = imag;
  30. }
  31.  
  32. complex::complex(const complex& z)
  33. {
  34.   re = z.re;
  35.   im = z.im;
  36. }
  37.  
  38. /* -------------------------------------------------------------------- */
  39.  
  40. double re(const complex& z)        // friend version
  41. {
  42.   return z.re;
  43. }
  44.  
  45. double im(const complex& z)        // friend version
  46. {
  47.   return z.im;
  48. }
  49.  
  50. double real(const complex& z)        // friend version
  51. {
  52.   return z.re;
  53. }
  54.  
  55. double imag(const complex& z)        // friend version
  56. {
  57.   return z.im;
  58. }
  59.  
  60. double mag(const complex& z)        // magnitude |z|
  61. {
  62.   return sqrt(z.re*z.re + z.im*z.im);
  63. }
  64.  
  65. double arg(const complex& z)        // argument (angle)
  66. {
  67.   if(z.re == 0 && z.im == 0)        // this is actually a domain error
  68.     return 0;
  69.   if(complex::zArgMode == Z_RADIANS)
  70.     return atan2(z.im, z.re);
  71.   return atan2(z.im, z.re)/M_PI*180;
  72. }
  73.  
  74. complex conj(const complex& z)        // complex conjugate
  75. {
  76.   return complex(z.re, -z.im);
  77. }
  78.  
  79. double norm(const complex& z)
  80. {
  81.   return z.re*z.re + z.im*z.im;
  82. }
  83.  
  84. complex ptor(double mag, double angle)  // polar-to-rectangular
  85. {
  86.   if(complex::zArgMode == Z_RADIANS)
  87.     return complex(mag*cos(angle), mag*sin(angle));
  88.   return complex(mag*cos(angle/180*M_PI), mag*sin(angle/180*M_PI));
  89. }
  90.  
  91. complex rtop(double x, double y)    // rectangular-to-polar
  92. {
  93.   if(x == 0 && y == 0)            // this is actually a domain error
  94.     return Z0;
  95.   if(complex::zArgMode == Z_RADIANS)
  96.     return complex(sqrt(x*x + y*y), atan2(y, x));
  97.   return complex(sqrt(x*x + y*y), atan2(y, x)*180/M_PI);
  98. }
  99.  
  100. complex& complex::topolar(void)
  101. {
  102.   double re_tmp = re;
  103.  
  104.   if(re != 0 || im != 0)        // z = (0,0) is a domain error
  105.   {
  106.     re = sqrt(re*re + im*im);
  107.     im = atan2(im, re_tmp);
  108.   }
  109.  
  110.   if(complex::zArgMode == Z_DEGREES)
  111.     im *= (180/M_PI);
  112.  
  113.   return *this;
  114. }
  115.  
  116. complex& complex::torect(void)
  117. {
  118.   double re_tmp = re;
  119.  
  120.   re = re_tmp*cos(im);
  121.   im = re_tmp*sin(im);
  122.  
  123.   return *this;
  124. }
  125.  
  126. /* ----- Operators ---------------------------------------------------- */
  127.  
  128. void complex::operator=(const complex& z)
  129. {
  130.   re = z.re;
  131.   im = z.im;
  132. }
  133.  
  134. complex& complex::operator+=(const complex& z)
  135. {
  136.   re += z.re;
  137.   im += z.im;
  138.   return *this;
  139. }
  140.  
  141. complex& complex::operator-=(const complex& z)
  142. {
  143.   re -= z.re;
  144.   im -= z.im;
  145.   return *this;
  146. }
  147.  
  148. complex& complex::operator*=(const complex& z)
  149. {
  150.   *this = *this * z;
  151.   return *this;
  152. }
  153.  
  154. complex& complex::operator/=(const complex& z)
  155. {
  156.   *this = *this / z;
  157.   return *this;
  158. }
  159.  
  160. complex complex::operator+() const
  161. {
  162.   return *this;
  163. }
  164.  
  165. complex complex::operator-() const
  166. {
  167.   return complex(-re, -im);
  168. }
  169.  
  170. complex operator+(const complex& z1, const complex& z2)
  171. {
  172.   return complex(z1.re + z2.re, z1.im + z2.im);
  173. }
  174.  
  175. complex operator+(const complex& z, const double x)
  176. {
  177.   return complex(z.re+x, z.im);
  178. }
  179.  
  180. complex operator+(const double x, const complex& z)
  181. {
  182.   return complex(z.re+x, z.im);
  183. }
  184.  
  185. complex operator-(const complex& z1, const complex& z2)
  186. {
  187.   return complex(z1.re - z2.re, z1.im - z2.im);
  188. }
  189.  
  190. complex operator-(const complex& z, const double x)
  191. {
  192.   return complex(z.re-x, z.im);
  193. }
  194.  
  195. complex operator-(const double x, const complex& z)
  196. {
  197.   return complex(x-z.re, -z.im);
  198. }
  199.  
  200. complex operator*(const complex& z1, const complex& z2)
  201. {
  202.   double re = z1.re*z2.re - z1.im*z2.im;
  203.   double im = z1.re*z2.im + z1.im*z2.re;
  204.   return complex(re, im);
  205. }
  206.  
  207. complex operator*(const complex& z, const double x)
  208. {
  209.   return complex(z.re*x, z.im*x);
  210. }
  211.  
  212. complex operator*(const double x, const complex& z)
  213. {
  214.   return complex(z.re*x, z.im*x);
  215. }
  216.  
  217. complex operator/(const complex& z1, const complex& z2)
  218. {
  219.   if(z2 == Z0)
  220.     return complex(Zinf);        // z2 = Z0 is an error!
  221.  
  222.   double denom = z2.re*z2.re + z2.im*z2.im;
  223.   double re = (z1.re*z2.re + z1.im*z2.im)/denom;
  224.   double im = (z2.re*z1.im - z2.im*z1.re)/denom;
  225.   return complex(re, im);
  226. }
  227.  
  228. complex operator/(const complex& z, const double x)
  229. {
  230.   return complex(z.re/x, z.im/x);
  231. }
  232.  
  233. complex operator/(const double x, const complex& z)
  234. {
  235.   if(z == Z0)
  236.     return complex(Zinf);        // z = Z0 is an error!
  237.  
  238.   double denom = z.re*z.re + z.im*z.im;
  239.   return complex(x*z.re/denom, -z.im*x/denom);
  240. }
  241.  
  242. complex operator^(const complex& z1, const complex& z2)
  243. {
  244.   return pow(z1, z2);
  245. }
  246.  
  247. int operator==(const complex& z1, const complex& z2)
  248. {
  249.   return (z1.re == z2.re) && (z1.im == z2.im);
  250. }
  251.  
  252. int operator!=(const complex& z1, const complex& z2)
  253. {
  254.   return (z1.re != z2.re) || (z1.im != z2.im);
  255. }
  256.  
  257. /* ----- Math functions ----------------------------------------------- */
  258.  
  259. double abs(const complex& z)
  260. {
  261.   return sqrt(z.re*z.re + z.im*z.im);
  262. }
  263.  
  264. complex sqrt(const complex& z)
  265. {
  266.   return ptor(sqrt(abs(z)), arg(z)/2);
  267. }
  268.  
  269. complex pow(const complex& base, const double exponent)
  270. {
  271.   if(base != Z0 && exponent == 0.0)
  272.     return complex(1,0);
  273.  
  274.   if (base == Z0 && exponent > 0)
  275.     return Z0;
  276.  
  277.   // base == Z0 && exponent == 0 is undefined!
  278.  
  279.   return ptor(pow(abs(base), exponent), exponent*arg(base));
  280. }
  281.  
  282. complex pow(const double base, const complex& exponent)
  283. {
  284.   if(base != 0.0 && exponent == Z0)
  285.     return complex(1,0);
  286.  
  287.   if (base == 0 && re(exponent) > 0)
  288.     return complex(0,0);
  289.  
  290.   // base == 0 && re(exponent) == 0 is undefined!
  291.  
  292.   if(base > 0.0)
  293.     return exp(exponent * log(fabs(base)));
  294.  
  295.   return exp(exponent * complex(log(fabs(base)), M_PI));
  296. }
  297.  
  298. complex pow(const complex& base, const complex& exponent)
  299. {
  300.   if(base != Z0 && exponent == Z0)
  301.     return complex(1,0);
  302.  
  303.   if(base == Z0 && re(exponent) > 0)
  304.     return complex(0,0);
  305.  
  306.   // base == Z0 && re(exponent) == 0 is undefined!
  307.  
  308.   return exp(exponent * log(base));
  309. }
  310.  
  311. /* ----- Trig functions ----------------------------------------------- */
  312.  
  313. complex exp(const complex& z)
  314. {
  315.   double mag = exp(z.re);
  316.   return complex(mag*cos(z.im), mag*sin(z.im));
  317. }
  318.  
  319. complex log(const complex& z)
  320. {
  321.   return complex(log(mag(z)), atan2(z.im, z.re));
  322. }
  323.  
  324. complex log10(const complex& z)
  325. {
  326.   return log(z) * M_LOG10E;
  327. }
  328.  
  329. complex sin(const complex& z)
  330. {
  331.   return (exp(Zi*z) - exp(-Zi*z))/(2*Zi);
  332. }
  333.  
  334. complex cos(const complex& z)
  335. {
  336.   return (exp(Zi*z) + exp(-Zi*z))/2;
  337. }
  338.  
  339. complex tan(const complex& z)
  340. {
  341.   return sin(z)/cos(z);
  342. }
  343.  
  344. complex asin(const complex& z)
  345. {
  346.   return -Zi*log(Zi*z+sqrt(1-z*z));
  347. }
  348.  
  349. complex acos(const complex& z)
  350. {
  351.   return -Zi*log(z+Zi*sqrt(1-z*z));
  352. }
  353.  
  354. complex atan(const complex& z)
  355. {
  356.   return -0.5*Zi * log((1+Zi*z)/(1-Zi*z));
  357. }
  358.  
  359. complex sinh(const complex& z)
  360. {
  361.   return (exp(z) - exp(-z))/2;
  362. }
  363.  
  364. complex cosh(const complex& z)
  365. {
  366.   return (exp(z) + exp(-z))/2;
  367. }
  368.  
  369. complex tanh(const complex& z)
  370. {
  371.   return sinh(z)/cosh(z);
  372. }
  373.  
  374. /* ----- Misc functions ----------------------------------------------- */
  375.  
  376. void complex::SetArgMode(unsigned char mode) const
  377. {
  378.   if(mode == Z_RADIANS || mode == Z_DEGREES)
  379.     zArgMode = mode;
  380. }
  381.  
  382. void complex::SetPrintMode(unsigned char mode) const
  383. {
  384.   if(mode == Z_COMMA || mode == Z_LETTER)
  385.     zPrintMode = mode;
  386. }
  387.  
  388. void complex::SetLetter(unsigned char letter) const
  389. {
  390.   zLetter = letter;
  391. }
  392.  
  393. /* ----- Stream I/O --------------------------------------------------- */
  394.  
  395. ostream& operator<<(ostream& s, const complex& z)
  396. {
  397.   char sign[] = "   ";
  398.  
  399.   if(complex::zPrintMode == Z_COMMA)
  400.     return s << "(" << z.re << ", " << z.im << ")";
  401.  
  402.   if(z.im == 0 || z.im/fabs(z.im) == 1)
  403.     sign[1] = '+';
  404.   else
  405.     sign[1] = '-';
  406.   return s << z.re << sign << complex::zLetter << fabs(z.im);
  407.  
  408. }
  409.  
  410. istream& operator>>(istream& s, complex& z)
  411. {
  412.   char ch;
  413.   double real=0, imag=0;
  414.  
  415.   s >> ch;
  416.   if(ch == '(')
  417.   {
  418.     s >> real >> ch;
  419.     if(ch == ',')
  420.       s >> imag >> ch;
  421.     if(ch != ')')
  422.       s.clear(ios::badbit | s.rdstate());
  423.   }
  424.   else
  425.   {
  426.     s.putback(ch);
  427.     s >> real;
  428.   }
  429.  
  430.   z = complex(real, imag);
  431.   return s;
  432. }
  433.  
  434.