Corrections needed on HP
[u/mrichter/AliRoot.git] / HBTAN / volya_complex.h
1 //This file is part of CRAB
2 //http://www.nscl.msu.edu/~pratt/freecodes/crab/home.html
3
4 /*--------------------------------------------------
5 complex.cxx 
6 Version 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
25
26 #include <TMath.h>
27 #include <Riostream.h>
28 //#include <math.h>
29 #include <stdlib.h>
30
31 /////////////////////////// Helper Error Function
32
33 void ComplexError (char *msg)
34 {
35         cout << endl <<  msg << endl;
36         exit (1);
37 }
38
39 void ComplexWarning (char *msg)
40 {
41         cout << endl <<  msg << endl;
42 }
43
44 class Complex
45 {
46 public:
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
170
171 inline Complex Complex::operator + () const
172 {
173         return *this;
174 }
175
176 inline Complex Complex::operator + ()
177 {
178         return *this;
179 }
180
181 inline Complex Complex::operator - () const
182 {
183         return Complex (-Re, -Im);
184 }
185
186 inline Complex Complex::operator - ()
187 {
188         return Complex (-Re, -Im);
189 }
190
191 inline Complex Complex::operator ! () const  //  complex conjugate  //
192 {
193         return Complex (Re, -Im);
194 }
195
196 inline Complex Complex::operator ! ()       //  complex conjugate  //
197 {
198         return Complex (Re, -Im);
199 }
200
201 /////////////////////// +=,-=,*=,/=
202
203 double operator += (double &a, const Complex &b)
204 {
205         if (TMath::Abs(b.Im) < small_epsilon) a+=b.Re;
206         else
207                 ComplexError ("Error in double+=Complex: Complex is not double number");
208         return a;
209 }
210
211 double operator -= (double &a, const Complex &b)
212 {
213         if (TMath::Abs(b.Im) < small_epsilon) a-=b.Re;
214         else
215                 ComplexError ("Error in double -=Complex: Complex is not double number");
216         return a;
217 }
218
219 double operator *= (double &a, const Complex &b)
220 {
221         if (TMath::Abs(b.Im) < small_epsilon) a*=b.Re;
222         else
223                 ComplexError ("Error in double *=Complex: Complex is not double number");
224         return a;
225 }
226
227 double operator /= (double &a, const Complex &b)
228 {
229         if (TMath::Abs(b.Im) < small_epsilon) a/=b.Re;
230         else
231                 ComplexError ("Error in double /=Complex: Complex is not double number");
232         return a;
233 }
234
235 inline Complex Complex::operator += (const Complex &a)
236 {
237         Re+=a.Re;
238         Im+=a.Im;
239         return *this;
240 }
241
242 inline Complex Complex::operator += (double a)
243 {
244         Re+=a;
245         return *this;
246 }
247
248 inline Complex Complex::operator -= (const Complex &a)
249 {
250         Re-=a.Re;
251         Im-=a.Im;
252         return *this;
253 }
254
255 inline Complex Complex::operator -= (double a)
256 {
257         Re-=a;
258         return *this;
259 }
260
261 Complex 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
269 inline Complex Complex::operator *= (double a)
270 {
271         Re*=a;
272         Im*=a;
273         return *this;
274 }
275
276 Complex Complex::operator /= (const Complex &a)
277 {
278     double t1, t2, temp;
279     if (TMath::Abs(a.Re) >= TMath::Abs(a.Im))
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
298 inline Complex Complex::operator /= (double a)
299 {
300         Re/=a;
301         Im/=a;
302         return *this;
303 }
304
305 /////////////////////////// +,-,*,/
306
307 inline Complex operator + (const Complex &a, const Complex &b)
308 {
309         return Complex (a.Re+b.Re, a.Im+b.Im);
310 }
311
312 inline Complex operator + (const Complex &a, double b)
313 {
314         return Complex (a.Re+b, a.Im);
315 }
316
317 inline Complex operator + (double b, const Complex &a)
318 {
319         return Complex (a.Re+b, a.Im);
320 }
321
322 inline Complex operator - (const Complex &a, const Complex &b)
323 {
324         return Complex (a.Re-b.Re, a.Im-b.Im);
325 }
326
327 inline Complex operator - (const Complex &a, double b)
328 {
329         return Complex (a.Re-b, a.Im);
330 }
331
332 inline Complex operator - (double b, const Complex &a)
333 {
334         return Complex (b-a.Re, -a.Im);
335 }
336
337 Complex 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
344 inline Complex operator * (const Complex &a, double b)
345 {
346         return Complex (a.Re*b, a.Im*b);
347 }
348
349 inline Complex operator * (double b, const Complex &a)
350 {
351         return Complex (a.Re*b, a.Im*b);
352 }
353
354 inline Complex operator / (const Complex &a, double b)
355 {
356         return Complex (a.Re/b, a.Im/b);
357 }
358
359 inline Complex operator / (const Complex &a, const Complex &b)
360 {
361     double t1, t2;
362     if (TMath::Abs(b.Re) >= TMath::Abs(b.Im))
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
376 inline Complex operator / (double b, const Complex &a)
377 {
378         return (Complex (b,0)) / a;
379 }
380
381 //////////////////////// <<,>>
382
383 ostream& operator << (ostream &stream, const Complex &a)
384 {
385         stream<<"   "<<a.Re<<"  "<<a.Im<<"  ";
386         return stream;
387 }
388
389 istream& operator >> (istream &stream, Complex &a)
390 {
391         stream>>a.Re>>a.Im;
392         return stream;
393 }
394
395 /////////////////////// ==,!=
396
397 inline int operator == (const Complex &a, double b)
398 {
399         return (a.Re==b) && (TMath::Abs(a.Im)<small_epsilon);
400 }
401
402 inline int operator == (double b, const Complex &a)
403 {
404         return (a.Re==b) && (TMath::Abs(a.Im)<small_epsilon);
405 }
406
407 inline int operator == (const Complex &a, const Complex &b)
408 {
409         return (a.Re==b.Re) && (a.Im==b.Im);
410 }
411
412 inline int operator != (const Complex &a, double b)
413 {
414         return (a.Re!=b) || (TMath::Abs(a.Im)>small_epsilon);
415 }
416
417 inline int operator != (double b, const Complex &a)
418 {
419         return (a.Re!=b) || (TMath::Abs(a.Im)>small_epsilon);
420 }
421
422 inline int operator != (const Complex &a, const Complex &b)
423 {
424         return (a.Re!=b.Re) || (a.Im!=b.Im);
425 }
426
427 /////////////////////// >=,<=
428
429 inline int operator >= (const Complex &a, double b)
430 {
431         return abs(a) >= b;
432 }
433
434 inline int operator >= (double b, const Complex &a)
435 {
436         return b >= abs(a);
437 }
438
439 inline int operator >= (const Complex &a, const Complex &b)
440 {
441         return abs(a) >= abs(b);
442 }
443
444 inline int operator <= (const Complex &a, double b)
445 {
446         return abs(a) <= b;
447 }
448
449 inline int operator <= (double b, const Complex &a)
450 {
451         return b <= abs(a);
452 }
453
454 inline int operator <= (const Complex &a, const Complex &b)
455 {
456         return abs(a) <= abs(b);
457 }
458
459 /////////////////////// >,<
460
461 inline int operator > (const Complex &a, double b)
462 {
463         return abs(a) > b;
464 }
465
466 inline int operator > (double b, const Complex &a)
467 {
468         return b > abs(a);
469 }
470
471 inline int operator > (const Complex &a, const Complex &b)
472 {
473         return abs(a) > abs(b);
474 }
475
476 inline int operator < (const Complex &a, double b)
477 {
478         return abs(a) < b;
479 }
480
481 inline int operator < (double b, const Complex &a)
482 {
483         return b < abs(a);
484 }
485
486 inline int operator < (const Complex &a, const Complex &b)
487 {
488         return abs(a) < abs(b);
489 }
490
491 ///////////////////////// Functions
492 double  real (const Complex &a)
493 {
494         return a.Re;
495 }
496 double  imag (const Complex &a)
497 {
498         return a.Im;
499 }
500 double abs (const Complex &a)
501 {
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);
505         return (R >= I) ?
506            (R * sqrt (1 + a.Im*a.Im/a.Re/a.Re)) :
507            (I * sqrt (1 + a.Re*a.Re/a.Im/a.Im)) ;
508 }
509
510 double  Arg (const Complex &a)
511 {
512         return atan2(a.Im,a.Re);
513 }
514
515 double  phase (const Complex &a)
516 {
517         return atan2(a.Im,a.Re);
518 }
519
520 inline void Zero (Complex &a) { a.Re=a.Im=0; }
521
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; }
526
527 inline Complex sqr (const Complex &a) { return a*a; }
528
529 Complex sqrt (const Complex &a, int flag)
530 {
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)));
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
551 Complex 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
561 Complex operator ^ (const Complex &a, double n)
562 {
563         return exp(n*log(a));
564 }
565
566 Complex operator ^ (const Complex &a, const Complex &b)
567 {
568         return exp(b*log(a));
569 }
570
571 Complex 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
578 Complex exp (const Complex &a)
579 {
580         double t=exp(a.Re);
581         return Complex (t*cos(a.Im), t*sin(a.Im));
582 }
583
584 Complex 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
591 Complex sin (const Complex &a)
592 {
593         return Complex (sin(a.Re)*cosh(a.Im), cos(a.Re)*sinh(a.Im));
594 }
595
596 Complex cos (const Complex &a)
597 {
598         return Complex (cos(a.Re)*cosh(a.Im), -sin(a.Re)*sinh(a.Im));
599 }
600
601 Complex tan (const Complex &a)
602 {
603         return sin(a)/cos(a);
604 }
605
606 Complex cot (const Complex &a)
607 {
608         return cos(a)/sin(a);
609 }
610
611 Complex sec (const Complex &a)
612 {
613         return 1/cos(a);
614 }
615
616 Complex csc (const Complex &a)
617 {
618         return 1/sin(a);
619 }
620
621 Complex sinh (const Complex &a)
622 {
623         return Complex (sinh(a.Re)*cos(a.Im), cosh(a.Re)*sin(a.Im));
624 }
625
626 Complex cosh (const Complex &a)
627 {
628         return Complex (cosh(a.Re)*cos(a.Im), sinh(a.Re)*sin(a.Im));
629 }
630
631 Complex tanh (const Complex &a)
632 {
633         return sinh(a)/cosh(a);
634 }
635
636 Complex coth (const Complex &a)
637 {
638         return cosh(a)/sinh(a);
639 }
640
641 Complex sech (const Complex &a)
642 {
643         return 1/cosh(a);
644 }
645
646 Complex csch (const Complex &a)
647 {
648         return 1/sinh(a);
649 }
650 //////////////////////// Inverce trigonometric functions
651
652 Complex asin (const Complex &a, int flag)
653 {
654         return -ImUnit * log(ImUnit*a + sqrt(1-sqr(a), flag));
655 }
656
657 Complex acos (const Complex &a, int flag)
658 {
659         return -ImUnit * log(a + ImUnit*sqrt(1-sqr(a), flag));
660 }
661
662 Complex atan (const Complex &a)
663 {
664         return ImUnit/2 * log((ImUnit+a)/(ImUnit-a));
665 }
666
667 Complex acot (const Complex &a)
668 {
669         return ImUnit/2 * log((a-ImUnit)/(a+ImUnit));
670 }
671
672 Complex asec (const Complex &a, int flag)
673 {
674         return acos(1/a, flag);
675 }
676
677 Complex acsc (const Complex &a, int flag)
678 {
679         return asin(1/a, flag);
680 }
681
682 Complex asinh (const Complex &a, int flag)
683 {
684         return log(a + sqrt(sqr(a)+1, flag));
685 }
686
687 Complex acosh (const Complex &a, int flag)
688 {
689         return log(a + sqrt(sqr(a)-1, flag));
690 }
691
692 Complex atanh (const Complex &a)
693 {
694         return log((1+a)/(1-a)) / 2;
695 }
696
697 Complex acoth (const Complex &a)
698 {
699          return log((a+1)/(a-1)) / 2;
700 }
701
702 Complex asech (const Complex &a, int flag)
703 {
704         return acosh(1/a, flag);
705 }
706
707 Complex acsch (const Complex &a, int flag)
708 {
709         return asinh(1/a, flag);
710 }
711
712 Complex Polar (double a, double b)
713 {
714         return Complex (a*cos(b), a*sin(b));
715 }
716
717 void Complex::SetPolar (double a, double b)
718 {
719         Re=a*cos(b);
720         Im=a*sin(b);
721 }
722
723 void Complex::SetAbs (double a)
724 {
725         double b=Arg(*this);
726         Re=a*cos(b);
727         Im=a*sin(b);
728 }
729
730 void Complex::SetArg (double b)
731 {
732         double a=abs(*this);
733         Re=a*cos(b);
734         Im=a*sin(b);
735 }
736
737 void 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
746 Complex 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     
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
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));
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;
766   z[2] = a + b - a2/3;
767 }
768
769 Complex 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
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
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
790 Complex 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__ //