home *** CD-ROM | disk | FTP | other *** search
/ GEMini Atari / GEMini_Atari_CD-ROM_Walnut_Creek_December_1993.iso / files / gnu / g__lib / complex.cc < prev    next >
Encoding:
C/C++ Source or Header  |  1993-07-23  |  4.7 KB  |  245 lines
  1. /* 
  2. Copyright (C) 1988 Free Software Foundation
  3.     written by Doug Lea (dl@rocky.oswego.edu)
  4.  
  5. This file is part of GNU CC.
  6.  
  7. GNU CC is distributed in the hope that it will be useful,
  8. but WITHOUT ANY WARRANTY.  No author or distributor
  9. accepts responsibility to anyone for the consequences of using it
  10. or for whether it serves any particular purpose or works at all,
  11. unless he says so in writing.  Refer to the GNU CC General Public
  12. License for full details.
  13.  
  14. Everyone is granted permission to copy, modify and redistribute
  15. GNU CC, but only under the conditions described in the
  16. GNU CC General Public License.   A copy of this license is
  17. supposed to have been given to you along with GNU CC so you
  18. can know your rights and responsibilities.  It should be in a
  19. file named COPYING.  Among other things, the copyright notice
  20. and this notice must be preserved on all copies.  
  21. */
  22.  
  23. #include <Complex.h>
  24. #include <std.h>
  25.  
  26. // error handling
  27.  
  28. void default_Complex_error_handler(char* msg)
  29. {
  30.   cerr << "Fatal Complex arithmetic error. " << msg << "\n";
  31.   exit(1);
  32. }
  33.  
  34. one_arg_error_handler_t Complex_error_handler = default_Complex_error_handler;
  35.  
  36. one_arg_error_handler_t set_Complex_error_handler(one_arg_error_handler_t f)
  37. {
  38.   one_arg_error_handler_t old = Complex_error_handler;
  39.   Complex_error_handler = f;
  40.   return old;
  41. }
  42.  
  43. void Complex::error(char* msg)
  44. {
  45.   (*Complex_error_handler)(msg);
  46. }
  47.  
  48. Complex operator / (Complex& x, Complex& y)
  49. {
  50.   double den = norm(y);
  51.   if (den == 0.0) x.error ("Attempted division by zero.");
  52.   return Complex((x.re * y.re + x.im * y.im) / den, 
  53.                  (x.im * y.re - x.re * y.im) / den);
  54. }
  55.  
  56. Complex operator / (double x, Complex& y)
  57. {
  58.   double den = norm(y);
  59.   if (den == 0.0) y.error ("Attempted division by zero.");
  60.   return Complex((x * y.re) / den, -(x * y.im) / den);
  61. }
  62.  
  63. Complex operator / (Complex& x, double y)
  64. {
  65.   if (y == 0.0) x.error ("Attempted division by zero.");
  66.   return Complex(x.re / y, x.im / y);
  67. }
  68.  
  69. Complex& Complex::operator /= (Complex& y)
  70. {
  71.   double den = norm(y);
  72.   if (den == 0.0) error ("Attempted division by zero.");
  73.   double r = (re * y.re + im * y.im) / den;
  74.   im = (im * y.re - re * y.im) / den;
  75.   re = r;
  76.   return *this;
  77. }
  78.  
  79. Complex& Complex::operator /= (double y)
  80. {
  81.   if (y == 0.0) error ("Attempted division by zero.");
  82.   re /= y;
  83.   im /= y;
  84.   return *this;
  85. }
  86.  
  87. Complex polar(double r, double t = 0.0)
  88. {
  89.   return Complex(r * cos(t), r * sin(t));
  90. }
  91.  
  92. Complex exp(Complex& x)
  93. {
  94.   double r = exp(x.re);
  95.   return Complex(r * cos(x.im), r * sin(x.im));
  96. }
  97.  
  98. Complex cosh(Complex& x)
  99. {
  100.   return Complex(cos(x.im) * cosh(x.re), sin(x.im) * sinh(x.re));
  101. }
  102.  
  103. Complex sinh(Complex& x)
  104. {
  105.   return Complex(cos(x.im) * sinh(x.re), sin(x.im) * cosh(x.re));
  106. }
  107.  
  108. Complex cos(Complex& x)
  109. {
  110.   return Complex(cos(x.re) * cosh(x.im), -sin(x.re) * sinh(x.im));
  111. }
  112.  
  113. Complex sin(Complex& x)
  114. {
  115.   return Complex(sin(x.re) * cosh(x.im), cos(x.re) * sinh(x.im));
  116. }
  117.  
  118. Complex log(Complex& x)
  119. {
  120.   double h = hypot(x.re, x.im);
  121.   if (h <= 0.0) x.error("attempted log of zero magnitude number.");
  122.   return Complex(log(h), atan2(x.im, x.re));
  123. }
  124.  
  125. Complex pow(Complex& x, Complex& p)
  126. {
  127.   Complex res;
  128.   if (p.re == 0.0 && p.im == 0.0)
  129.   {
  130.     res.re = 1.0; res.im = 0.0;
  131.   }
  132.   else if (x.re == 0.0 && x.im == 0.0)
  133.   {
  134.     res.re = 0.0; res.im = 0.0;
  135.   }
  136.   else
  137.   {
  138.     double h = hypot(x.re, x.im);
  139.     if (h <= 0.0) x.error("attempted power of zero magnitude number.");
  140.     double lr = exp(log(h) * p.re);
  141.     double li = atan2(x.im, x.re) * p.im;
  142.     res.re = lr * cos(li);
  143.     res.im = lr * sin(li);
  144.   }
  145.   return res;
  146. }
  147.  
  148. Complex sqrt(Complex& x)
  149. {
  150.   double r, i;
  151.   if (x.re == 0.0 && x.im == 0.0)
  152.     r = i = 0.0;
  153.   else
  154.   {
  155.     double s = sqrt((abs(x.re) + hypot(x.re, x.im)) * 0.5);
  156.     double d = (x.im / s) * 0.5;
  157.     if (x.re > 0.0)
  158.     {
  159.       r = s;
  160.       i = d;
  161.     }
  162.     else if (x.im >= 0.0)
  163.     {
  164.       r = d;
  165.       i = s;
  166.     }
  167.     else
  168.     {
  169.       r = -d;
  170.       i = -s;
  171.     }
  172.   }
  173.   return Complex(r, i);
  174. }
  175.  
  176.  
  177. Complex pow(Complex& x, long p)
  178. {
  179.   Complex res;
  180.   if (p == 0)
  181.   {
  182.     res.re = 1.0; res.im = 0.0;
  183.   }
  184.   else if (x == 0.0)
  185.   {
  186.     res.re = 0.0; res.im = 0.0;
  187.   }
  188.   else
  189.   {
  190.     res.re = 1.0; res.im = 0.0;
  191.     Complex b = x;
  192.     if (p < 0)
  193.     {
  194.       p = -p;
  195.       b = 1.0 / b;
  196.     }
  197.     for(;;)
  198.     {
  199.       if (p & 1)
  200.         res *= b;
  201.       if ((p >>= 1) == 0)
  202.         break;
  203.       else
  204.         b *= b;
  205.     }
  206.   }
  207.   return res;
  208. }
  209.  
  210. ostream& operator << (ostream& s, Complex& x)
  211. {
  212.   return s << "(" << x.re << ", " << x.im << ")" ;
  213. }
  214.  
  215. istream& operator >> (istream& s, Complex& x)
  216. {
  217.   char ch;
  218.   s >> WS;
  219.   s.get(ch);
  220.   if (ch == '(')
  221.   {
  222.     s >> x.re;
  223.     s >> WS;
  224.     s.get(ch);
  225.     if (ch == ',')
  226.     {
  227.       s >> x.im;
  228.       s >> WS;
  229.       s .get(ch);
  230.     }
  231.     else
  232.       x.im = 0;
  233.     if (ch != ')')
  234.       s.clear(_bad);
  235.   }
  236.   else
  237.   {
  238.     s.unget(ch);
  239.     s >> x.re;
  240.     x.im = 0;
  241.   }
  242.   return s;
  243. }
  244.  
  245.