Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- '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];
- };
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement