CommitLineData
dd82cadc 1//This file is part of CRAB
2//http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html
3
4/*--------------------------------------------------
5complex.cxx
6Version 1.0.1
7
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--------------------------------------------------*/
13
14#ifndef __COMPLEX_CXX__
15#define __COMPLEX_CXX__
16
17#ifndef small_epsilon
18#define small_epsilon 1e-12
19#endif
20
21#ifndef PI
22#define PI 3.141592653589793238462643
23#endif
24
f93c53d8 25
26#include <TMath.h>
27#include <Riostream.h>
28//#include <math.h>
3294e9ff 29#include <stdlib.h>
dd82cadc 30
31/////////////////////////// Helper Error Function
32
8e8eae84 33void ComplexError (const char *msg)
dd82cadc 34{
35 cout << endl << msg << endl;
36 exit (1);
37}
38
8e8eae84 39void ComplexWarning (const char *msg)
dd82cadc 40{
41 cout << endl << msg << endl;
42}
43
44class Complex
45{
46public:
47 double Re,Im;
48
49 Complex (double a=0, double b=0): Re(a),Im(b) {}
50
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);
57 void SetAbs (double);
58 void SetArg (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);
73
74 friend ostream& operator << (ostream&, const Complex &);
75 friend istream& operator >> (istream&, Complex &);
76};
77
78 const Complex ImUnit (0,1);
79
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);
126
127 Complex operator ^ (const Complex &, int);
128 Complex operator ^ (const Complex &, double);
129 Complex operator ^ (const Complex &, const Complex &);
130
131 double operator += (double&, const Complex &);
132 double operator -= (double&, const Complex &);
133 double operator *= (double&, const Complex &);
134 double operator /= (double&, const Complex &);
135
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 &);
148
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 &);
167
168////////////////////////////////////////////////
169
dd82cadc 170
171inline Complex Complex::operator + () const
172{
173 return *this;
174}
175
176inline Complex Complex::operator + ()
177{
178 return *this;
179}
180
181inline Complex Complex::operator - () const
182{
183 return Complex (-Re, -Im);
184}
185
186inline Complex Complex::operator - ()
187{
188 return Complex (-Re, -Im);
189}
190
191inline Complex Complex::operator ! () const // complex conjugate //
192{
193 return Complex (Re, -Im);
194}
195
196inline Complex Complex::operator ! () // complex conjugate //
197{
198 return Complex (Re, -Im);
199}
200
201/////////////////////// +=,-=,*=,/=
202
203double operator += (double &a, const Complex &b)
204{
3294e9ff 205 if (TMath::Abs(b.Im) < small_epsilon) a+=b.Re;
dd82cadc 206 else
207 ComplexError ("Error in double+=Complex: Complex is not double number");
208 return a;
209}
210
211double operator -= (double &a, const Complex &b)
212{
3294e9ff 213 if (TMath::Abs(b.Im) < small_epsilon) a-=b.Re;
dd82cadc 214 else
215 ComplexError ("Error in double -=Complex: Complex is not double number");
216 return a;
217}
218
219double operator *= (double &a, const Complex &b)
220{
3294e9ff 221 if (TMath::Abs(b.Im) < small_epsilon) a*=b.Re;
dd82cadc 222 else
223 ComplexError ("Error in double *=Complex: Complex is not double number");
224 return a;
225}
226
227double operator /= (double &a, const Complex &b)
228{
3294e9ff 229 if (TMath::Abs(b.Im) < small_epsilon) a/=b.Re;
dd82cadc 230 else
231 ComplexError ("Error in double /=Complex: Complex is not double number");
232 return a;
233}
234
235inline Complex Complex::operator += (const Complex &a)
236{
237 Re+=a.Re;
238 Im+=a.Im;
239 return *this;
240}
241
242inline Complex Complex::operator += (double a)
243{
244 Re+=a;
245 return *this;
246}
247
248inline Complex Complex::operator -= (const Complex &a)
249{
250 Re-=a.Re;
251 Im-=a.Im;
252 return *this;
253}
254
255inline Complex Complex::operator -= (double a)
256{
257 Re-=a;
258 return *this;
259}
260
261Complex Complex::operator *= (const Complex &a)
262{
263 double t1=Re*a.Re, t2=Im*a.Im;
264 Im = (Re+Im) * (a.Re+a.Im) - t1 - t2;
265 Re = t1 - t2;
266 return *this;
267}
268
269inline Complex Complex::operator *= (double a)
270{
271 Re*=a;
272 Im*=a;
273 return *this;
274}
275
276Complex Complex::operator /= (const Complex &a)
277{
278 double t1, t2, temp;
3294e9ff 279 if (TMath::Abs(a.Re) >= TMath::Abs(a.Im))
dd82cadc 280 {
281 t1 = a.Im / a.Re;
282 t2 = a.Re + a.Im * t1;
283 temp = (Re + Im * t1) / t2;
284 Im = (Im - Re * t1) / t2;
285 Re = temp;
286 }
287 else
288 {
289 t1 = a.Re / a.Im;
290 t2 = a.Re * t1 + a.Im;
291 temp = (Re * t1 + Im) / t2;
292 Im = (Im * t1 - Re) / t2;
293 Re = temp;
294 }
295 return *this;
296}
297
298inline Complex Complex::operator /= (double a)
299{
300 Re/=a;
301 Im/=a;
302 return *this;
303}
304
305/////////////////////////// +,-,*,/
306
307inline Complex operator + (const Complex &a, const Complex &b)
308{
309 return Complex (a.Re+b.Re, a.Im+b.Im);
310}
311
312inline Complex operator + (const Complex &a, double b)
313{
314 return Complex (a.Re+b, a.Im);
315}
316
317inline Complex operator + (double b, const Complex &a)
318{
319 return Complex (a.Re+b, a.Im);
320}
321
322inline Complex operator - (const Complex &a, const Complex &b)
323{
324 return Complex (a.Re-b.Re, a.Im-b.Im);
325}
326
327inline Complex operator - (const Complex &a, double b)
328{
329 return Complex (a.Re-b, a.Im);
330}
331
332inline Complex operator - (double b, const Complex &a)
333{
334 return Complex (b-a.Re, -a.Im);
335}
336
337Complex operator * (const Complex &a, const Complex &b)
338{
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);
342}
343
344inline Complex operator * (const Complex &a, double b)
345{
346 return Complex (a.Re*b, a.Im*b);
347}
348
349inline Complex operator * (double b, const Complex &a)
350{
351 return Complex (a.Re*b, a.Im*b);
352}
353
354inline Complex operator / (const Complex &a, double b)
355{
356 return Complex (a.Re/b, a.Im/b);
357}
358
359inline Complex operator / (const Complex &a, const Complex &b)
360{
361 double t1, t2;
3294e9ff 362 if (TMath::Abs(b.Re) >= TMath::Abs(b.Im))
dd82cadc 363 {
364 t1 = b.Im / b.Re;
365 t2 = b.Re + b.Im * t1;
366 return Complex ((a.Re + a.Im * t1) / t2, (a.Im - a.Re * t1) / t2);
367 }
368 else
369 {
370 t1 = b.Re / b.Im;
371 t2 = b.Re * t1 + b.Im;
372 return Complex ((a.Re * t1 + a.Im) / t2, (a.Im * t1 - a.Re) / t2);
373 }
374}
375
376inline Complex operator / (double b, const Complex &a)
377{
378 return (Complex (b,0)) / a;
379}
380
381//////////////////////// <<,>>
382
383ostream& operator << (ostream &stream, const Complex &a)
384{
385 stream<<" "<<a.Re<<" "<<a.Im<<" ";
386 return stream;
387}
388
389istream& operator >> (istream &stream, Complex &a)
390{
391 stream>>a.Re>>a.Im;
392 return stream;
393}
394
395/////////////////////// ==,!=
396
397inline int operator == (const Complex &a, double b)
398{
3294e9ff 399 return (a.Re==b) && (TMath::Abs(a.Im)<small_epsilon);
dd82cadc 400}
401
402inline int operator == (double b, const Complex &a)
403{
3294e9ff 404 return (a.Re==b) && (TMath::Abs(a.Im)<small_epsilon);
dd82cadc 405}
406
407inline int operator == (const Complex &a, const Complex &b)
408{
409 return (a.Re==b.Re) && (a.Im==b.Im);
410}
411
412inline int operator != (const Complex &a, double b)
413{
3294e9ff 414 return (a.Re!=b) || (TMath::Abs(a.Im)>small_epsilon);
dd82cadc 415}
416
417inline int operator != (double b, const Complex &a)
418{
3294e9ff 419 return (a.Re!=b) || (TMath::Abs(a.Im)>small_epsilon);
dd82cadc 420}
421
422inline int operator != (const Complex &a, const Complex &b)
423{
424 return (a.Re!=b.Re) || (a.Im!=b.Im);
425}
426
427/////////////////////// >=,<=
428
429inline int operator >= (const Complex &a, double b)
430{
431 return abs(a) >= b;
432}
433
434inline int operator >= (double b, const Complex &a)
435{
436 return b >= abs(a);
437}
438
439inline int operator >= (const Complex &a, const Complex &b)
440{
441 return abs(a) >= abs(b);
442}
443
444inline int operator <= (const Complex &a, double b)
445{
446 return abs(a) <= b;
447}
448
449inline int operator <= (double b, const Complex &a)
450{
451 return b <= abs(a);
452}
453
454inline int operator <= (const Complex &a, const Complex &b)
455{
456 return abs(a) <= abs(b);
457}
458
459/////////////////////// >,<
460
461inline int operator > (const Complex &a, double b)
462{
463 return abs(a) > b;
464}
465
466inline int operator > (double b, const Complex &a)
467{
468 return b > abs(a);
469}
470
471inline int operator > (const Complex &a, const Complex &b)
472{
473 return abs(a) > abs(b);
474}
475
476inline int operator < (const Complex &a, double b)
477{
478 return abs(a) < b;
479}
480
481inline int operator < (double b, const Complex &a)
482{
483 return b < abs(a);
484}
485
486inline int operator < (const Complex &a, const Complex &b)
487{
488 return abs(a) < abs(b);
489}
490
491///////////////////////// Functions
492double real (const Complex &a)
493{
494 return a.Re;
495}
496double imag (const Complex &a)
497{
498 return a.Im;
499}
500double abs (const Complex &a)
501{
3294e9ff 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);
dd82cadc 505 return (R >= I) ?
ec094451 506 (R * sqrt (1 + a.Im*a.Im/a.Re/a.Re)) :
507 (I * sqrt (1 + a.Re*a.Re/a.Im/a.Im)) ;
dd82cadc 508}
509
510double Arg (const Complex &a)
511{
512 return atan2(a.Im,a.Re);
513}
514
515double phase (const Complex &a)
516{
517 return atan2(a.Im,a.Re);
518}
519
520inline void Zero (Complex &a) { a.Re=a.Im=0; }
521
522Complex conj (const Complex &a) { return !a; }
523Complex pow (Complex a, int n) { return a^n; }
524Complex pow (Complex a, double n) { return a^n; }
525Complex pow (Complex a, Complex b) { return a^b; }
526
527inline Complex sqr (const Complex &a) { return a*a; }
528
529Complex sqrt (const Complex &a, int flag)
530{
3294e9ff 531 if ((a.Re>=0) && (TMath::Abs(a.Im)<small_epsilon))
dd82cadc 532 return flag ? -sqrt(a.Re) : sqrt(a.Re);
3294e9ff 533 double R = TMath::Abs(a.Re), I = TMath::Abs(a.Im);
dd82cadc 534 double w = (R >= I) ?
ec094451 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)));
dd82cadc 537 Complex c;
538 if (a.Re >= 0)
539 {
540 c.Re = w;
541 c.Im = a.Im / (2*w);
542 }
543 else
544 {
545 c.Re = I / (2*w);
546 c.Im = (a.Im >= 0) ? w : -w;
547 }
548 return ((flag && (c.Re<0)) || (!flag && (c.Re>=0))) ? c : -c;
549}
550
551Complex operator ^ (const Complex &a, int n)
552{
553 Complex c(1,0);
554 if (n==0) return c;
555 if (n>0) for (int i=0;i<n;i++) c*=a;
556 else for (int j=0;j>n;j--) c*=a;
557 if (n>0) return c;
558 else return 1/c;
559}
560
561Complex operator ^ (const Complex &a, double n)
562{
563 return exp(n*log(a));
564}
565
566Complex operator ^ (const Complex &a, const Complex &b)
567{
568 return exp(b*log(a));
569}
570
571Complex root (const Complex &z, int n, int k)
572{
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));
576}
577
578Complex exp (const Complex &a)
579{
580 double t=exp(a.Re);
581 return Complex (t*cos(a.Im), t*sin(a.Im));
582}
583
584Complex log (const Complex &a)
585{
586 if (a==0)
587 ComplexError("Error in function log(Complex): argument is 0");
588 return Complex (log(abs(a)), Arg(a));
589}
590
591Complex sin (const Complex &a)
592{
593 return Complex (sin(a.Re)*cosh(a.Im), cos(a.Re)*sinh(a.Im));
594}
595
596Complex cos (const Complex &a)
597{
598 return Complex (cos(a.Re)*cosh(a.Im), -sin(a.Re)*sinh(a.Im));
599}
600
601Complex tan (const Complex &a)
602{
603 return sin(a)/cos(a);
604}
605
606Complex cot (const Complex &a)
607{
608 return cos(a)/sin(a);
609}
610
611Complex sec (const Complex &a)
612{
613 return 1/cos(a);
614}
615
616Complex csc (const Complex &a)
617{
618 return 1/sin(a);
619}
620
621Complex sinh (const Complex &a)
622{
623 return Complex (sinh(a.Re)*cos(a.Im), cosh(a.Re)*sin(a.Im));
624}
625
626Complex cosh (const Complex &a)
627{
628 return Complex (cosh(a.Re)*cos(a.Im), sinh(a.Re)*sin(a.Im));
629}
630
631Complex tanh (const Complex &a)
632{
633 return sinh(a)/cosh(a);
634}
635
636Complex coth (const Complex &a)
637{
638 return cosh(a)/sinh(a);
639}
640
641Complex sech (const Complex &a)
642{
643 return 1/cosh(a);
644}
645
646Complex csch (const Complex &a)
647{
648 return 1/sinh(a);
649}
650//////////////////////// Inverce trigonometric functions
651
652Complex asin (const Complex &a, int flag)
653{
654 return -ImUnit * log(ImUnit*a + sqrt(1-sqr(a), flag));
655}
656
657Complex acos (const Complex &a, int flag)
658{
659 return -ImUnit * log(a + ImUnit*sqrt(1-sqr(a), flag));
660}
661
662Complex atan (const Complex &a)
663{
664 return ImUnit/2 * log((ImUnit+a)/(ImUnit-a));
665}
666
667Complex acot (const Complex &a)
668{
669 return ImUnit/2 * log((a-ImUnit)/(a+ImUnit));
670}
671
672Complex asec (const Complex &a, int flag)
673{
674 return acos(1/a, flag);
675}
676
677Complex acsc (const Complex &a, int flag)
678{
679 return asin(1/a, flag);
680}
681
682Complex asinh (const Complex &a, int flag)
683{
684 return log(a + sqrt(sqr(a)+1, flag));
685}
686
687Complex acosh (const Complex &a, int flag)
688{
689 return log(a + sqrt(sqr(a)-1, flag));
690}
691
692Complex atanh (const Complex &a)
693{
694 return log((1+a)/(1-a)) / 2;
695}
696
697Complex acoth (const Complex &a)
698{
699 return log((a+1)/(a-1)) / 2;
700}
701
702Complex asech (const Complex &a, int flag)
703{
704 return acosh(1/a, flag);
705}
706
707Complex acsch (const Complex &a, int flag)
708{
709 return asinh(1/a, flag);
710}
711
712Complex Polar (double a, double b)
713{
714 return Complex (a*cos(b), a*sin(b));
715}
716
717void Complex::SetPolar (double a, double b)
718{
719 Re=a*cos(b);
720 Im=a*sin(b);
721}
722
723void Complex::SetAbs (double a)
724{
725 double b=Arg(*this);
726 Re=a*cos(b);
727 Im=a*sin(b);
728}
729
730void Complex::SetArg (double b)
731{
732 double a=abs(*this);
733 Re=a*cos(b);
734 Im=a*sin(b);
735}
736
737void Solve2 (Complex* z, const Complex &b, const Complex &c)
738 // finding z of equation: z^2 + b*z + c = 0
739{
740 Complex t = sqrt(sqr(b)-4*c);
741 Complex q = ((!b * t).Re >= 0) ? (-(b + t) / 2) : (-(b - t) / 2);
742 z[0] = q;
743 z[1] = c/q;
744}
745
746Complex Solve2 (const Complex &b, const Complex &c, int RootNumber)
747{
748 if ((RootNumber < 0) || (RootNumber > 1))
749 ComplexError ("Error in Solve2: wrong root number");
750 Complex z[2];
751 Solve2 (z,b,c);
752 return z[RootNumber];
753}
754
755void 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
757{
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));
f93c53d8 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;
dd82cadc 766 z[2] = a + b - a2/3;
767}
768
769Complex Solve3 (const Complex &a2, const Complex &a1, const Complex &a0, int RootNumber)
770{
771 if ((RootNumber < 0) || (RootNumber > 2))
772 ComplexError ("Error in Solve3: wrong root number");
773 Complex z[3];
774 Solve3 (z,a2,a1,a0);
775 return z[RootNumber];
776}
777
778void 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
780{
781 Complex u[3], t1, t2, t;
782 Solve3 (u, -a2, (a1*a3 - 4*a0), -((a1^2) + a0*(a3^2) - 4*a0*a2));
783 t = *u;
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));
788}
789
790Complex Solve4 (const Complex &a3, const Complex &a2, const Complex &a1, const Complex &a0, int RootNumber)
791{
792 if ((RootNumber < 0) || (RootNumber > 3))
793 ComplexError ("Error in Solve4: wrong root number");
794 Complex z[4];
795 Solve4 (z,a3,a2,a1,a0);
796 return z[RootNumber];
797}
798
799
800#endif // __COMPLEX_CXX__ //