complex.cpp 6.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318
  1. //#include "stdafx.h"
  2. #include <stdlib.h>
  3. #include <math.h>
  4. #include "complex.h"
  5. /*********************** constructors *****************************/
  6. Complex::Complex(double r, double i)
  7. : re(r), im(i)
  8. {}
  9. Complex::Complex(const Complex& c) {
  10. re = c.re;
  11. im = c.im;
  12. }
  13. /*********************** members *********************************/
  14. /********************* friends *****************************/
  15. Complex operator*(const double& r, const Complex& c) {
  16. Complex com;
  17. com.re = r * c.re;
  18. com.im = r * c.im;
  19. return(com);
  20. }
  21. Complex operator+(const double& r, const Complex& c) {
  22. Complex com;
  23. com.re = r + c.re;
  24. com.im = c.im;
  25. return(com);
  26. }
  27. Complex operator-(const double& r, const Complex& c) {
  28. Complex com;
  29. com.re = r - c.re;
  30. com.im = - c.im;
  31. return(com);
  32. }
  33. Complex operator/(const double& r, const Complex& c) {
  34. Complex com;
  35. double t, d;
  36. if(fabs(c.re) <= fabs(c.im)) {
  37. t = c.re / c.im;
  38. d = c.im * (1 + t*t);
  39. com.re = r * t / d;
  40. com.im = -r / d;
  41. } else {
  42. t = c.im / c.re;
  43. d = c.re * (1 + t*t);
  44. com.re = r / d;
  45. com.im = -r * t / d;
  46. }
  47. return(com);
  48. }
  49. double real(const Complex& c) {
  50. return(c.re);
  51. }
  52. double imag(const Complex& c) {
  53. return(c.im);
  54. }
  55. double arg(const Complex& c) {
  56. return(c == 0 ? 0 : atan2(c.im, c.re));
  57. }
  58. double cabs(const Complex& c) {
  59. double x, y, ans, temp;
  60. x = fabs(c.re);
  61. y = fabs(c.im);
  62. if (x == 0.0)
  63. ans = y;
  64. else if (y == 0.0)
  65. ans = x;
  66. else if (x > y) {
  67. temp = y/x;
  68. ans = x*sqrt(1.0+temp*temp);
  69. } else {
  70. temp = x/y;
  71. ans = y*sqrt(1.0+temp*temp);
  72. }
  73. return ans;
  74. }
  75. double cabsm(const Complex& c) {
  76. return(sqrt(c.re*c.re + c.im*c.im));
  77. }
  78. Complex conjg(const Complex& c) {
  79. Complex com;
  80. com.re = c.re;
  81. com.im = -c.im;
  82. return(com);
  83. }
  84. Complex cmplx(const double re,const double im) {
  85. Complex com;
  86. com.re = re;
  87. com.im = im;
  88. return(com);
  89. }
  90. Complex polar(const double& r, const double& theta) {
  91. return(r*cmplx(cos(theta),sin(theta)));
  92. }
  93. /****************** cosine of Complex number c *******************/
  94. /*** cos(c) = cos(c.re)*cosh(c.im)-i*sin(c.re)*sinh(c.im) **/
  95. Complex cos(const Complex& c) {
  96. Complex com;
  97. com = cmplx(cos(c.re)*cosh(c.im),-sin(c.re)*sinh(c.im));
  98. return(com);
  99. }
  100. /************* hyperbolic cosine of Complex number c *************/
  101. /** cosh(c) = cosh(c.re)*cos(c.im)+i*sinh(c.re)*sin(c.im) **/
  102. Complex cosh(const Complex& c) {
  103. Complex com;
  104. com = cmplx(cosh(c.re)*cos(c.im),sinh(c.re)*sin(c.im));
  105. return(com);
  106. }
  107. /******************* sine of Complex number c ********************/
  108. /*** sin(c) = sin(c.re)*cosh*c.im)+i*cos(c.re)*sinh(c.im) **/
  109. Complex sin(const Complex& c) {
  110. Complex com;
  111. com = cmplx(sin(c.re)*cosh(c.im),cos(c.re)*sinh(c.im));
  112. return(com);
  113. }
  114. /************** hyperbolic sine of Complex number c **************/
  115. /** sinh(c) = sinh(c.re)*cos(c.im)+i*cosh(c.re)*sin(c.im) **/
  116. Complex sinh(const Complex& c) {
  117. Complex com;
  118. com = cmplx(sinh(c.re)*cos(c.im),cosh(c.re)*sin(c.im));
  119. return(com);
  120. }
  121. /******************** tan of Complex number c ********************/
  122. Complex tan(const Complex& c) {
  123. return(sin(c)/cos(c));
  124. }
  125. /*************** hyperbolic tan of Complex number c **************/
  126. Complex tanh(const Complex& c) {
  127. return(sinh(c)/cosh(c));
  128. }
  129. /*************** square root of Complex number c *****************/
  130. Complex sqrt(const Complex& c) {
  131. double r = cabs(c);
  132. double nr, ni;
  133. if(r == 0.0) nr = ni = r;
  134. else if(c.re > 0) {
  135. nr = sqrt(0.5*(r+c.re));
  136. ni = c.im/nr/2.0;
  137. }
  138. else {
  139. ni = sqrt(0.5*(r-c.re));
  140. if(c.im < 0) ni = -ni;
  141. nr = c.im/ni/2.0;
  142. }
  143. return(cmplx(nr,ni));
  144. }
  145. /*************** square root of Complex number c *****************/
  146. /*** sqrtm(c) = sqrt(abs(c))*(cos(arg/2) + i*sin(arg/2)) ***/
  147. Complex sqrtm(const Complex& c) {
  148. Complex com;
  149. double mod, theta;
  150. mod = cabs(c);
  151. theta = arg(c)/2.0; /*****argument [-Pi, +Pi]*********/
  152. com = sqrt(mod) * cmplx(cos(theta),sin(theta));
  153. return(com);
  154. }
  155. /*************** exponential of Complex number c *****************/
  156. /******* exp(c) = exp(c.re)*(cos(c.re)+i*sin(c.im)) ********/
  157. Complex exp(const Complex& c) {
  158. Complex com;
  159. com = exp(c.re) * cmplx(cos(c.im),sin(c.im));
  160. return(com);
  161. }
  162. /*************** natural log of Complex number c *****************/
  163. /************** log(c) = log(abs(c))+i*arg(c)***************/
  164. Complex log(const Complex& c) {
  165. Complex com;
  166. int cs;
  167. if(c == 0) cs = 1;
  168. else cs = 2;
  169. switch(cs) {
  170. case 1: cerr << "Problem with log !!!" << endl; com = 0.0; break;
  171. case 2: com = cmplx(log(cabs(c)),arg(c)); break;
  172. default: cerr << "Should nevere get here!" << "\n";
  173. }
  174. return(com);
  175. }
  176. /****************** log10 of Complex number c ********************/
  177. Complex log10(const Complex& c) {
  178. return(log10(exp(1.0))*log(c));
  179. }
  180. /********* Complex number c raised to a Complex power y **********/
  181. Complex pow(const Complex& c, const Complex& y) {
  182. if(c == 0) {
  183. if(y == 0) return(cmplx(1,0));
  184. else return(cmplx(0,0));
  185. }
  186. double logr = log(cabs(c));
  187. double t = arg(c);
  188. return(polar(exp(logr*y.re-t*y.im),logr*y.im+t*y.re));
  189. }
  190. /********* Complex number c raised to a double power y **********/
  191. /***************** c**y = exp(y * log(c)) ******************/
  192. Complex pow(const Complex& c, const double& y) {
  193. if(c == 0) {
  194. if(y == 0) return(cmplx(1,0));
  195. else return(cmplx(0,0));
  196. }
  197. return(exp(y*log(c)));
  198. }
  199. /******** Complex number c raised to an integer power y *********/
  200. /**** c**y = abs(c)**y * (cos(y*theta)+i*sin(y*theta)) *****/
  201. Complex pow(const Complex& c, int& y) {
  202. if(c == 0) {
  203. if(y == 0) return(cmplx(1,0));
  204. else return(cmplx(0,0));
  205. }
  206. if(c.im == 0) {
  207. if(c.re < 0) return(pow(c,cmplx(y,0)));
  208. else return(cmplx(pow(c.re,double(y)),0.0));
  209. }
  210. register double r = pow(cabs(c),double(y));
  211. register double th = y*arg(c);
  212. return(cmplx(r*cos(th),r*sin(th)));
  213. }
  214. /********* double number y raised to a Complex power c **********/
  215. /****************** y**c = exp(c * log(y)) *****************/
  216. Complex pow(const double& y, const Complex& c) {
  217. if(y == 0) {
  218. if(c == 0) return(cmplx(1,0));
  219. else return(cmplx(0,0));
  220. }
  221. if(y < 0) return(pow(cmplx(y,0),c));
  222. if(c.im == 0.0) return(cmplx(pow(y,c.re),0.0));
  223. return(exp(c*log(y)));
  224. }
  225. istream& operator>>(istream& in, Complex& c) {
  226. in >> c.re >> c.im;
  227. return(in);
  228. }
  229. ostream& operator<<(ostream& out, Complex& c) {
  230. out << " ( " << c.re << ',' << " " << c.im << " )";
  231. return(out);
  232. }