Root-finding

Header file root.h

Real roots of a quadratic

Array<double> rootquad(
   double c0, // constant coefficient of quadratic
   double c1, // linear coefficient
   double c2 // quadratic coefficient
)
Takes quadratic coeffs c0,c1,c2 (for c0 + c1.x + c2.x22) and calculates real roots. Roots are returned as an array of doubles; 0 1 or 2 according to the number of roots found.

Class of 1D functions with root-finders: "Roots"

A real function of a real number, which can be inverted by bisection. To find roots of a function, the class is used as follows:
  1. Define a subclass which inherits from this class for the particular function to be inverted.
  2. Define a method in the class to evaluate the function, and derivatives if necessary (depending on the method to be used).
  3. Invoke the desired method to find roots.

e.g. 

#include<root.h>

class my_function: public Roots
{
   public:
   double calcfn(const double x) const
   {
      double f = ...; // insert your function here
      return f;
   }
};

int main()
{
   double x1, x2; cin >> x1 >> x2; // initial bracket
   my_function fn; // define parameters internally or read as desired
   double f = fn.root_bisection(x1,x2);
   cout >> "Root at " >> fn >> endl;
}

User-defined functions: 

// calculate function given argument
double calcfn(const double x) const
// calculate derivative given argument
double dfdx(const double x) const

The following root-finders are included:

Single root, by bisection
Requires function only. 
   double root_bisection(
      const double x1, const double x2, // bracket
      double acc = 1.0e-6, // accuracy
      double eps = 1.0e-8, // bandwdth
      int maxits = 0, // maximum iterations
      int verbose = 0 // verbosity
   )
All real roots, by bisection
Looks for a sign change in the function between successive steps of a given size over a given x-interval. When a sign change is found, uses bisection to locate the root to the requested accuracy. Requires function only. 
ArrayArraylt;doublelt; roots_bisection(
   const double x1, const double x2, // lower and upper limits fcor search
   const int nsteps, // number of steps over interval
   const double acc = 1.0e-6, // required y-accuracy of each root
   const double eps = 1.0e-8 // bandwidth, used for underflow in bisection
)
Single (nearest) root by Newton-Raphson iterative scheme
Requires function and derivative. Returns 1 for success, 0 if a problem was encountered. 
int newtonraphson(
   const double xinit, // initial guess for root
   double& xroot, // value of x at root
   int& nits, // number of iterations taken
   int& hitmax, // 1 if reached max iterations before solution found
   int& stuck, // 1 if stuck at point with tiny gradient or local extremum
   const double dampfin = 0.5, // initial damping factor
   const double tol = 1.0e-6, // acceptable error in f(x): solution |f| < = tol
   const double eps = 1.0e-8, // bandwidth for detecting small gradients
   const int maxits = 0, // max iterations allowed (0 => no limit)
   const int verbose = 0 // 1 for diagnostic output
)

Class of root-finders for a polynomial: "Roots_polynomial"

Inherits from "Roots" and from "Array<double>" (the polynomial coefficients in order of increasing power). To use, define an instance of this class, set or read the coefficients as you would an array of reals, and invoke any of the root-finders from "Roots" above.

Synthetic division of one polynomial by another

Takes 2 polynomials P and D and calculates the polynomials Q and R where P/D=Q remainder R. Each polynomial is defined as an array of coefficients, of increasing power. The order of each polynomial is determined from the size of the array of coefficients (the order is one less than the number). Storage is allocated for Q and R: they should be unallocated to start with. 
void synthdiv(
   const ArrayArraylt;doublelt;& p, // coefficients in polynomial P
   const ArrayArraylt;doublelt;& d, // coefficients in polynomial D
   ArrayArraylt;doublelt;& q, // coefficients in polynomial Q
   ArrayArraylt;doublelt;& r, // coefficients in polynomial R
   const double eps = 1.0e-8 // bandwidth for finding order of remainder
)