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