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