// rpoly_ak1.cpp - Program for calculating the roots of a polynomial of real coefficients.
// Written in Visual C++ 2005 Express Edition
// 14 July 2007
//
// The sub-routines listed below are translations of the FORTRAN routines included in RPOLY.FOR,
// posted off the NETLIB site as TOMS/493:
//
// http://www.netlib.org/toms/493
//
// TOMS/493 is based on the Jenkins-Traub algorithm.
//
// To distinguish the routines posted below from others, an _ak1 suffix has been appended to them.
//
// Following is a list of the major changes made in the course of translating the TOMS/493 routines
// to the C++ versions posted below:
// 1) All global variables have been eliminated.
// 2) The "FAIL" parameter passed into RPOLY.FOR has been eliminated.
// 3) RPOLY.FOR solves polynomials of degree up to 100, but does not explicitly state this limit.
// rpoly_ak1 explicitly states this limit; uses the macro name MAXDEGREE to specify this limit;
// and does a check to ensure that the user input variable Degree is not greater than MAXDEGREE
// (if it is, an error message is output and rpoly_ak1 terminates). If a user wishes to compute
// roots of polynomials of degree greater than MAXDEGREE, using a macro name like MAXDEGREE provides
// the simplest way of offering this capability.
// 4) All "GO TO" statements have been eliminated.
//
// A small main program is included also, to provide an example of how to use rpoly_ak1. In this
// example, data is input from a file to eliminate the need for a user to type data in via
// the console.
#include "stdafx.h"
#include <iostream>
#include <fstream>
#include <cctype>
#include <cmath>
#include <cfloat>
using namespace std;
#define MAXDEGREE 100
#define MDP1 MAXDEGREE+1
void rpoly_ak1(double op[MDP1], int* Degree, double zeror[MAXDEGREE], double zeroi[MAXDEGREE]);
void Fxshfr_ak1(int L2, int* NZ, double sr, double v, double K[MDP1], int N, double p[MDP1], int NN, double qp[MDP1], double u, double* lzi, double* lzr, double* szi, double* szr);
void QuadSD_ak1(int NN, double u, double v, double p[MDP1], double q[MDP1], double* a, double* b);
int calcSC_ak1(int N, double a, double b, double* a1, double* a3, double* a7, double* c, double* d, double* e, double* f, double* g, double* h, double K[MDP1], double u, double v, double qk[MDP1]);
void nextK_ak1(int N, int tFlag, double a, double b, double a1, double* a3, double* a7, double K[MDP1], double qk[MDP1], double qp[MDP1]);
void newest_ak1(int tFlag, double* uu, double* vv, double a, double a1, double a3, double a7, double b, double c, double d, double f, double g, double h, double u, double v, double K[MDP1], int N, double p[MDP1]);
void QuadIT_ak1(int N, int* NZ, double uu, double vv, double* szr, double* szi, double* lzr, double* lzi, double qp[MDP1], int NN, double* a, double* b, double p[MDP1], double qk[MDP1], double* a1, double* a3, double* a7, double* c, double* d, double* e, double* f, double* g, double* h, double K[MDP1]);
void RealIT_ak1(int* iFlag, int* NZ, double* sss, int N, double p[MDP1], int NN, double qp[MDP1], double* szr, double* szi, double K[MDP1], double qk[MDP1]);
void Quad_ak1(double a, double b1, double c, double* sr, double* si, double* lr, double* li);
void rpoly_ak1(double op[MDP1], int* Degree, double zeror[MAXDEGREE], double zeroi[MAXDEGREE]){
int i, j, jj, l, N, NM1, NN, NZ, zerok;
double K[MDP1], p[MDP1], pt[MDP1], qp[MDP1], temp[MDP1];
double bnd, df, dx, factor, ff, moduli_max, moduli_min, sc, x, xm;
double aa, bb, cc, lzi, lzr, sr, szi, szr, t, u, xx, xxx, yy;
const double RADFAC = 3.14159265358979323846/180; // Degrees-to-radians conversion factor = pi/180
const double lb2 = log(2.0); // Dummy variable to avoid re-calculating this value in loop below
const double lo = DBL_MIN/DBL_EPSILON;
const double cosr = cos(94.0*RADFAC); // = -0.069756474
const double sinr = sin(94.0*RADFAC); // = 0.99756405
if ((*Degree) > MAXDEGREE){
cout << "\nThe entered Degree is greater than MAXDEGREE. Exiting rpoly. No further action taken.\n";
*Degree = -1;
return;
} // End ((*Degree) > MAXDEGREE)
//Do a quick check to see if leading coefficient is 0
if (op[0] != 0){
N = *Degree;
xx = sqrt(0.5); // = 0.70710678
yy = -xx;
// Remove zeros at the origin, if any
j = 0;
while (op[N] == 0){
zeror[j] = zeroi[j] = 0.0;
N--;
j++;
} // End while (op[N] == 0)
NN = N + 1;
// Make a copy of the coefficients
for (i = 0; i < NN; i++) p[i] = op[i];
while (N >= 1){ // Main loop
// Start the algorithm for one zero
if (N <= 2){
// Calculate the final zero or pair of zeros
if (N < 2){
zeror[(*Degree) - 1] = -(p[1]/p[0]);
zeroi[(*Degree) - 1] = 0.0;
} // End if (N < 2)
else { // else N == 2
Quad_ak1(p[0], p[1], p[2], &zeror[(*Degree) - 2], &zeroi[(*Degree) - 2], &zeror[(*Degree) - 1], &zeroi[(*Degree) - 1]);
} // End else N == 2
break;
} // End if (N <= 2)
// Find the largest and smallest moduli of the coefficients
moduli_max = 0.0;
moduli_min = DBL_MAX;
for (i = 0; i < NN; i++){
x = fabs(p[i]);
if (x > moduli_max) moduli_max = x;
if ((x != 0) && (x < moduli_min)) moduli_min = x;
} // End for i
// Scale if there are large or very small coefficients
// Computes a scale factor to multiply the coefficients of the polynomial. The scaling
// is done to avoid overflow and to avoid undetected underflow interfering with the
// convergence criterion.
// The factor is a power of the base.
sc = lo/moduli_min;
if (((sc <= 1.0) && (moduli_max >= 10)) || ((sc > 1.0) && (DBL_MAX/sc >= moduli_max))){
sc = ((sc == 0) ? DBL_MIN : sc);
l = (int)(log(sc)/lb2 + 0.5);
factor = pow(2.0, l);
if (factor != 1.0){
for (i = 0; i < NN; i++) p[i] *= factor;
} // End if (factor != 1.0)
} // End if (((sc <= 1.0) && (moduli_max >= 10)) || ((sc > 1.0) && (DBL_MAX/sc >= moduli_max)))
// Compute lower bound on moduli of zeros
for (i = 0; i < NN; i++) pt[i] = fabs(p[i]);
pt[N] = -(pt[N]);
NM1 = N - 1;
// Compute upper estimate of bound
x = exp((log(-pt[N]) - log(pt[0]))/(double)N);
if (pt[NM1] != 0) {
// If Newton step at the origin is better, use it
xm = -pt[N]/pt[NM1];
x = ((xm < x) ? xm : x);
} // End if (pt[NM1] != 0)
// Chop the interval (0, x) until ff <= 0
xm = x;
do {
x = xm;
xm = 0.1*x;
ff = pt[0];
for (i = 1; i < NN; i++) ff = ff *xm + pt[i];
} while (ff > 0); // End do-while loop
dx = x;
// Do Newton iteration until x converges to two decimal places
do {
df = ff = pt[0];
for (i = 1; i < N; i++){
ff = x*ff + pt[i];
df = x*df + ff;
} // End for i
ff = x*ff + pt[N];
dx = ff/df;
x -= dx;
} while (fabs(dx/x) > 0.005); // End do-while loop
bnd = x;
// Compute the derivative as the initial K polynomial and do 5 steps with no shift
for (i = 1; i < N; i++) K[i] = (double)(N - i)*p[i]/((double)N);
K[0] = p[0];
aa = p[N];
bb = p[NM1];
zerok = ((K[NM1] == 0) ? 1 : 0);
for (jj = 0; jj < 5; jj++) {
cc = K[NM1];
if (zerok){
// Use unscaled form of recurrence
for (i = 0; i < NM1; i++){
j = NM1 - i;
K[j] = K[j - 1];
} // End for i
K[0] = 0;
zerok = ((K[NM1] == 0) ? 1 : 0);
} // End if (zerok)
else { // else !zerok
// Used scaled form of recurrence if value of K at 0 is nonzero
t = -aa/cc;
for (i = 0; i < NM1; i++){
j = NM1 - i;
K[j] = t*K[j - 1] + p[j];
} // End for i
K[0] = p[0];
zerok = ((fabs(K[NM1]) <= fabs(bb)*DBL_EPSILON*10.0) ? 1 : 0);
} // End else !zerok
} // End for jj
// Save K for restarts with new shifts
for (i = 0; i < N; i++) temp[i] = K[i];
// Loop to select the quadratic corresponding to each new shift
for (jj = 1; jj <= 20; jj++){
// Quadratic corresponds to a double shift to a non-real point and its
// complex conjugate. The point has modulus BND and amplitude rotated
// by 94 degrees from the previous shift.
xxx = -(sinr*yy) + cosr*xx;
yy = sinr*xx + cosr*yy;
xx = xxx;
sr = bnd*xx;
u = -(2.0*sr);
// Second stage calculation, fixed quadratic
Fxshfr_ak1(20*jj, &NZ, sr, bnd, K, N, p, NN, qp, u, &lzi, &lzr, &szi, &szr);
if (NZ != 0){
// The second stage jumps directly to one of the third stage iterations and
// returns here if successful. Deflate the polynomial, store the zero or
// zeros, and return to the main algorithm.
j = (*Degree) - N;
zeror[j] = szr;
zeroi[j] = szi;
NN = NN - NZ;
N = NN - 1;
for (i = 0; i < NN; i++) p[i] = qp[i];
if (NZ != 1){
zeror[j + 1] = lzr;
zeroi[j + 1] = lzi;
} // End if (NZ != 1)
break;
} // End if (NZ != 0)
else { // Else (NZ == 0)
// If the iteration is unsuccessful, another quadratic is chosen after restoring K
for (i = 0; i < N; i++) K[i] = temp[i];
} // End else (NZ == 0)
} // End for jj
// Return with failure if no convergence with 20 shifts
if (jj > 20) {
cout << "\nFailure. No convergence after 20 shifts. Program terminated.\n";
*Degree -= N;
break;
} // End if (jj > 20)
} // End while (N >= 1)
} // End if op[0] != 0
else { // else op[0] == 0
cout << "\nThe leading coefficient is zero. No further action taken. Program terminated.\n";
*Degree = 0;
} // End else op[0] == 0
return;
} // End rpoly_ak1
void Fxshfr_ak1(int L2, int* NZ, double sr, double v, double K[MDP1], int N, double p[MDP1], int NN, double qp[MDP1], double u, double* lzi, double* lzr, double* szi, double* szr){
// Computes up to L2 fixed shift K-polynomials, testing for convergence in the linear or
// quadratic case. Initiates one of the variable shift iterations and returns with the
// number of zeros found.
// L2 limit of fixed shift steps
// NZ number of zeros found
int fflag, i, iFlag = 1, j, spass, stry, tFlag, vpass, vtry;
double a, a1, a3, a7, b, betas, betav, c, d, e, f, g, h, oss, ots, otv, ovv, s, ss, ts, tss, tv, tvv, ui, vi, vv;
double qk[MDP1], svk[MDP1];
*NZ = 0;
betav = betas = 0.25;
oss = sr;
ovv = v;
//Evaluate polynomial by synthetic division
QuadSD_ak1(NN, u, v, p, qp, &a, &b);
tFlag = calcSC_ak1(N, a, b, &a1, &a3, &a7, &c, &d, &e, &f, &g, &h, K, u, v, qk);
for (j = 0; j < L2; j++){
fflag = 1;
//Calculate next K polynomial and estimate v
nextK_ak1(N, tFlag, a, b, a1, &a3, &a7, K, qk, qp);
tFlag = calcSC_ak1(N, a, b, &a1, &a3, &a7, &c, &d, &e, &f, &g, &h, K, u, v, qk);
newest_ak1(tFlag, &ui, &vi, a, a1, a3, a7, b, c, d, f, g, h, u, v, K, N, p);
vv = vi;
// Estimate s
ss = ((K[N - 1] != 0.0) ? -(p[N]/K[N - 1]) : 0.0);
ts = tv = 1.0;
if ((j != 0) && (tFlag != 3)){
// Compute relative measures of convergence of s and v sequences
tv = ((vv != 0.0) ? fabs((vv - ovv)/vv) : tv);
ts = ((ss != 0.0) ? fabs((ss - oss)/ss) : ts);
// If decreasing, multiply the two most recent convergence measures
tvv = ((tv < otv) ? tv*otv : 1.0);
tss = ((ts < ots) ? ts*ots : 1.0);
// Compare with convergence criteria
vpass = ((tvv < betav) ? 1 : 0);
spass = ((tss < betas) ? 1 : 0);
if ((spass) || (vpass)){
// At least one sequence has passed the convergence test.
// Store variables before iterating
for (i = 0; i < N; i++) svk[i] = K[i];
s = ss;
// Choose iteration according to the fastest converging sequence
stry = vtry = 0;
for ( ; ; ) {
if ((fflag && ((fflag = 0) == 0)) && ((spass) && (!vpass || (tss < tvv)))){
; // Do nothing. Provides a quick "short circuit".
} // End if (fflag)
else { // else !fflag
QuadIT_ak1(N, NZ, ui, vi, szr, szi, lzr, lzi, qp, NN, &a, &b, p, qk, &a1, &a3, &a7, &c, &d, &e, &f, &g, &h, K);
if ((*NZ) > 0) return;
// Quadratic iteration has failed. Flag that it has been tried and decrease the
// convergence criterion
iFlag = vtry = 1;
betav *= 0.25;
// Try linear iteration if it has not been tried and the s sequence is converging
if (stry || (!spass)){
iFlag = 0;
} // End if (stry || (!spass))
else {
for (i = 0; i < N; i++) K[i] = svk[i];
} // End if (stry || !spass)
} // End else fflag
if (iFlag != 0){
RealIT_ak1(&iFlag, NZ, &s, N, p, NN, qp, szr, szi, K, qk);
if ((*NZ) > 0) return;
// Linear iteration has failed. Flag that it has been tried and decrease the
// convergence criterion
stry = 1;
betas *= 0.25;
if (iFlag != 0){
// If linear iteration signals an almost double real zero, attempt quadratic iteration
ui = -(s + s);
vi = s*s;
continue;
} // End if (iFlag != 0)
} // End if (iFlag != 0)
// Restore variables
for (i = 0; i < N; i++) K[i] = svk[i];
// Try quadratic iteration if it has not been tried and the v sequence is converging
if (!vpass || vtry) break; // Break out of infinite for loop
} // End infinite for loop
// Re-compute qp and scalar values to continue the second stage
QuadSD_ak1(NN, u, v, p, qp, &a, &b);
tFlag = calcSC_ak1(N, a, b, &a1, &a3, &a7, &c, &d, &e, &f, &g, &h, K, u, v, qk);
} // End if ((spass) || (vpass))
} // End if ((j != 0) && (tFlag != 3))
ovv = vv;
oss = ss;
otv = tv;
ots = ts;
} // End for j
return;
} // End Fxshfr_ak1
void QuadSD_ak1(int NN, double u, double v, double p[MDP1], double q[MDP1], double* a, double* b){
// Divides p by the quadratic 1, u, v placing the quotient in q and the remainder in a, b
int i;
q[0] = *b = p[0];
q[1] = *a = -((*b)*u) + p[1];
for (i = 2; i < NN; i++){
q[i] = -((*a)*u + (*b)*v) + p[i];
*b = (*a);
*a = q[i];
} // End for i
return;
} // End QuadSD_ak1
int calcSC_ak1(int N, double a, double b, double* a1, double* a3, double* a7, double* c, double* d, double* e, double* f, double* g, double* h, double K[MDP1], double u, double v, double qk[MDP1]){
// This routine calculates scalar quantities used to compute the next K polynomial and
// new estimates of the quadratic coefficients.
// calcSC - integer variable set here indicating how the calculations are normalized
// to avoid overflow.
int dumFlag = 3; // TYPE = 3 indicates the quadratic is almost a factor of K
// Synthetic division of K by the quadratic 1, u, v
QuadSD_ak1(N, u, v, K, qk, c, d);
if (fabs((*c)) <= (100.0*DBL_EPSILON*fabs(K[N - 1]))) {
if (fabs((*d)) <= (100.0*DBL_EPSILON*fabs(K[N - 2]))) return dumFlag;
} // End if (fabs(c) <= (100.0*DBL_EPSILON*fabs(K[N - 1])))
*h = v*b;
if (fabs((*d)) >= fabs((*c))){
dumFlag = 2; // TYPE = 2 indicates that all formulas are divided by d
*e = a/(*d);
*f = (*c)/(*d);
*g = u*b;
*a3 = (*e)*((*g) + a) + (*h)*(b/(*d));
*a1 = -a + (*f)*b;
*a7 = (*h) + ((*f) + u)*a;
} // End if(fabs(d) >= fabs(c))
else {
dumFlag = 1; // TYPE = 1 indicates that all formulas are divided by c;
*e = a/(*c);
*f = (*d)/(*c);
*g = (*e)*u;
*a3 = (*e)*a + ((*g) + (*h)/(*c))*b;
*a1 = -(a*((*d)/(*c))) + b;
*a7 = (*g)*(*d) + (*h)*(*f) + a;
} // End else
return dumFlag;
} // End calcSC_ak1
void nextK_ak1(int N, int tFlag, double a, double b, double a1, double* a3, double* a7, double K[MDP1], double qk[MDP1], double qp[MDP1]){
// Computes the next K polynomials using the scalars computed in calcSC_ak1
int i;
double temp;
if (tFlag == 3){ // Use unscaled form of the recurrence
K[1] = K[0] = 0.0;
for (i = 2; i < N; i++) K[i] = qk[i - 2];
return;
} // End if (tFlag == 3)
temp = ((tFlag == 1) ? b : a);
if (fabs(a1) > (10.0*DBL_EPSILON*fabs(temp))){
// Use scaled form of the recurrence
(*a7) /= a1;
(*a3) /= a1;
K[0] = qp[0];
K[1] = -((*a7)*qp[0]) + qp[1];
for (i = 2; i < N; i++) K[i] = -((*a7)*qp[i - 1]) + (*a3)*qk[i - 2] + qp[i];
} // End if (fabs(a1) > (10.0*DBL_EPSILON*fabs(temp)))
else {
// If a1 is nearly zero, then use a special form of the recurrence
K[0] = 0.0;
K[1] = -(*a7)*qp[0];
for (i = 2; i < N; i++) K[i] = -((*a7)*qp[i - 1]) + (*a3)*qk[i - 2];
} // End else
return;
} // End nextK_ak1
void newest_ak1(int tFlag, double* uu, double* vv, double a, double a1, double a3, double a7, double b, double c, double d, double f, double g, double h, double u, double v, double K[MDP1], int N, double p[MDP1]){
// Compute new estimates of the quadratic coefficients using the scalars computed in calcSC_ak1
double a4, a5, b1, b2, c1, c2, c3, c4, temp;
(*vv) = (*uu) = 0.0; // The quadratic is zeroed
if (tFlag != 3){
if (tFlag != 2){
a4 = a + u*b + h*f;
a5 = c + (u + v*f)*d;
} // End if (tFlag != 2)
else { // else tFlag == 2
a4 = (a + g)*f + h;
a5 = (f + u)*c + v*d;
} // End else tFlag == 2
// Evaluate new quadratic coefficients
b1 = -K[N - 1]/p[N];
b2 = -(K[N - 2] + b1*p[N - 1])/p[N];
c1 = v*b2*a1;
c2 = b1*a7;
c3 = b1*b1*a3;
c4 = -(c2 + c3) + c1;
temp = -c4 + a5 + b1*a4;
if (temp != 0.0) {
*uu= -((u*(c3 + c2) + v*(b1*a1 + b2*a7))/temp) + u;
*vv = v*(1.0 + c4/temp);
} // End if (temp != 0)
} // End if (tFlag != 3)
return;
} // End newest_ak1
void QuadIT_ak1(int N, int* NZ, double uu, double vv, double* szr, double* szi, double* lzr, double* lzi, double qp[MDP1], int NN, double* a, double* b, double p[MDP1], double qk[MDP1], double* a1, double* a3, double* a7, double* c, double* d, double* e, double* f, double* g, double* h, double K[MDP1]){
// Variable-shift K-polynomial iteration for a quadratic factor converges only if the
// zeros are equimodular or nearly so.
int i, j = 0, tFlag, triedFlag = 0;
double ee, mp, omp, relstp, t, u, ui, v, vi, zm;
*NZ = 0; // Number of zeros found
u = uu; // uu and vv are coefficients of the starting quadratic
v = vv;
do {
Quad_ak1(1.0, u, v, szr, szi, lzr, lzi);
// Return if roots of the quadratic are real and not close to multiple or nearly
// equal and of opposite sign.
if (fabs(fabs(*szr) - fabs(*lzr)) > 0.01*fabs(*lzr)) break;
// Evaluate polynomial by quadratic synthetic division
QuadSD_ak1(NN, u, v, p, qp, a, b);
mp = fabs(-((*szr)*(*b)) + (*a)) + fabs((*szi)*(*b));
// Compute a rigorous bound on the rounding error in evaluating p
zm = sqrt(fabs(v));
ee = 2.0*fabs(qp[0]);
t = -((*szr)*(*b));
for (i = 1; i < N; i++) ee = ee*zm + fabs(qp[i]);
ee = ee*zm + fabs((*a) + t);
ee = (9.0*ee + 2.0*fabs(t) - 7.0*(fabs((*a) + t) + zm*fabs((*b))))*DBL_EPSILON;
// Iteration has converged sufficiently if the polynomial value is less than 20 times this bound
if (mp <= 20.0*ee){
*NZ = 2;
break;
} // End if (mp <= 20.0*ee)
j++;
// Stop iteration after 20 steps
if (j > 20) break;
if (j >= 2){
if ((relstp <= 0.01) && (mp >= omp) && (!triedFlag)){
// A cluster appears to be stalling the convergence. Five fixed shift
// steps are taken with a u, v close to the cluster.
relstp = ((relstp < DBL_EPSILON) ? sqrt(DBL_EPSILON) : sqrt(relstp));
u -= u*relstp;
v += v*relstp;
QuadSD_ak1(NN, u, v, p, qp, a, b);
for (i = 0; i < 5; i++){
tFlag = calcSC_ak1(N, *a, *b, a1, a3, a7, c, d, e, f, g, h, K, u, v, qk);
nextK_ak1(N, tFlag, *a, *b, *a1, a3, a7, K, qk, qp);
} // End for i
triedFlag = 1;
j = 0;
} // End if ((relstp <= 0.01) && (mp >= omp) && (!triedFlag))
} // End if (j >= 2)
omp = mp;
// Calculate next K polynomial and new u and v
tFlag = calcSC_ak1(N, *a, *b, a1, a3, a7, c, d, e, f, g, h, K, u, v, qk);
nextK_ak1(N, tFlag, *a, *b, *a1, a3, a7, K, qk, qp);
tFlag = calcSC_ak1(N, *a, *b, a1, a3, a7, c, d, e, f, g, h, K, u, v, qk);
newest_ak1(tFlag, &ui, &vi, *a, *a1, *a3, *a7, *b, *c, *d, *f, *g, *h, u, v, K, N, p);
// If vi is zero, the iteration is not converging
if (vi != 0){
relstp = fabs((-v + vi)/vi);
u = ui;
v = vi;
} // End if (vi != 0)
} while (vi != 0); // End do-while loop
return;
} //End QuadIT_ak1
void RealIT_ak1(int* iFlag, int* NZ, double* sss, int N, double p[MDP1], int NN, double qp[MDP1], double* szr, double* szi, double K[MDP1], double qk[MDP1]){
// Variable-shift H-polynomial iteration for a real zero
// sss - starting iterate
// NZ - number of zeros found
// iFlag - flag to indicate a pair of zeros near real axis
int i, j = 0, nm1 = N - 1;
double ee, kv, mp, ms, omp, pv, s, t;
*iFlag = *NZ = 0;
s = *sss;
for ( ; ; ) {
pv = p[0];
// Evaluate p at s
qp[0] = pv;
for (i = 1; i < NN; i++) qp[i] = pv = pv*s + p[i];
mp = fabs(pv);
// Compute a rigorous bound on the error in evaluating p
ms = fabs(s);
ee = 0.5*fabs(qp[0]);
for (i = 1; i < NN; i++) ee = ee*ms + fabs(qp[i]);
// Iteration has converged sufficiently if the polynomial value is less than
// 20 times this bound
if (mp <= 20.0*DBL_EPSILON*(2.0*ee - mp)){
*NZ = 1;
*szr = s;
*szi = 0.0;
break;
} // End if (mp <= 20.0*DBL_EPSILON*(2.0*ee - mp))
j++;
// Stop iteration after 10 steps
if (j > 10) break;
if (j >= 2){
if ((fabs(t) <= 0.001*fabs(-t + s)) && (mp > omp)){
// A cluster of zeros near the real axis has been encountered;
// Return with iFlag set to initiate a quadratic iteration
*iFlag = 1;
*sss = s;
break;
} // End if ((fabs(t) <= 0.001*fabs(s - t)) && (mp > omp))
} //End if (j >= 2)
// Return if the polynomial value has increased significantly
omp = mp;
// Compute t, the next polynomial and the new iterate
qk[0] = kv = K[0];
for (i = 1; i < N; i++) qk[i] = kv = kv*s + K[i];
if (fabs(kv) > fabs(K[nm1])*10.0*DBL_EPSILON){
// Use the scaled form of the recurrence if the value of K at s is non-zero
t = -(pv/kv);
K[0] = qp[0];
for (i = 1; i < N; i++) K[i] = t*qk[i - 1] + qp[i];
} // End if (fabs(kv) > fabs(K[nm1])*10.0*DBL_EPSILON)
else { // else (fabs(kv) <= fabs(K[nm1])*10.0*DBL_EPSILON)
// Use unscaled form
K[0] = 0.0;
for (i = 1; i < N; i++) K[i] = qk[i - 1];
} // End else (fabs(kv) <= fabs(K[nm1])*10.0*DBL_EPSILON)
kv = K[0];
for (i = 1; i < N; i++) kv = kv*s + K[i];
t = ((fabs(kv) > (fabs(K[nm1])*10.0*DBL_EPSILON)) ? -(pv/kv) : 0.0);
s += t;
} // End infinite for loop
return;
} // End RealIT_ak1
void Quad_ak1(double a, double b1, double c, double* sr, double* si, double* lr, double* li) {
// Calculates the zeros of the quadratic a*Z^2 + b1*Z + c
// The quadratic formula, modified to avoid overflow, is used to find the larger zero if the
// zeros are real and both zeros are complex. The smaller real zero is found directly from
// the product of the zeros c/a.
double b, d, e;
*sr = *si = *lr = *li = 0.0;
if (a == 0) {
*sr = ((b1 != 0) ? -(c/b1) : *sr);
return;
} // End if (a == 0))
if (c == 0){
*lr = -(b1/a);
return;
} // End if (c == 0)
// Compute discriminant avoiding overflow
b = b1/2.0;
if (fabs(b) < fabs(c)){
e = ((c >= 0) ? a : -a);
e = -e + b*(b/fabs(c));
d = sqrt(fabs(e))*sqrt(fabs(c));
} // End if (fabs(b) < fabs(c))
else { // Else (fabs(b) >= fabs(c))
e = -((a/b)*(c/b)) + 1.0;
d = sqrt(fabs(e))*(fabs(b));
} // End else (fabs(b) >= fabs(c))
if (e >= 0) {
// Real zeros
d = ((b >= 0) ? -d : d);
*lr = (-b + d)/a;
*sr = ((*lr != 0) ? (c/(*lr))/a : *sr);
} // End if (e >= 0)
else { // Else (e < 0)
// Complex conjugate zeros
*lr = *sr = -(b/a);
*si = fabs(d/a);
*li = -(*si);
} // End else (e < 0)
return;
} // End Quad_ak1
int main()
{char rflag = 0; //Readiness flag
cout << " rpoly_ak1 (14 July 2007)\n";
cout << "=========================================================================== \n";
cout << "This program calculates the roots of a polynomial of real coefficients:\n";
cout << "\nop[0]*x^N + op[1]*x^(N-1) + op[2]*x^(N-2) + . . . + op[N]*x^0 = 0 \n";
cout << "\n--------------------------------------------------------------------------- \n";
cout << "\nThis program can accept polynomials of degree 100 or less, specified by the\n";
cout << "constant MAXDEGREE. If a higher order polynomial is to be input, redefine\n";
cout << "MAXDEGREE and re-compile the program.\n";
cout << "\n--------------------------------------------------------------------------- \n";
cout << "\nAll relevant data for the polynomial whose roots will be sought should have\n";
cout << "been saved beforehand in a file named rpoly_ak1dat.txt.\n";
cout << "rpoly_ak1dat.txt should be in the same folder as the rpoly_ak1 executable. \n";
cout << "--------------------------------------------------------------------------- \n";
cout << "\nThe first entry of this file must be the degree, N, of the polynomial for\n";
cout << "which the roots are to be calculated.\n";
cout << "Entries for the coefficients of the polynomial should follow, starting with\n";
cout << "the coefficient for the highest power of x and working down to the coefficient\n";
cout << "for the x^0 term.\n";
cout << "\nThe data is assumed to be of type double. Variables used within this program\n";
cout << "are type double.\n";
cout << "\n--------------------------------------------------------------------------- \n";
cout << "\nThe output is written to the file rpoly_ak1out.txt.\n";
cout << "\nNote the returned value of the variable Degree.\n";
cout << "If Degree > 0, it specifies the number of zeros found.\n";
cout << "If Degree = 0, the leading coefficient of the input polynomial was 0.\n";
cout << "If Degree = -1, the input value of Degree was greater than MAXDEGREE.\n";
cout << "\n--------------------------------------------------------------------------- \n";
cout << "\nAdditional information is posted at the following URL:\n";
cout << "http://www.akiti.ca/rpoly_ak1_Intro.html\n";
cout << "--------------------------------------------------------------------------- \n";
cout << "\nIs everything ready (are you ready to continue?)? If yes, Enter y. \n";
cout << "Otherwise Enter any other key. \n";
cin >> rflag;
if (toupper(rflag) == 'Y') {
int Degree; // The degree of the polynomial to be solved
cout << "Appear to be ready. \n";
double op[MDP1], zeroi[MAXDEGREE], zeror[MAXDEGREE]; // Coefficient vectors
int i; // vector index
Degree = 3;
for (i = 0; i < (Degree+1); i++){
if (i<(Degree+1))
{
// I have inserted the equation 2*x^3 + 1*x^2 + 2*x + 1 = 0 below
if (i % 2 == 0 )
{
op[i] = 2;
}
else
{
op[i] =1;
}
}
else
{
op[i]=0;
}
}//End for i
rpoly_ak1(op, &Degree, zeror, zeroi);
cout << "Degree = " << Degree << ".\n";
cout << "\n";
if (Degree <= 0){
cout << "\nReturned from rpoly_ak1 and Degree had a value <= 0.\n";
} // End if (Degree <= 0)
else { // else Degree > 0
cout.precision(DBL_DIG);
cout << "The roots follow:\n";
cout << "\n";
for (i = 0; i < Degree; i++){
cout << zeror[i] << " + " << zeroi[i] << "i" << " \n";
}//End for i
} // End else Degree > 0
cin >> rflag;
} //End if rflag = 'Y'
else cout << "\nNot ready. Try again when ready with information. \n";
cout << "\nEnter any key to continue. \n";
cin >> rflag;
return 0;
} // End main program.