PolynomComplexRoots.h

Go to the documentation of this file.

/* Copyright (C) 2020 Facility for Antiproton and Ion Research in Europe,
   Darmstadt SPDX-License-Identifier: GPL-3.0-only
   Authors: Pierre-Alain Loizeau [committer] */
 
#include <iomanip>
#include <iostream>
#include <math.h>
#include <stdlib.h>
 
#define maxiter 500
 
//
// Extract individual real or complex roots from list of quadratic factors
//
int roots(float* a, int n, float* wr, float* wi)
{
    float sq, b2, c, disc;
    int m, numroots;
 
    m = n;
    numroots = 0;
    while (m > 1) {
        b2 = -0.5 * a[m - 2];
        c = a[m - 1];
        disc = b2 * b2 - c;
        if (disc < 0.0) {   // complex roots
            sq = sqrt(-disc);
            wr[m - 2] = b2;
            wi[m - 2] = sq;
            wr[m - 1] = b2;
            wi[m - 1] = -sq;
            numroots += 2;
        } else {   // real roots
            sq = sqrt(disc);
            wr[m - 2] = fabs(b2) + sq;
            if (b2 < 0.0)
                wr[m - 2] = -wr[m - 2];
            if (wr[m - 2] == 0)
                wr[m - 1] = 0;
            else {
                wr[m - 1] = c / wr[m - 2];
                numroots += 2;
            }
            wi[m - 2] = 0.0;
            wi[m - 1] = 0.0;
        }
        m -= 2;
    }
    if (m == 1) {
        wr[0] = -a[0];
        wi[0] = 0.0;
        numroots++;
    }
    return numroots;
}
//
// Deflate polynomial 'a' by dividing out 'quad'. Return quotient
// polynomial in 'b' and error metric based on remainder in 'err'.
//
void deflate(float* a, int n, float* b, float* quad, float* err)
{
    float r, s;
    int i;
 
    r = quad[1];
    s = quad[0];
 
    b[1] = a[1] - r;
 
    for (i = 2; i <= n; i++) {
        b[i] = a[i] - r * b[i - 1] - s * b[i - 2];
    }
    *err = fabs(b[n]) + fabs(b[n - 1]);
}
//
// Find quadratic factor using Bairstow's method (quadratic Newton method).
// A number of ad hoc safeguards are incorporated to prevent stalls due
// to common difficulties, such as zero slope at iteration point, and
// convergence problems.
//
// Bairstow's method is sensitive to the starting estimate. It is possible
// for convergence to fail or for 'wild' values to trigger an overflow.
//
// It is advisable to institute traps for these problems. (To do!)
//
void find_quad(float* a, int n, float* b, float* quad, float* err, int* iter)
{
    float *c, dn, dr, ds, drn, dsn, eps, r, s;
    int i;
 
    c = new float[n + 1];
    c[0] = 1.0;
    r = quad[1];
    s = quad[0];
    eps = 1e-15;
    *iter = 1;
 
    do {
        if (*iter > maxiter)
            break;
        if (((*iter) % 200) == 0) {
            eps *= 10.0;
        }
        b[1] = a[1] - r;
        c[1] = b[1] - r;
 
        for (i = 2; i <= n; i++) {
            b[i] = a[i] - r * b[i - 1] - s * b[i - 2];
            c[i] = b[i] - r * c[i - 1] - s * c[i - 2];
        }
        dn = c[n - 1] * c[n - 3] - c[n - 2] * c[n - 2];
        drn = b[n] * c[n - 3] - b[n - 1] * c[n - 2];
        dsn = b[n - 1] * c[n - 1] - b[n] * c[n - 2];
 
        if (fabs(dn) < 1e-10) {
            if (dn < 0.0)
                dn = -1e-8;
            else
                dn = 1e-8;
        }
        dr = drn / dn;
        ds = dsn / dn;
        r += dr;
        s += ds;
        (*iter)++;
    } while ((fabs(dr) + fabs(ds)) > eps);
    quad[0] = s;
    quad[1] = r;
    *err = fabs(ds) + fabs(dr);
    delete[] c;
}
//
// Differentiate polynomial 'a' returning result in 'b'.
//
void diff_poly(float* a, int n, float* b)
{
    float coef;
    int i;
 
    coef = (float)n;
    b[0] = 1.0;
    for (i = 1; i < n; i++) {
        b[i] = a[i] * ((float)(n - i)) / coef;
    }
}
//
// Attempt to find a reliable estimate of a quadratic factor using modified
// Bairstow's method with provisions for 'digging out' factors associated
// with multiple roots.
//
// This resursive routine operates on the principal that differentiation of
// a polynomial reduces the order of all multiple roots by one, and has no
// other roots in common with it. If a root of the differentiated polynomial
// is a root of the original polynomial, there must be multiple roots at
// that location. The differentiated polynomial, however, has lower order
// and is easier to solve.
//
// When the original polynomial exhibits convergence problems in the
// neighborhood of some potential root, a best guess is obtained and tried
// on the differentiated polynomial. The new best guess is applied
// recursively on continually differentiated polynomials until failure
// occurs. At this point, the previous polynomial is accepted as that with
// the least number of roots at this location, and its estimate is
// accepted as the root.
//
void recurse(float* a, int n, float* b, int m, float* quad, float* err, int* iter)
{
    float *c, *x, rs[2], tst;
 
    if (fabs(b[m]) < 1e-16)
        m--;   // this bypasses roots at zero
    if (m == 2) {
        quad[0] = b[2];
        quad[1] = b[1];
        *err = 0;
        *iter = 0;
        return;
    }
    c = new float[m + 1];
    x = new float[n + 1];
    c[0] = x[0] = 1.0;
    rs[0] = quad[0];
    rs[1] = quad[1];
    *iter = 0;
    find_quad(b, m, c, rs, err, iter);
    tst = fabs(rs[0] - quad[0]) + fabs(rs[1] - quad[1]);
    if (*err < 1e-12) {
        quad[0] = rs[0];
        quad[1] = rs[1];
    }
    // tst will be 'large' if we converge to wrong root
    if (((*iter > 5) && (tst < 1e-4)) || ((*iter > 20) && (tst < 1e-1))) {
        diff_poly(b, m, c);
        recurse(a, n, c, m - 1, rs, err, iter);
        quad[0] = rs[0];
        quad[1] = rs[1];
    }
    delete[] x;
    delete[] c;
}
//
// Top level routine to manage the determination of all roots of the given
// polynomial 'a', returning the quadratic factors (and possibly one linear
// factor) in 'x'.
//
void get_quads(float* a, int n, float* quad, float* x)
{
    float *b, *z, err, tmp;
    int iter, i, m;
 
    if ((tmp = a[0]) != 1.0) {
        a[0] = 1.0;
        for (i = 1; i <= n; i++) {
            a[i] /= tmp;
        }
    }
    if (n == 2) {
        x[0] = a[1];
        x[1] = a[2];
        return;
    } else if (n == 1) {
        x[0] = a[1];
        return;
    }
    m = n;
    b = new float[n + 1];
    z = new float[n + 1];
    b[0] = 1.0;
    for (i = 0; i <= n; i++) {
        z[i] = a[i];
        x[i] = 0.0;
    }
    do {
        if (n > m) {
            quad[0] = 3.14159e-1;
            quad[1] = 2.78127e-1;
        }
        do {   // This loop tries to assure convergence
            for (i = 0; i < 5; i++) {
                find_quad(z, m, b, quad, &err, &iter);
                if ((err > 1e-7) || (iter > maxiter)) {
                    diff_poly(z, m, b);
                    recurse(z, m, b, m - 1, quad, &err, &iter);
                }
                deflate(z, m, b, quad, &err);
                if (err < 0.001)
                    break;
                quad[0] = 9999.;
                quad[1] = 9999.;
            }
            if (err > 0.01) {
                printf("Error! Convergence failure in quadratic x^2 + r*x + s.\n");
                exit(1);
            }
        } while (err > 0.01);
        x[m - 2] = quad[1];
        x[m - 1] = quad[0];
        m -= 2;
        for (i = 0; i <= m; i++) {
            z[i] = b[i];
        }
    } while (m > 2);
    if (m == 2) {
        x[0] = b[1];
        x[1] = b[2];
    } else
        x[0] = b[1];
    delete[] z;
    delete[] b;
}
 
void polynomComplexRoots(float* wr, float* wi, int n, float* a, int& numr)
{
    float* quad = new float[2];
    float* x = new float[n];
 
    // initialize estimate for 1st root pair
    quad[0] = 2.71828e-1;
    quad[1] = 3.14159e-1;
 
    // get roots
    get_quads(a, n, quad, x);
    numr = roots(x, n, wr, wi);
 
    delete[] quad;
    delete[] x;
}