complex.c

来自「开放源码的编译器open watcom 1.6.0版的源代码」· C语言 代码 · 共 434 行 · 第 1/2 页

  if (c2.r == 0.0 && c2.i == 0.0)               // If both num and den zero
    this->divide_by_zero ("operator/");         // Raise exception
  if (c2.r == 0.0)                              // If zero real part
    if (this->r < 0.0 && this->i >= 0.0)        // If negative complex
      this->minus_infinity ("operator/");       // Raise exception
    else
      this->plus_infinity ("operator/");        // Raise exception
  if (c2.i == 0.0)                              // If zero real part
    if (this->i < 0.0 && this->r >= 0.0)        // If negative complex
      this->minus_infinity ("operator/");       // Raise exception
    else
      this->plus_infinity ("operator/");        // Raise exception
  double normalize = (c2.r * c2.r) + (c2.i * c2.i);
  double new_r = (this->r * c2.r) + (this->i * c2.i); // multiply with conjugate
  double new_i = (this->i * c2.r) - (this->r * c2.i);
  this->r = new_r / normalize;
  this->i = new_i / normalize;
  return *this;
}


// operator<< -- Overload the output operator for a reference to a CoolComplex
// Input:        Ostream reference, reference to a CoolComplex object
// Output:       Ostream reference

ostream& operator<< (ostream& os, const CoolComplex& c) {
  os << "(" << c.r << "," << c.i << ")";
  return os;
}


// print --  terse print function for CoolComplex
// Input:    reference to output stream
// Output:   none

void CoolComplex::print(ostream& os) {
  os << "                                       /* CoolComplex %lx */" << (unsigned long) this;
}


// curt -- Cubic root of a double

double curt (double d) {
  if (d == 0.0) 
    return 0.0;
  else if (d < 0) 
    return -pow(-d, 1.0/3.0);
  else
    return pow(d, 1.0/3.0);
}

// roots of linear polynomial: a*x + b = 0.

int CoolComplex::roots_of_linear (const double& a, const double& b, 
                                  CoolComplex& r) {
  if (a == 0) {
    if (b == 0) {
      cout << "Homogenous has infinite number of roots." << endl;
      r = 0;                                    // least norm solution
      return 1;
    } else {
      cout << "Inconsistent equation. No root." << endl;
      return 0;
    }
  } else {
    r = -b / a;                                 // normal case, degree=1
    return 1;                                   
  }
}

// roots of quadratic -- a*x^2 + b*x + c = 0.
//         Return two roots, largest magnitude first.

int CoolComplex::roots_of_quadratic (const double& a, const double& b, 
                                     const double& c, 
                                     CoolComplex& r1, CoolComplex& r2) {
  if (a < 0) {                                  // a is make positive, try again.
    return roots_of_quadratic(-a, -b, -c, r1, r2);
  } else if (a == 0) {                          // special cases
    return roots_of_linear(b, c, r1);
  } else if (c == 0) {
    r2 = 0;
    return roots_of_linear(a, b, r1) + 1;
  } else {                                      // normal case
    double discriminant = (b * b) - (4.0 * a * c);
    double a2 = 2 * a;
    if (discriminant < 0) {
      double real = -b / a2;
      double imag = sqrt(-discriminant) / a2;
      r1 = CoolComplex(real, imag);             // two complex roots
      r2 = CoolComplex(real, -imag);
    } else if (discriminant == 0) {             // one double root
      r1 = r2 = -b / a2;
    } else {                                    // two real roots
      double n;
      if (b < 0) n = sqrt(discriminant) - b;
      else       n = -(sqrt(discriminant) + b);
      r1 = n / a2;                              // root with largest 
      r2 = (2 * c) / n;                         // magnitude first.
    }
    return 2;                                   // number of roots=2
  }
}

// roots of cubic -- a*x^3 + b*x^2 + c*x + d = 0.
//         Return three roots, largest magnitude first.


int CoolComplex::roots_of_cubic (const double& a, const double& b, 
                                 const double& c, const double& d,
                                 CoolComplex& r1, CoolComplex& r2, 
                                 CoolComplex& r3) {
  if (a < 0) {
    return roots_of_cubic(-a, -b, -c, -d, r1, r2, r3);
  } else if (a == 0) {                          // special cases
    return roots_of_quadratic(b, c, d, r1, r2);
  } else if (d == 0) {
    r3 = 0;
    return roots_of_quadratic(a, b, c, r1, r2) + 1;
  } else {                                      // normal case
    CoolComplex rs1, rs2;                       // roots of resolvent
    roots_of_quadratic(1, 
                       ((2*b*b*b) + (9 * a * ((3*a*d) - (b*c)))),
                       pow((b*b) - (3*a*c), 3),
                       rs1, rs2);
    if (rs1.imaginary() == 0) {                 // resolvent roots are real
      double r = curt(rs1.real());              // find cube roots of resolvents
      double s = curt(rs2.real());
      double a3 = 3 * a;
      r1 = (r + s - b) / a3;                    // real root first
      double real = (((r + s) / -2) - b) / a3;
      double imag = fabs(((r - s) * (sqrt(3.0) / 2)) / a3);
      r2 = CoolComplex(real, imag);             // two complex conjugate
      r3 = CoolComplex(real, -imag);            // roots last.
    } else {                                    // resolvent roots are complex
      double rho_3 = curt(rs1.modulus());
      double theta_3 = rs1.argument() / 3.0;
      double rd = 2 * rho_3;
      double cd = ::cos(theta_3) / -2.0;
      double sd = ::sin(theta_3) * sqrt(3.0) / 2.0;
      double a3 = 3 * a;
      if (b < 0) {                              // root with largest magnitude
        r1 = CoolComplex(((-2*rd*cd) - b) / a3, 0); // first
        r2 = CoolComplex((rd * (cd + sd) - b) / a3, 0);
        r3 = CoolComplex((rd * (cd - sd) - b) / a3, 0);
      } else {
        r1 = CoolComplex((rd * (cd - sd) - b) / a3, 0);
        r2 = CoolComplex((rd * (cd + sd) - b) / a3, 0);
        r3 = CoolComplex(((-2*rd*cd) - b) / a3, 0);
      }
    }
    return 3;                                   // number of roots=3
  }
}


// roots of quartic -- a*x^4 + b*x^3 + c*x^2 + d*x + e = 0.
//         Return four roots, largest magnitude first.
//         Decompose quartic into two quadratics for better numerical accuracy
//         than directly solving the 4 roots using Ferrari's formula.

int CoolComplex::roots_of_quartic (const double& a, const double& b, 
                                   const double& c, const double& d, 
                                   const double& e,
                                   CoolComplex& r1, CoolComplex& r2, 
                                   CoolComplex& r3, CoolComplex& r4) {
  if (a < 0) {
    return roots_of_quartic(-a, -b, -c, -d, -e, r1, r2, r3, r4);
  } else if (a == 0) {                          // special cases
    return roots_of_cubic(b, c, d, e, r1, r2, r3);
  } else if (e == 0) {
    r4 = 0;
    return roots_of_cubic(a, b, c, d, r1, r2, r3) + 1;
  } else {                                      // normal case
    double s;                                   // most pos real root of resolvent
    {
      CoolComplex rs1, rs2, rs3;                // roots of resolvent
      roots_of_cubic(1,
                     -c,
                     (b * d) - (4 * a * e),
                     (4 * a * c * e) - ((a * d * d) + (b * b * e)),
                     rs1, rs2, rs3);
      if (rs3.imaginary() != 0)                 // 2 resolvent roots are imaginary
        s = rs1;
      else
      if (rs1.real() > rs3.real())              // most positive/negative real
        s = rs1;                                // root is either rs1 or rs3.
      else
        s = rs3;
    }
    double s1 = sqrt((b * b) - (4 * a * (c - s)));
    double s2 = sqrt((s * s) - (4 * a * e));
    double s1s2 = (b * s) - (2 * a * d);
    double a2 = 2 * a;
    if ((s1 * s2 * s1s2) < 0) {                 // same sign?
      roots_of_quadratic(a2,
                         (b - s1),
                         (s + s2),
                         r1, r2);
      roots_of_quadratic(a2,
                         (b + s1),
                         (s - s2),
                         r3, r4);
    } else {
      roots_of_quadratic(a2,
                         (b - s1),
                         (s - s2),
                         r1, r2);
      roots_of_quadratic(a2,
                         (b + s1),
                         (s + s2),
                         r3, r4);
    }
    return 4;                                   // number of roots=4
  }
}