home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 7 / 07.iso / c / c065 / 2.ddi / CLIB2.ZIP / CPLX2.CPP < prev    next >
Encoding:
C/C++ Source or Header  |  1990-06-07  |  5.4 KB  |  221 lines
  1. /*-----------------------------------------------------------------------*
  2.  * filename - cplx2.cpp
  3.  * C++ Complex Library Routines
  4.  *-----------------------------------------------------------------------*/
  5.  
  6. /*[]------------------------------------------------------------[]*/
  7. /*|                                                              |*/
  8. /*|     Turbo C++ Run Time Library - Version 1.0                 |*/
  9. /*|                                                              |*/
  10. /*|                                                              |*/
  11. /*|     Copyright (c) 1990 by Borland International              |*/
  12. /*|     All Rights Reserved.                                     |*/
  13. /*|                                                              |*/
  14. /*[]------------------------------------------------------------[]*/
  15.  
  16. #include <complex.h>
  17.  
  18. static const complex _Cdecl complex_zero(0, 0);
  19. inline long double sqr(double x) { return x*x;}
  20.  
  21. #pragma warn -lvc
  22.  
  23. complex& _Cdecl complex::operator/=(complex& z2)
  24. {
  25.     long double sum_sqrs = norm(z2);
  26.     *this *= conj(z2);
  27.     re /= sum_sqrs;
  28.     im /= sum_sqrs;
  29.     return *this;
  30. }
  31.  
  32. complex _Cdecl operator/(complex& z1, complex& z2)
  33. {
  34.     return z1 * conj(z2) / norm(z2);
  35. }
  36.  
  37. complex _Cdecl operator/(double re_val1, complex& z2)
  38. {
  39.     return re_val1 * conj(z2) / norm(z2);
  40. }
  41.  
  42. complex _Cdecl acos(complex& z)
  43. {
  44.     // complex i(0,1);
  45.     // return -i*log(z+i*sqrt(1-z*z));
  46.     complex temp1(1 - sqr(z.re) + sqr(z.im), -2*z.re*z.im);
  47.     double phi = arg(temp1)/2;
  48.     double rp = sqrt(abs(temp1));
  49.     complex temp2(z.re - rp*sin(phi), z.im + rp*cos(phi));
  50.     return complex(arg(temp2), -log(abs(temp2)));
  51. }
  52.  
  53. double _Cdecl arg(complex& z)
  54. {
  55.     return z == complex_zero ? 0.0 : atan2(z.im, z.re);
  56. }
  57.  
  58. complex _Cdecl asin(complex& z)
  59. {
  60.     // complex i(0,1);
  61.     // return -i*log(i*z+sqrt(1-z*z));
  62.  
  63.     complex temp1(1 - sqr(z.re) + sqr(z.im), -2*z.re*z.im);
  64.     double phi = arg(temp1) / 2;
  65.     double rp = sqrt(abs(temp1));
  66.     complex temp2(-z.im + rp*cos(phi), z.re + rp*sin(phi));
  67.     return complex(arg(temp2), -log(abs(temp2)));
  68. }
  69.  
  70. complex _Cdecl atan(complex& z)
  71. {
  72.     // complex i(0, 1);
  73.     // return -0.5*i * log((1+i*z)/(1-i*z));
  74.     // if z=i then a floating-point exception may occur
  75.  
  76.     double opb = 1 + z.im;
  77.     double a2 = sqr(z.re);
  78.     double den = opb*opb + a2;
  79.     complex temp(((1-z.im)*opb - a2)/den, 2*z.re/den);
  80.     return complex(arg(temp)/2, -log(abs(temp))/2);
  81. }
  82.  
  83. complex _Cdecl cos(complex& z)
  84. {
  85.     // complex i(0, 1);
  86.     // return (exp(i*z) + exp(-i*z))/2;
  87.  
  88.     long double eb  = exp(z.im);
  89.     long double emb = 1 / eb;
  90.     return complex(cos(z.re)*(emb+eb)/2, sin(z.re)*(emb-eb)/2);
  91. }
  92.  
  93. complex _Cdecl cosh(complex& z)
  94. {
  95.     // return (exp(z) + exp(-z))/2;
  96.     long double ea = exp(z.re);
  97.     long double eainv = 1 / ea;
  98.     return complex(cos(z.im)*(ea + eainv)/2, sin(z.im)*(ea - eainv)/2);
  99. }
  100.  
  101. complex _Cdecl exp(complex& z)
  102. {
  103.     return polar(exp(z.re),z.im);
  104. }
  105.  
  106. complex _Cdecl log(complex& z)
  107. {
  108.     // exception if z = 0
  109.     return complex(log(abs(z)), arg(z));
  110. }
  111.  
  112. complex _Cdecl log10(complex& z)
  113. {
  114.     // return log10e*log(z);
  115.     return log(z) * M_LOG10E;
  116. }
  117.  
  118. complex _Cdecl pow(complex& base, double expon)
  119. {
  120.     // return exp(expon*log(base));
  121.  
  122.     if (base == complex_zero && expon > 0) return complex_zero;
  123.     return polar(pow(abs(base), expon), expon*arg(base));
  124. }
  125.  
  126. complex _Cdecl pow(double base, complex& expon)
  127. {
  128.     // return exp(expon*log(base));
  129.     if (base == 0 && real(expon) > 0) return complex_zero;
  130.     double lnx = log(fabs(base));
  131.     if (base > 0.0)
  132.         return exp(expon * lnx);
  133.     return exp(expon * complex(lnx,M_PI));
  134. }
  135.  
  136. complex _Cdecl pow(complex& base, complex& expon)
  137. {
  138.     if (base == complex_zero && real(expon) > 0) return complex_zero;
  139.     return exp(expon * log(base));
  140. }
  141.  
  142. complex _Cdecl sin(complex& z)
  143. {
  144.     // complex i(0,1);
  145.     // return (exp(i*z) - exp(-i*z))/(2*i);
  146.  
  147.     long double eb  = exp(z.im);
  148.     long double emb = 1 / eb;
  149.     return complex(sin(z.re)*(emb+eb)/2, -0.5*cos(z.re)*(emb-eb));
  150. }
  151.  
  152. complex _Cdecl sinh(complex& z)
  153. {
  154.     // return (exp(z) - exp(-z))/2
  155.     long double ea = exp(z.re);
  156.     long double eainv = 1 / ea;
  157.     return complex(cos(z.im)*(ea - eainv)/2, sin(z.im)*(ea + eainv)/2);
  158. }
  159.  
  160. complex _Cdecl sqrt(complex& z)
  161. {
  162.     // return pow(z, 0.5);
  163.     return polar(sqrt(abs(z)), arg(z) / 2);
  164. }
  165.  
  166. complex _Cdecl tan(complex& z)
  167. {
  168.     // return sin(z)/cos(z)
  169.     double sina = sin(z.re);
  170.     double cosa = cos(z.re);
  171.     long double eb = exp(z.im);
  172.     long double emb = 1 / eb;
  173.     double emin = emb - eb;
  174.     double eplus = emb + eb;
  175.     return complex(4*sina*cosa, -emin*eplus)
  176.         / (sqr(cosa*eplus) + sqr(sina*emin));
  177. }
  178.  
  179. complex _Cdecl tanh(complex& z)
  180. {
  181.     // return sinh(z)/cosh(z);
  182.     double sinb = sin(z.im);
  183.     double cosb = cos(z.im);
  184.     long double ea = exp(z.re);
  185.     long double eainv = 1 / ea;
  186.     double eamin = ea - eainv;
  187.     double eapls = ea + eainv;
  188.     return complex(eamin*eapls, 4*sinb*cosb)
  189.         / (sqr(cosb*eapls) + sqr(sinb*eamin));
  190. }
  191.  
  192.  
  193. // Stream I/O function definitions
  194.  
  195. ostream& _Cdecl operator<<(ostream& s, complex& z)
  196. {
  197.     return s << "(" << real(z) << ", " << imag(z) << ")";
  198. }
  199.  
  200. istream& _Cdecl operator>>(istream& s, complex& z)
  201. {
  202.     double re_val = 0.0, im_val = 0.0;
  203.     char c = 0;
  204.  
  205.     s >> c;
  206.     if (c == '(') {
  207.         s >> re_val >> c;
  208.          if (c == ',') s >> im_val >> c;
  209.          if (c != ')') s.clear(ios::badbit | s.rdstate());
  210.     }
  211.     else {
  212.         s.putback(c);
  213.         s >> re_val;
  214.     }
  215.     if (s) z = complex(re_val, im_val);
  216.     return s;
  217. }
  218.  
  219. #pragma warn .lvc
  220.  
  221.