Skip to content

Instantly share code, notes, and snippets.

@lehni

lehni/poly-roots.js

Created January 11, 2017 20:51
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save lehni/5941121e4e4bdd54540e75420e105a6b to your computer and use it in GitHub Desktop.
Save lehni/5941121e4e4bdd54540e75420e105a6b to your computer and use it in GitHub Desktop.
Download ZIP
Raw
poly-roots.js
'use strict';
/*function show(a,b) {
var i,nums;
nums = [];
for(i=0; i<a.length; i++) {
nums.push('('+a[i]+'+'+b[i]+'i)');
}
return '['+nums.join(', ')+']';
}
function showPoly(n,a,b) {
var i, terms;
if( b===undefined ) {
terms = [];
for(i=0; i<n; i++) {
terms.push(''+a[i]+' * z^' + (n-i-1));
}
return terms.join(' + ');
} else {
terms = [];
for(i=0; i<n; i++) {
terms.push('( '+a[i]+' + '+b[i]+'*%i ) * z^' + (n-i-1));
}
return terms.join(' + ');
}
}
function showPolyPython(n,a,b) {
var i, terms;
if( b===undefined ) {
terms = [];
for(i=0; i<n; i++) {
terms.push(''+a[i]+' * z^' + (n-i-1));
}
return terms.join(' + ');
} else {
terms = [];
for(i=0; i<n; i++) {
terms.push(a[i]+' + '+b[i]+'j');
}
return '[' + terms.join(', ') + ']';
}
}*/
function polyev ( nn, sr, si, pr, pi, qr, qi, pvri ) {
//
// Evaluates a polynomial P at s by the Horner recurrence, placing the
// partial sums in Q and the computed value in pvri, which for JS is
// obviously an array since we can't pass by ref.
//
var i,t,pvr,pvi;
pvr = qr[0] = pr[0];
pvi = qi[0] = pi[0];
for(i=1; i<nn; i++) {
t = pvr * sr - pvi * si + pr[i];
pvi = pvr * si + pvi * sr + pi[i];
pvr = t;
qr[i] = pvr;
qi[i] = pvi;
}
pvri[0] = pvr;
pvri[1] = pvi;
}
function calct( nn, sr, si, hr, hi, qhr, qhi, pvri, tri ) {
// Computes t = -P(s)/H(s)
// Returns
var bool, tmp;
var n = nn - 1;
var hvri = new Float64Array(2);
polyev( n, sr, si, hr, hi, qhr, qhi, hvri );
var m1 = Math.sqrt(hvri[0]*hvri[0]+hvri[0]*hvri[0]);
var m2 = Math.sqrt(hr[n-1]*hr[n-1]+hi[n-1]*hi[n-1]);
bool = (m1 <= 1e-15 * m2);
if( ! bool ) {
tmp = hvri[0]*hvri[0] + hvri[1]*hvri[1];
tri[0] = - ( pvri[0]*hvri[0] + pvri[1]*hvri[1] ) / tmp;
tri[1] = - ( pvri[1]*hvri[0] - pvri[0]*hvri[1] ) / tmp;
return bool;
}
tri[0] = 0;
tri[1] = 0;
return bool;
}
function nexth( nn, bool, qhr, qhi, qpr, qpi, hr, hi, tri ) {
// Calculates the next shifted H polynomial
// return true if H(s) is essentially zero
var t1, t2, j;
var n = nn - 1;
if( ! bool ) {
for(j=1; j<n; j++) {
t1 = qhr[j-1];
t2 = qhi[j-1];
hr[j] = tri[0] * t1 - tri[1] * t2 + qpr[j];
hi[j] = tri[0] * t2 + tri[1] * t1 + qpi[j];
}
hr[0] = qpr[0];
hi[0] = qpi[0];
return;
}
// If H(s) is zero, replace h with qh:
for(j=1; j<n; j++) {
hr[j] = qhr[j-1];
hi[j] = qhi[j-1];
}
hr[0] = 0;
hi[0] = 0;
}
function errev ( nn, qr, qi, ms, mp ) {
// Bounds the error in evaluating the polynomial by Horner recurrence
// qr, qi: the partial sums
// ms: modulus of the point
// mp: modulus of the polynomial value
// are, mre: error bounds on complex addition and multiplication
var e, i;
var are = 1.1e-16;
var mre = 3.11e-16;
e = Math.sqrt(qr[0]*qr[0]+qi[0]*qi[0]) * mre / (are + mre);
for(i=0;i<nn;i++) {
e = e * ms + Math.sqrt(qr[i]*qr[i]+qi[i]*qi[i]);
}
return e * (are + mre) - mp * mre;
}
function cauchy (nn, pt, q ) {
// Cauchy computs a lower bound on the moduli ofthe zeros of a polynomial.
// pt is the modulus of the coefficients.
//
// The lower bound of the modulus of the zeros is given by the roots of the polynomial:
//
// n n-1
// z + |a | z + ... + |a | z - |a |
// 1 n-1 n
//
var x, dx, df, i, n, xm, f;
// The sign of the last coefficient is reversed:
pt[nn-1] = -pt[nn-1];
// Compute upper estimate of bound:
n = nn - 1;
x = Math.exp((Math.log(-pt[nn-1]) - Math.log(pt[0])) / n);
if( pt[n-1] !== 0 ) {
xm = - pt[nn-1] / pt[n-1];
if( xm < x ) {
x = xm;
}
}
// Chop the interval (0,x) until f <= 0
while( true ) {
xm = x * 0.1;
f = pt[0];
for(i=1; i<nn; i++) {
f = f * xm + pt[i];
}
if( f > 0 ) {
x = xm;
} else {
break;
}
}
dx = x;
// Do newton iteration until x converges to two decimal places
while( Math.abs(dx/x) > 0.005 ) {
q[0] = pt[0];
for(i=1; i<nn; i++) {
q[i] = q[i-1] * x + pt[i];
}
f = q[nn-1];
df = q[0];
for(i=1; i<n; i++) {
df = df * x + q[i];
}
dx = f / df;
x -= dx;
}
return x;
}
function noshft (l1, nn, hr, hi, pr, pi) {
var i, j, jj, nm1, n, tr, ti;
var xni, t1, t2, tmp;
//console.log('begin stage 1');
n = nn - 1;
nm1 = n - 1;
// From Eqn. 9.3 for the 'scaled recurrence', calculate
//
// _ (0) 1
// H (z) = - P' (z)
// n
//
// In my copy of the paper, the 'P' appears to be missing a prime. Since the
// leading coefficient is constrained to be 1, this makes H a monic polynomial.
for(i=0; i<n; i++) {
xni = nn - i - 1;
hr[i] = xni * pr[i] / n;
hi[i] = xni * pi[i] / n;
}
for(jj=0; jj<l1; jj++) {
//console.log('No shift iteration #',jj+1,'H =',showPolyPython(n,hr,hi));
//console.log('No shift iteration #',jj+1,'P =',showPolyPython(n,pr,pi));
// Compare the trailing coefficient of H to that of P:
//console.log( Math.sqrt(hr[nm1]*hr[nm1]+hi[nm1]*hi[nm1]));
//console.log( Math.sqrt(pr[n]*pr[n]+pi[n]*pi[n]));
if( (hr[nm1]*hr[nm1]+hi[nm1]*hi[nm1]) > 1e-30 * (pr[n]*pr[n]+pi[n]*pi[n]) ) {
// t = - p[n] / h[nm1]
tmp = - (hr[nm1]*hr[nm1]+hi[nm1]*hi[nm1]);
tr = ( pr[n]*hr[nm1]+pi[n]*hi[nm1] ) / tmp;
ti = ( pi[n]*hr[nm1]-pr[n]*hi[nm1] ) / tmp;
// Calculate the new polynomial using the horner recurrence:
for(j=nm1; j>0; j--) {
t1 = hr[j-1];
t2 = hi[j-1];
hr[j] = tr * t1 - ti * t2 + pr[j];
hi[j] = tr * t2 + ti * t1 + pi[j];
}
hr[0] = pr[0];
hi[0] = pi[0];
} else {
//console.log('edge case');
// If the constant term is essentially zero, shift h coefficients
for(i=n-1; i>0; i--) {
//console.log("shifting",i);
hr[i] = hr[i-1];
hi[i] = hi[i-1];
}
hr[0] = 0;
hi[0] = 0;
//console.log('shifted H =',showPolyPython(nm1,hr,hi));
}
}
}
function vrshft (l3, nn, zri, sri, hr, hi, pr, pi, qpr, qpi, qhr, qhi, shr, shi, pvri, tri) {
var mp, ms, omp, relstp, r1, r2, tp, bool, i, j;
//console.log('begin stage 3');
var conv = false;
var b = false;
sri[0] = zri[0];
sri[1] = zri[1];
var eta = 1.1e-16;
// Main loop for stage three
for(i=0; i<l3; i++) {
//console.log('variable shift iteration #',i+1, ', H =',showPolyPython(nn-1,hr,hi));
// Evaluate P at s and test for convergence
polyev( nn, sri[0], sri[1], pr, pi, qpr, qpi, pvri );
mp = Math.sqrt( pvri[0]*pvri[0] + pvri[1]*pvri[1] );
ms = Math.sqrt( sri[0]*sri[0] + sri[1]*sri[1] );
var err = errev(nn, qpr, qpi, ms, mp);
//console.log('compare mp=',mp,' to err=',err);
if( mp < 20 * err ) {
// Polynomial value is smaller in value than a bound on the error in evaluating P,
// terminate the iteration
conv = true;
zri[0] = sri[0];
zri[1] = sri[1];
//console.log('converged to',zri[0],'+',zri[1]+'i');
return conv;
}
if( i !== 0 ) {
if( ! ( b || mp < omp || relstp >= 0.5 ) ) {
// Iteration has stalled. Probably a cluster of zeros. Do 5 fixed
// shift steps into the cluster to force one zero to dominate.
tp = relstp;
b = true;
if( relstp < eta ) {
tp = eta;
}
r1 = Math.sqrt(tp);
r2 = sri[0] * (1+r1) - sri[1] * r1;
sri[1] = sri[0] * r1 + sri[1] * (1+r1);
sri[0] = r2;
polyev( nn, sri[0], sri[1], pr, pi, qpr, qpi, pvri );
for(j=1; j<5; j++) {
calct( nn, sri[0], sri[1], hr, hi, qhr, qhi, pvri, tri );
nexth( nn, bool, qhr, qhi, qpr, qpi, hr, hi, tri );
}
omp = Infinity;
}
if( mp * 0.1 > omp ) {
return conv;
}
}
omp = mp;
bool = calct( nn, sri[0], sri[1], hr, hi, qhr, qhi, pvri, tri );
bool = nexth( nn, bool, qhr, qhi, qpr, qpi, hr, hi, tri );
bool = calct( nn, sri[0], sri[1], hr, hi, qhr, qhi, pvri, tri );
if( ! bool ) {
relstp = Math.sqrt(tri[0]*tri[0]+tri[1]*tri[1]) / Math.sqrt(sri[0]*sri[0]+sri[1]*sri[1]);
sri[0] += tri[0];
sri[1] += tri[1];
}
}
return conv;
}
function fxshft (l2, nn, zri, sri, hr, hi, pr, pi, qpr, qpi, qhr, qhi, shr, shi, pvri, tri) {
var i, j, n, conv;
var otr, oti, bool, svsr, svsi;
n = nn - 1;
//console.log('begin stage 2');
//console.log('np.polyval(',showPolyPython(nn,pr,pi)+', '+sr+'+'+si+'j)');
//console.log('np.polyval(',showPolyPython(n,hr,hi)+', '+sr+'+'+si+'j)');
// Evaluate P at s:
polyev( nn, sri[0], sri[1], pr, pi, qpr, qpi, pvri );
var test = true;
var pasd = false;
// Calculate first t = -P(s)/H(s):
bool = calct( nn, sri[0], sri[1], hr, hi, qhr, qhi, pvri, tri );
zri[0] = sri[0] + tri[0];
zri[1] = sri[1] + tri[1];
// Main loop for one second stage step
for(j=0; j<l2; j++) {
//console.log('fixed shift iteration #',j+1,' H =',showPolyPython(n,hr,hi));
otr = tri[0];
oti = tri[1];
// Compute next h polynomial and new t:
nexth( nn, bool, qhr, qhi, qpr, qpi, hr, hi, tri );
bool = calct( nn, sri[0], sri[1], hr, hi, qhr, qhi, pvri, tri );
zri[0] = sri[0] + tri[0];
zri[1] = sri[1] + tri[1];
// Test for convergence unless stage 3 has failed once or this is the last H polynomial
if( ! ( bool || !test || j === l2-1) ) {
var tmpi = tri[0] - otr;
var tmpr = tri[1] - oti;
var m1 = Math.sqrt(tmpr*tmpr + tmpi*tmpi);
var m2 = Math.sqrt(zri[0]*zri[0] + zri[1]*zri[1]);
if( m1 < 0.5 * m2 ) {
if( pasd ) {
//console.log('weak convergence passed twice');
// The weak convergence test has been passed twice, start the third stage
// iteration, after saving the current H polynomial and shift
for(i=0; i<n; i++) {
shr[i] = hr[i];
shi[i] = hi[i];
}
svsr = sri[0];
svsi = sri[1];
conv = vrshft( 10, nn, zri, sri, hr, hi, pr, pi, qpr, qpi, qhr, qhi, shr, shi, pvri, tri );
if( conv ) {
return conv;
}
// The iteration failed to converge. Turn off testing and restore H, S, PV, and T.
//console.log('iteration failed to converge. Turn off testing & restore H, S, PV, T');
test = false;
for(i=0; i<n; i++) {
hr[i] = shr[i];
hi[i] = shi[i];
}
sri[0] = svsr;
sri[1] = svsi;
polyev( nn, sri[0], sri[1], pr, pi, qpr, qpi, pvri );
bool = calct( nn, sri[0], sri[1], hr, hi, qhr, qhi, pvri, tri );
continue;
}
pasd = true;
} else {
pasd = false;
}
}
}
conv = vrshft( 10, nn, zri, sri, hr, hi, pr, pi, qpr, qpi, qhr, qhi, shr, shi, pvri, tri );
return conv;
}
var cpoly = function cpoly ( opr, opi ) {
var i, bound, xxx, conv, cnt1, cnt2, tmp, idnn2;
if( opi === undefined ) {
opi = new Float64Array( opr.length );
}
if( opr.length !== opi.length ) {
throw new TypeError('cpoly: error: real/complex polynomial coefficient input length mismatch');
}
var degree = opr.length - 1;
// Initialization of constants
var xx = 0.7071067811865475,
yy = -xx,
cosr = -0.06975647374412534,
sinr = 0.9975640502598242,
nn = degree + 1;
// Output:
var zeror = new Float64Array(degree),
zeroi = new Float64Array(degree);
// Allocate arrays
var pr = new Float64Array(nn),
pi = new Float64Array(nn),
hr = new Float64Array(degree),
hi = new Float64Array(degree),
qpr = new Float64Array(nn),
qpi = new Float64Array(nn),
qhr = new Float64Array(degree),
qhi = new Float64Array(degree),
shr = new Float64Array(nn),
shi = new Float64Array(nn),
zri = new Float64Array(2), // for pasing around zeros since js doesn't do by reference;
pvri = new Float64Array(2), // for passing around the polynomial value, reason = ditto
tri = new Float64Array(2), // for passing around T = -P(S)/H(S)
sri = new Float64Array(2); // for passing around current s
//console.log('degree =',degree);
//console.log('nn =',nn);
// Remove the zeros at the origin if any
while( opr[nn] === 0 && opi[nn] === 0 ) {
idnn2 = degree - nn + 1;
zeror[idnn2] = 0;
zeroi[idnn2] = 0;
nn--;
}
// Make a copy of the coefficients
for(i=0; i<nn; i++) {
pr[i] = opr[i];
pi[i] = opi[i];
shr[i] = Math.sqrt(pr[i]*pr[i]+pi[i]*pi[i]);
}
// Skip scaling the polynomial for unusually large
// or small coefficients. Caveat emptor.
// Start the algorithm for one zero
while( nn > 2 ) {
// Calculate bound, a lower bound on the modulus of the zeros:
for(i=0; i<nn; i++) {
shr[i] = Math.sqrt(pr[i]*pr[i] + pi[i]*pi[i]);
}
bound = cauchy( nn, shr, shi);
// Outer loop to control two major passes with different sequences of shifts
for(cnt1=0; cnt1<2; cnt1++) {
//console.log("BEGIN OUTER LOOP");
noshft(5, nn, hr, hi, pr, pi);
// Inner loop to select a shift
for(cnt2=0; cnt2<9; cnt2++) {
//console.log("BEGIN INNER LOOP");
// rotate shift angle xx and yy:
xxx = cosr * xx - sinr * yy;
yy = sinr * xx + cosr * yy;
xx = xxx;
// Calculate the new shift:
sri[0] = bound * xx;
sri[1] = bound * yy;
conv = fxshft( 10*cnt2, nn, zri, sri, hr, hi, pr, pi, qpr, qpi, qhr, qhi, shr, shi, pvri, tri );
if( conv ) {
//console.log('MAIN LOOP CONVERGENCE. STORE ZERO');
idnn2 = degree - nn + 1;
zeror[idnn2] = zri[0];
zeroi[idnn2] = zri[1];
nn--;
//console.log('nn-- = ',nn);
for(i=0; i<nn; i++) {
pr[i] = qpr[i];
pi[i] = qpi[i];
}
cnt1 = 3; // exit from the cnt2 loop also
break;
}
}
}
}
if( nn <= 2 ) {
//console.log('END LOOP CONVERGENCE. STORE ZERO');
// Calculate the final zero and return ( - p[1] / p[0] ):
tmp = pr[0]*pr[0] + pi[0]*pi[0];
zeror[degree-1] = - ( pr[1]*pr[0] + pi[1]*pi[0] ) / tmp;
zeroi[degree-1] = - ( pi[1]*pr[0] - pr[1]*pi[0] ) / tmp;
}
return [zeror, zeroi];
};
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment