1 //This file is part of CRAB
2 //http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html
4 /*--------------------------------------------------
8 -- Some modification for speed
9 -- Some protections from overflow and underflow
10 -- Solving of equations of order 2,3,4:
11 Solve2, Solve3, Solve4;
12 --------------------------------------------------*/
14 #ifndef __COMPLEX_CXX__
15 #define __COMPLEX_CXX__
18 #define small_epsilon 1e-12
22 #define PI 3.141592653589793238462643
27 #include <Riostream.h>
31 /////////////////////////// Helper Error Function
33 void ComplexError (const char *msg)
35 cout << endl << msg << endl;
39 void ComplexWarning (const char *msg)
41 cout << endl << msg << endl;
49 Complex (double a=0, double b=0): Re(a),Im(b) {}
51 double GetRe () { return Re; }
52 double GetIm () { return Im; }
53 void SetRe (double a) { Re = a; }
54 void SetIm (double a) { Im = a; }
55 void Set (double a, double b) { Re = a; Im = b; }
56 void SetPolar (double, double);
59 Complex operator ! () const;
60 Complex operator ! ();
61 Complex operator + () const;
62 Complex operator + ();
63 Complex operator - () const;
64 Complex operator - ();
65 Complex operator += (const Complex &);
66 Complex operator += (double);
67 Complex operator -= (const Complex &);
68 Complex operator -= (double);
69 Complex operator *= (const Complex &);
70 Complex operator *= (double);
71 Complex operator /= (const Complex &);
72 Complex operator /= (double);
74 friend ostream& operator << (ostream&, const Complex &);
75 friend istream& operator >> (istream&, Complex &);
78 const Complex ImUnit (0,1);
80 void Zero (Complex &);
81 Complex Polar (double, double);
82 double real (const Complex &);
83 double imag (const Complex &);
84 double Arg (const Complex &);
85 double phase(const Complex &);
86 double abs (const Complex &);
87 void Solve2 (Complex*, const Complex&, const Complex&);
88 void Solve3 (Complex*, const Complex&, const Complex&, const Complex&);
89 void Solve4 (Complex*, const Complex&, const Complex&, const Complex&, const Complex&);
90 Complex Solve2 (const Complex&, const Complex&, int RootNumber);
91 Complex Solve3 (const Complex&, const Complex&, const Complex&, int RootNumber);
92 Complex Solve4 (const Complex&, const Complex&, const Complex&, const Complex&, int RootNumber);
93 Complex conj (const Complex &);
94 Complex sqr (const Complex &);
95 Complex sqrt (const Complex &, int=0);
96 Complex pow (Complex, int);
97 Complex pow (Complex, double);
98 Complex pow (Complex, Complex);
99 Complex root (const Complex &, int, int=0);
100 Complex exp (const Complex &);
101 Complex log (const Complex &);
102 Complex sin (const Complex &);
103 Complex cos (const Complex &);
104 Complex tan (const Complex &);
105 Complex cot (const Complex &);
106 Complex sec (const Complex &);
107 Complex csc (const Complex &);
108 Complex sinh (const Complex &);
109 Complex cosh (const Complex &);
110 Complex tanh (const Complex &);
111 Complex coth (const Complex &);
112 Complex sech (const Complex &);
113 Complex csch (const Complex &);
114 Complex asin (const Complex &, int=0);
115 Complex acos (const Complex &, int=0);
116 Complex atan (const Complex &);
117 Complex acot (const Complex &);
118 Complex asec (const Complex &, int=0);
119 Complex acsc (const Complex &, int=0);
120 Complex asinh(const Complex &, int=0);
121 Complex acosh(const Complex &, int=0);
122 Complex atanh(const Complex &);
123 Complex acoth(const Complex &);
124 Complex asech(const Complex &, int=0);
125 Complex acsch(const Complex &, int=0);
127 Complex operator ^ (const Complex &, int);
128 Complex operator ^ (const Complex &, double);
129 Complex operator ^ (const Complex &, const Complex &);
131 double operator += (double&, const Complex &);
132 double operator -= (double&, const Complex &);
133 double operator *= (double&, const Complex &);
134 double operator /= (double&, const Complex &);
136 Complex operator + (const Complex &, const Complex &);
137 Complex operator + (const Complex &, double);
138 Complex operator + (double, const Complex &);
139 Complex operator - (const Complex &, const Complex &);
140 Complex operator - (const Complex &, double);
141 Complex operator - (double, const Complex &);
142 Complex operator * (const Complex &, const Complex &);
143 Complex operator * (const Complex &, double);
144 Complex operator * (double, const Complex &);
145 Complex operator / (const Complex &, const Complex &);
146 Complex operator / (const Complex &, double);
147 Complex operator / (double, const Complex &);
149 int operator == (const Complex &, double);
150 int operator == (double, const Complex &);
151 int operator == (const Complex &, const Complex &);
152 int operator != (const Complex &, double);
153 int operator != (double, const Complex &);
154 int operator != (const Complex &, const Complex &);
155 int operator >= (const Complex &, double);
156 int operator >= (double, const Complex &);
157 int operator >= (const Complex &, const Complex &);
158 int operator <= (const Complex &, double);
159 int operator <= (double, const Complex &);
160 int operator <= (const Complex &, const Complex &);
161 int operator > (const Complex &, double);
162 int operator > (double, const Complex &);
163 int operator > (const Complex &, const Complex &);
164 int operator < (const Complex &, double);
165 int operator < (double, const Complex &);
166 int operator < (const Complex &, const Complex &);
168 ////////////////////////////////////////////////
171 inline Complex Complex::operator + () const
176 inline Complex Complex::operator + ()
181 inline Complex Complex::operator - () const
183 return Complex (-Re, -Im);
186 inline Complex Complex::operator - ()
188 return Complex (-Re, -Im);
191 inline Complex Complex::operator ! () const // complex conjugate //
193 return Complex (Re, -Im);
196 inline Complex Complex::operator ! () // complex conjugate //
198 return Complex (Re, -Im);
201 /////////////////////// +=,-=,*=,/=
203 double operator += (double &a, const Complex &b)
205 if (TMath::Abs(b.Im) < small_epsilon) a+=b.Re;
207 ComplexError ("Error in double+=Complex: Complex is not double number");
211 double operator -= (double &a, const Complex &b)
213 if (TMath::Abs(b.Im) < small_epsilon) a-=b.Re;
215 ComplexError ("Error in double -=Complex: Complex is not double number");
219 double operator *= (double &a, const Complex &b)
221 if (TMath::Abs(b.Im) < small_epsilon) a*=b.Re;
223 ComplexError ("Error in double *=Complex: Complex is not double number");
227 double operator /= (double &a, const Complex &b)
229 if (TMath::Abs(b.Im) < small_epsilon) a/=b.Re;
231 ComplexError ("Error in double /=Complex: Complex is not double number");
235 inline Complex Complex::operator += (const Complex &a)
242 inline Complex Complex::operator += (double a)
248 inline Complex Complex::operator -= (const Complex &a)
255 inline Complex Complex::operator -= (double a)
261 Complex Complex::operator *= (const Complex &a)
263 double t1=Re*a.Re, t2=Im*a.Im;
264 Im = (Re+Im) * (a.Re+a.Im) - t1 - t2;
269 inline Complex Complex::operator *= (double a)
276 Complex Complex::operator /= (const Complex &a)
279 if (TMath::Abs(a.Re) >= TMath::Abs(a.Im))
282 t2 = a.Re + a.Im * t1;
283 temp = (Re + Im * t1) / t2;
284 Im = (Im - Re * t1) / t2;
290 t2 = a.Re * t1 + a.Im;
291 temp = (Re * t1 + Im) / t2;
292 Im = (Im * t1 - Re) / t2;
298 inline Complex Complex::operator /= (double a)
305 /////////////////////////// +,-,*,/
307 inline Complex operator + (const Complex &a, const Complex &b)
309 return Complex (a.Re+b.Re, a.Im+b.Im);
312 inline Complex operator + (const Complex &a, double b)
314 return Complex (a.Re+b, a.Im);
317 inline Complex operator + (double b, const Complex &a)
319 return Complex (a.Re+b, a.Im);
322 inline Complex operator - (const Complex &a, const Complex &b)
324 return Complex (a.Re-b.Re, a.Im-b.Im);
327 inline Complex operator - (const Complex &a, double b)
329 return Complex (a.Re-b, a.Im);
332 inline Complex operator - (double b, const Complex &a)
334 return Complex (b-a.Re, -a.Im);
337 Complex operator * (const Complex &a, const Complex &b)
339 double t1 = a.Re * b.Re;
340 double t2 = a.Im * b.Im;
341 return Complex (t1 - t2, (a.Re+a.Im) * (b.Re+b.Im) - t1 - t2);
344 inline Complex operator * (const Complex &a, double b)
346 return Complex (a.Re*b, a.Im*b);
349 inline Complex operator * (double b, const Complex &a)
351 return Complex (a.Re*b, a.Im*b);
354 inline Complex operator / (const Complex &a, double b)
356 return Complex (a.Re/b, a.Im/b);
359 inline Complex operator / (const Complex &a, const Complex &b)
362 if (TMath::Abs(b.Re) >= TMath::Abs(b.Im))
365 t2 = b.Re + b.Im * t1;
366 return Complex ((a.Re + a.Im * t1) / t2, (a.Im - a.Re * t1) / t2);
371 t2 = b.Re * t1 + b.Im;
372 return Complex ((a.Re * t1 + a.Im) / t2, (a.Im * t1 - a.Re) / t2);
376 inline Complex operator / (double b, const Complex &a)
378 return (Complex (b,0)) / a;
381 //////////////////////// <<,>>
383 ostream& operator << (ostream &stream, const Complex &a)
385 stream<<" "<<a.Re<<" "<<a.Im<<" ";
389 istream& operator >> (istream &stream, Complex &a)
395 /////////////////////// ==,!=
397 inline int operator == (const Complex &a, double b)
399 return (a.Re==b) && (TMath::Abs(a.Im)<small_epsilon);
402 inline int operator == (double b, const Complex &a)
404 return (a.Re==b) && (TMath::Abs(a.Im)<small_epsilon);
407 inline int operator == (const Complex &a, const Complex &b)
409 return (a.Re==b.Re) && (a.Im==b.Im);
412 inline int operator != (const Complex &a, double b)
414 return (a.Re!=b) || (TMath::Abs(a.Im)>small_epsilon);
417 inline int operator != (double b, const Complex &a)
419 return (a.Re!=b) || (TMath::Abs(a.Im)>small_epsilon);
422 inline int operator != (const Complex &a, const Complex &b)
424 return (a.Re!=b.Re) || (a.Im!=b.Im);
427 /////////////////////// >=,<=
429 inline int operator >= (const Complex &a, double b)
434 inline int operator >= (double b, const Complex &a)
439 inline int operator >= (const Complex &a, const Complex &b)
441 return abs(a) >= abs(b);
444 inline int operator <= (const Complex &a, double b)
449 inline int operator <= (double b, const Complex &a)
454 inline int operator <= (const Complex &a, const Complex &b)
456 return abs(a) <= abs(b);
459 /////////////////////// >,<
461 inline int operator > (const Complex &a, double b)
466 inline int operator > (double b, const Complex &a)
471 inline int operator > (const Complex &a, const Complex &b)
473 return abs(a) > abs(b);
476 inline int operator < (const Complex &a, double b)
481 inline int operator < (double b, const Complex &a)
486 inline int operator < (const Complex &a, const Complex &b)
488 return abs(a) < abs(b);
491 ///////////////////////// Functions
492 double real (const Complex &a)
496 double imag (const Complex &a)
500 double abs (const Complex &a)
502 if (a.Im == 0) return TMath::Abs(a.Re);
503 if (a.Re == 0) return TMath::Abs(a.Im);
504 double R = TMath::Abs(a.Re), I = TMath::Abs(a.Im);
506 (R * sqrt (1 + a.Im*a.Im/a.Re/a.Re)) :
507 (I * sqrt (1 + a.Re*a.Re/a.Im/a.Im)) ;
510 double Arg (const Complex &a)
512 return atan2(a.Im,a.Re);
515 double phase (const Complex &a)
517 return atan2(a.Im,a.Re);
520 inline void Zero (Complex &a) { a.Re=a.Im=0; }
522 Complex conj (const Complex &a) { return !a; }
523 Complex pow (Complex a, int n) { return a^n; }
524 Complex pow (Complex a, double n) { return a^n; }
525 Complex pow (Complex a, Complex b) { return a^b; }
527 inline Complex sqr (const Complex &a) { return a*a; }
529 Complex sqrt (const Complex &a, int flag)
531 if ((a.Re>=0) && (TMath::Abs(a.Im)<small_epsilon))
532 return flag ? -sqrt(a.Re) : sqrt(a.Re);
533 double R = TMath::Abs(a.Re), I = TMath::Abs(a.Im);
534 double w = (R >= I) ?
535 sqrt (R/2 * ( 1 + sqrt (1 + a.Im*a.Im/a.Re/a.Re))):
536 sqrt (I/2 * (R/I + sqrt (1 + a.Re*a.Re/a.Im/a.Im)));
546 c.Im = (a.Im >= 0) ? w : -w;
548 return ((flag && (c.Re<0)) || (!flag && (c.Re>=0))) ? c : -c;
551 Complex operator ^ (const Complex &a, int n)
555 if (n>0) for (int i=0;i<n;i++) c*=a;
556 else for (int j=0;j>n;j--) c*=a;
561 Complex operator ^ (const Complex &a, double n)
563 return exp(n*log(a));
566 Complex operator ^ (const Complex &a, const Complex &b)
568 return exp(b*log(a));
571 Complex root (const Complex &z, int n, int k)
573 double c=exp(log(abs(z))/n);
574 double t=(Arg(z)+2*PI*k)/n;
575 return Complex (c*cos(t), c*sin(t));
578 Complex exp (const Complex &a)
581 return Complex (t*cos(a.Im), t*sin(a.Im));
584 Complex log (const Complex &a)
587 ComplexError("Error in function log(Complex): argument is 0");
588 return Complex (log(abs(a)), Arg(a));
591 Complex sin (const Complex &a)
593 return Complex (sin(a.Re)*cosh(a.Im), cos(a.Re)*sinh(a.Im));
596 Complex cos (const Complex &a)
598 return Complex (cos(a.Re)*cosh(a.Im), -sin(a.Re)*sinh(a.Im));
601 Complex tan (const Complex &a)
603 return sin(a)/cos(a);
606 Complex cot (const Complex &a)
608 return cos(a)/sin(a);
611 Complex sec (const Complex &a)
616 Complex csc (const Complex &a)
621 Complex sinh (const Complex &a)
623 return Complex (sinh(a.Re)*cos(a.Im), cosh(a.Re)*sin(a.Im));
626 Complex cosh (const Complex &a)
628 return Complex (cosh(a.Re)*cos(a.Im), sinh(a.Re)*sin(a.Im));
631 Complex tanh (const Complex &a)
633 return sinh(a)/cosh(a);
636 Complex coth (const Complex &a)
638 return cosh(a)/sinh(a);
641 Complex sech (const Complex &a)
646 Complex csch (const Complex &a)
650 //////////////////////// Inverce trigonometric functions
652 Complex asin (const Complex &a, int flag)
654 return -ImUnit * log(ImUnit*a + sqrt(1-sqr(a), flag));
657 Complex acos (const Complex &a, int flag)
659 return -ImUnit * log(a + ImUnit*sqrt(1-sqr(a), flag));
662 Complex atan (const Complex &a)
664 return ImUnit/2 * log((ImUnit+a)/(ImUnit-a));
667 Complex acot (const Complex &a)
669 return ImUnit/2 * log((a-ImUnit)/(a+ImUnit));
672 Complex asec (const Complex &a, int flag)
674 return acos(1/a, flag);
677 Complex acsc (const Complex &a, int flag)
679 return asin(1/a, flag);
682 Complex asinh (const Complex &a, int flag)
684 return log(a + sqrt(sqr(a)+1, flag));
687 Complex acosh (const Complex &a, int flag)
689 return log(a + sqrt(sqr(a)-1, flag));
692 Complex atanh (const Complex &a)
694 return log((1+a)/(1-a)) / 2;
697 Complex acoth (const Complex &a)
699 return log((a+1)/(a-1)) / 2;
702 Complex asech (const Complex &a, int flag)
704 return acosh(1/a, flag);
707 Complex acsch (const Complex &a, int flag)
709 return asinh(1/a, flag);
712 Complex Polar (double a, double b)
714 return Complex (a*cos(b), a*sin(b));
717 void Complex::SetPolar (double a, double b)
723 void Complex::SetAbs (double a)
730 void Complex::SetArg (double b)
737 void Solve2 (Complex* z, const Complex &b, const Complex &c)
738 // finding z of equation: z^2 + b*z + c = 0
740 Complex t = sqrt(sqr(b)-4*c);
741 Complex q = ((!b * t).Re >= 0) ? (-(b + t) / 2) : (-(b - t) / 2);
746 Complex Solve2 (const Complex &b, const Complex &c, int RootNumber)
748 if ((RootNumber < 0) || (RootNumber > 1))
749 ComplexError ("Error in Solve2: wrong root number");
752 return z[RootNumber];
755 void Solve3 (Complex* z, const Complex &a2, const Complex &a1, const Complex &a0)
756 // finding z of equation: z^3 + a2*z^2 + a1*z + a0 = 0
758 Complex q, r, t, a, b, zero(0,0);
759 q = (sqr(a2) - 3*a1) / 9;
760 r = (2*(a2^3) - 9*a2*a1 + 27*a0) / 54;
761 t = sqrt((r^2) - (q^3));
762 a = ((!r * t) >=0.0) ? -((r+t)^(1.0/3)) : -((r-t)^(1.0/3));
763 b = ((a == zero) ? zero : (q/a));
764 z[0] = -(a+b)/2 - a2/3 + ImUnit*sqrt(3.0)*(a-b)/2;
765 z[1] = -(a+b)/2 - a2/3 - ImUnit*sqrt(3.0)*(a-b)/2;
769 Complex Solve3 (const Complex &a2, const Complex &a1, const Complex &a0, int RootNumber)
771 if ((RootNumber < 0) || (RootNumber > 2))
772 ComplexError ("Error in Solve3: wrong root number");
775 return z[RootNumber];
778 void Solve4 (Complex* z, const Complex &a3, const Complex &a2, const Complex &a1, const Complex &a0)
779 // finding z of equation: z^4 + a3*z^3 + a2*z^2 + a1*z + a0 = 0
781 Complex u[3], t1, t2, t;
782 Solve3 (u, -a2, (a1*a3 - 4*a0), -((a1^2) + a0*(a3^2) - 4*a0*a2));
784 t1 = sqrt((a3^2)/4 - a2 + t);
785 t2 = sqrt((t^2)/4 - a0);
786 Solve2 ( z, (a3/2 - t1), (t/2 + t2));
787 Solve2 (&(z[2]), (a3/2 + t1), (t/2 - t2));
790 Complex Solve4 (const Complex &a3, const Complex &a2, const Complex &a1, const Complex &a0, int RootNumber)
792 if ((RootNumber < 0) || (RootNumber > 3))
793 ComplexError ("Error in Solve4: wrong root number");
795 Solve4 (z,a3,a2,a1,a0);
796 return z[RootNumber];
800 #endif // __COMPLEX_CXX__ //