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 (const char *msg)
  34 {
  35         cout << endl <<  msg << endl;
  36         exit (1);
  37 }
  38
  39 void ComplexWarning (const 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__ //