Exercise 6.1: polynomial-4.0.cpp

//polynomial-4.0.cpp

//Program to solve polynomial equations up to 3rd order.
//Adapts polynomial-1.0.cpp to cope with cubics.
//Reads coefficients from keyboard and outputs to screen.
//Uses complex numbers for quadratic and cubic solutions.
//Uses functions to organise code.
//Steven Bamford

#include <cstdlib>
#include <iostream>
#include <cmath>
#include <complex>
#include <iomanip>
#include <vector>
using namespace std;

typedef vector< complex<double> > cdvec;

//function declarations
double solve_linear(const double a, const double b);
cdvec solve_quadratic(const double a, const double b, 
                      const double c);
cdvec solve_cubic(const double a, const double b, 
                  const double c, const double d);

int main()
{
  double c0, c1, c2, c3;  //coefficients of polynomial
                   //c0 + c1*x + c2*x*x + c3*x*x*x = 0
  const double epsilon = 1.0e-300; //a few magnitudes
             //larger than the smallest double number
  
  int width = 10;    //width of quad output fields 
  int precision = 6; //precision of quad output values

  //set output format
  cout.setf(ios::fixed, ios::floatfield); //set format
  cout.precision(precision);           //set precision

  //read in coefficients
  cout << "Enter coefficients c0, c1, c2, and c3:" << endl;
  cin >> c0 >> c1 >> c2 >> c3;  //read coeffs from keyboard
 
  //solve appropriate equation and output solutions
  if (abs(c3) > epsilon)
    {
      cdvec x_sol(3);
      x_sol = solve_cubic(c0, c1, c2, c3);
      cout << "Solutions are: " << endl;
      for (int i = 0; i < 3; i++)
	{
	  cout << "(" << setw(width) << x_sol[i].real() << ","
	              << setw(width) << x_sol[i].imag() << ")" 
	              << endl;
	}
    }
  else if (abs(c2) > epsilon)
    {
      cdvec x_sol(2);
      x_sol = solve_quadratic(c0, c1, c2);
      cout << "Solutions are: " << endl;
      for (int i = 0; i < 2; i++)
	{
	  cout << "(" << setw(width) << x_sol[i].real() << ","
	              << setw(width) << x_sol[i].imag() << ")" 
	              << endl;
	}
    }
  else if (abs(c1) > epsilon)
    {
      double x_sol;
      x_sol = solve_linear(c0, c1);
      cout << "Solution is: " << x_sol << endl;
    }
  else
    {
      cout << "No solution. One of c3, c2, or c1 must be non-zero." 
	   << endl;
    }
  
  system("PAUSE");
  return(0);
}

double solve_linear(const double b, const double a)
{
  return -b/a;
}

cdvec solve_quadratic(const double c, const double b,
		      const double a)
{
  complex<double> des;      //variable for the discriminant
  cdvec x(2);               //two complex solutions
  des = b*b - 4*a*c;        //calc discriminant
  x[0] = (-b + sqrt(des))/(2*a);  //solution 1
  x[1] = (-b - sqrt(des))/(2*a);  //solution 2
  return x;
}

cdvec solve_cubic(const double d, const double c,
		   const double b, const double a)
{
  //normalised coefficients
  double a0 = d/a;
  double a1 = c/a;
  double a2 = b/a;

  //working variables
  complex<double> q = a1 / 3.0 - a2*a2 / 9.0;
  complex<double> r = (a1*a2 - 3.0*a0)/6.0 - a2*a2*a2 / 27.0;
  complex<double> s = sqrt(q*q*q + r*r);
  complex<double> s1 = pow((r + s), (1.0/3.0));
  complex<double> s2 = pow((r - s), (1.0/3.0));

  //solutions
  cdvec x(3);
  x[0] = (s1 + s2) - a2 / 3.0;
  complex<double> xa = -0.5 * (s1 + s2) - (a2 / 3.0);
  complex<double> xb = sqrt(3.0) * (s1 - s2) / 2.0;
  complex<double> i(0, 1);
  x[1] = xa + xb * i;
  x[2] = xa - xb * i;

  return x;
}