SHARE
TWEET

Untitled

a guest Jan 11th, 2017 61 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. 'use strict';
  2.  
  3. /*function show(a,b) {
  4.   var i,nums;
  5.   nums = [];
  6.   for(i=0; i<a.length; i++) {
  7.     nums.push('('+a[i]+'+'+b[i]+'i)');
  8.   }
  9.   return '['+nums.join(', ')+']';
  10. }
  11.  
  12. function showPoly(n,a,b) {
  13.   var i, terms;
  14.   if( b===undefined ) {
  15.     terms = [];
  16.     for(i=0; i<n; i++) {
  17.       terms.push(''+a[i]+' * z^' + (n-i-1));
  18.     }
  19.     return terms.join(' + ');
  20.   } else {
  21.     terms = [];
  22.     for(i=0; i<n; i++) {
  23.       terms.push('( '+a[i]+' + '+b[i]+'*%i ) * z^' + (n-i-1));
  24.     }
  25.     return terms.join(' + ');
  26.   }
  27. }
  28. function showPolyPython(n,a,b) {
  29.   var i, terms;
  30.   if( b===undefined ) {
  31.     terms = [];
  32.     for(i=0; i<n; i++) {
  33.       terms.push(''+a[i]+' * z^' + (n-i-1));
  34.     }
  35.     return terms.join(' + ');
  36.   } else {
  37.     terms = [];
  38.     for(i=0; i<n; i++) {
  39.       terms.push(a[i]+' + '+b[i]+'j');
  40.     }
  41.     return '[' + terms.join(', ') + ']';
  42.   }
  43. }*/
  44.  
  45. function polyev ( nn, sr, si, pr, pi, qr, qi, pvri ) {
  46.   //
  47.   // Evaluates a polynomial P at s by the Horner recurrence, placing the
  48.   // partial sums in Q and the computed value in pvri, which for JS is
  49.   // obviously an array since we can't pass by ref.
  50.   //
  51.   var i,t,pvr,pvi;
  52.  
  53.   pvr = qr[0] = pr[0];
  54.   pvi = qi[0] = pi[0];
  55.  
  56.   for(i=1; i<nn; i++) {
  57.     t = pvr * sr - pvi * si + pr[i];
  58.     pvi = pvr * si + pvi * sr + pi[i];
  59.     pvr = t;
  60.     qr[i] = pvr;
  61.     qi[i] = pvi;
  62.   }
  63.   pvri[0] = pvr;
  64.   pvri[1] = pvi;
  65. }
  66.  
  67. function calct( nn, sr, si, hr, hi, qhr, qhi, pvri, tri ) {
  68.   // Computes t = -P(s)/H(s)
  69.   // Returns
  70.  
  71.   var bool, tmp;
  72.   var n = nn - 1;
  73.   var hvri = new Float64Array(2);
  74.  
  75.   polyev( n, sr, si, hr, hi, qhr, qhi, hvri );
  76.  
  77.   var m1 = Math.sqrt(hvri[0]*hvri[0]+hvri[0]*hvri[0]);
  78.   var m2 = Math.sqrt(hr[n-1]*hr[n-1]+hi[n-1]*hi[n-1]);
  79.   bool = (m1 <= 1e-15 * m2);
  80.  
  81.   if( ! bool ) {
  82.     tmp = hvri[0]*hvri[0] + hvri[1]*hvri[1];
  83.     tri[0] = - ( pvri[0]*hvri[0] + pvri[1]*hvri[1] ) / tmp;
  84.     tri[1] = - ( pvri[1]*hvri[0] - pvri[0]*hvri[1] ) / tmp;
  85.     return bool;
  86.   }
  87.  
  88.   tri[0] = 0;
  89.   tri[1] = 0;
  90.  
  91.   return bool;
  92. }
  93.  
  94. function nexth( nn, bool, qhr, qhi, qpr, qpi, hr, hi, tri ) {
  95.   // Calculates the next shifted H polynomial
  96.   // return true if H(s) is essentially zero
  97.   var t1, t2, j;
  98.   var n = nn - 1;
  99.  
  100.   if( ! bool ) {
  101.     for(j=1; j<n; j++) {
  102.       t1 = qhr[j-1];
  103.       t2 = qhi[j-1];
  104.       hr[j] = tri[0] * t1 - tri[1] * t2 + qpr[j];
  105.       hi[j] = tri[0] * t2 + tri[1] * t1 + qpi[j];
  106.     }
  107.     hr[0] = qpr[0];
  108.     hi[0] = qpi[0];
  109.     return;
  110.   }
  111.  
  112.   // If H(s) is zero, replace h with qh:
  113.   for(j=1; j<n; j++) {
  114.     hr[j] = qhr[j-1];
  115.     hi[j] = qhi[j-1];
  116.   }
  117.   hr[0] = 0;
  118.   hi[0] = 0;
  119. }
  120.  
  121. function errev ( nn, qr, qi, ms, mp ) {
  122.   // Bounds the error in evaluating the polynomial by Horner recurrence
  123.   // qr, qi: the partial sums
  124.   // ms: modulus of the point
  125.   // mp: modulus of the polynomial value
  126.   // are, mre: error bounds on complex addition and multiplication
  127.  
  128.   var e, i;
  129.   var are = 1.1e-16;
  130.   var mre = 3.11e-16;
  131.  
  132.   e = Math.sqrt(qr[0]*qr[0]+qi[0]*qi[0]) * mre / (are + mre);
  133.   for(i=0;i<nn;i++) {
  134.     e = e * ms + Math.sqrt(qr[i]*qr[i]+qi[i]*qi[i]);
  135.   }
  136.   return e * (are + mre) - mp * mre;
  137. }
  138.  
  139. function cauchy (nn, pt, q ) {
  140.   // Cauchy computs a lower bound on the moduli ofthe zeros of a polynomial.
  141.   // pt is the modulus of the coefficients.
  142.   //
  143.   // The lower bound of the modulus of the zeros is given by the roots of the polynomial:
  144.   //
  145.   //   n         n-1
  146.   //  z  + |a | z    + ... + |a   | z - |a |
  147.   //         1                 n-1        n
  148.   //
  149.   var x, dx, df, i, n, xm, f;
  150.  
  151.   // The sign of the last coefficient is reversed:
  152.   pt[nn-1] = -pt[nn-1];
  153.  
  154.   // Compute upper estimate of bound:
  155.   n = nn - 1;
  156.   x = Math.exp((Math.log(-pt[nn-1]) - Math.log(pt[0])) / n);
  157.  
  158.   if( pt[n-1] !== 0 ) {
  159.     xm = - pt[nn-1] / pt[n-1];
  160.     if( xm < x ) {
  161.       x = xm;
  162.     }
  163.   }
  164.  
  165.   // Chop the interval (0,x) until f <= 0
  166.   while( true ) {
  167.     xm = x * 0.1;
  168.     f = pt[0];
  169.     for(i=1; i<nn; i++) {
  170.       f = f * xm + pt[i];
  171.     }
  172.     if( f > 0 ) {
  173.       x = xm;
  174.     } else {
  175.       break;
  176.     }
  177.   }
  178.   dx = x;
  179.  
  180.   // Do newton iteration until x converges to two decimal places
  181.   while( Math.abs(dx/x) > 0.005 ) {
  182.     q[0] = pt[0];
  183.     for(i=1; i<nn; i++) {
  184.       q[i] = q[i-1] * x + pt[i];
  185.     }
  186.     f = q[nn-1];
  187.     df = q[0];
  188.     for(i=1; i<n; i++) {
  189.       df = df * x + q[i];
  190.     }
  191.     dx = f / df;
  192.     x -= dx;
  193.   }
  194.  
  195.   return x;
  196. }
  197.  
  198. function noshft (l1, nn, hr, hi, pr, pi) {
  199.   var i, j, jj, nm1, n, tr, ti;
  200.   var xni, t1, t2, tmp;
  201.  
  202.   //console.log('begin stage 1');
  203.  
  204.   n = nn - 1;
  205.   nm1 = n - 1;
  206.  
  207.   // From Eqn. 9.3 for the 'scaled recurrence', calculate
  208.   //
  209.   //  _ (0)        1
  210.   //  H    (z)  =  - P' (z)
  211.   //               n
  212.   //
  213.   // In my copy of the paper, the 'P' appears to be missing a prime. Since the
  214.   // leading coefficient is constrained to be 1, this makes H a monic polynomial.
  215.   for(i=0; i<n; i++) {
  216.     xni = nn - i - 1;
  217.     hr[i] = xni * pr[i] / n;
  218.     hi[i] = xni * pi[i] / n;
  219.   }
  220.  
  221.   for(jj=0; jj<l1; jj++) {
  222.     //console.log('No shift iteration #',jj+1,'H =',showPolyPython(n,hr,hi));
  223.     //console.log('No shift iteration #',jj+1,'P =',showPolyPython(n,pr,pi));
  224.  
  225.     // Compare the trailing coefficient of H to that of P:
  226.     //console.log( Math.sqrt(hr[nm1]*hr[nm1]+hi[nm1]*hi[nm1]));
  227.     //console.log( Math.sqrt(pr[n]*pr[n]+pi[n]*pi[n]));
  228.     if( (hr[nm1]*hr[nm1]+hi[nm1]*hi[nm1]) > 1e-30 * (pr[n]*pr[n]+pi[n]*pi[n]) ) {
  229.  
  230.       // t = - p[n] / h[nm1]
  231.       tmp = - (hr[nm1]*hr[nm1]+hi[nm1]*hi[nm1]);
  232.       tr = ( pr[n]*hr[nm1]+pi[n]*hi[nm1] ) / tmp;
  233.       ti = ( pi[n]*hr[nm1]-pr[n]*hi[nm1] ) / tmp;
  234.  
  235.       // Calculate the new polynomial using the horner recurrence:
  236.       for(j=nm1; j>0; j--) {
  237.         t1 = hr[j-1];
  238.         t2 = hi[j-1];
  239.         hr[j] = tr * t1 - ti * t2 + pr[j];
  240.         hi[j] = tr * t2 + ti * t1 + pi[j];
  241.       }
  242.       hr[0] = pr[0];
  243.       hi[0] = pi[0];
  244.  
  245.     } else {
  246.       //console.log('edge case');
  247.  
  248.       // If the constant term is essentially zero, shift h coefficients
  249.       for(i=n-1; i>0; i--) {
  250.         //console.log("shifting",i);
  251.         hr[i] = hr[i-1];
  252.         hi[i] = hi[i-1];
  253.       }
  254.       hr[0] = 0;
  255.       hi[0] = 0;
  256.       //console.log('shifted H =',showPolyPython(nm1,hr,hi));
  257.     }
  258.   }
  259. }
  260.  
  261. function vrshft (l3, nn, zri, sri, hr, hi, pr, pi, qpr, qpi, qhr, qhi, shr, shi, pvri, tri) {
  262.  
  263.   var mp, ms, omp, relstp, r1, r2, tp, bool, i, j;
  264.  
  265.   //console.log('begin stage 3');
  266.  
  267.   var conv = false;
  268.   var b = false;
  269.  
  270.   sri[0] = zri[0];
  271.   sri[1] = zri[1];
  272.   var eta = 1.1e-16;
  273.  
  274.  
  275.   // Main loop for stage three
  276.   for(i=0; i<l3; i++) {
  277.     //console.log('variable shift iteration #',i+1, ', H =',showPolyPython(nn-1,hr,hi));
  278.  
  279.     // Evaluate P at s and test for convergence
  280.     polyev( nn, sri[0], sri[1], pr, pi, qpr, qpi, pvri );
  281.     mp = Math.sqrt( pvri[0]*pvri[0] + pvri[1]*pvri[1] );
  282.     ms = Math.sqrt( sri[0]*sri[0] + sri[1]*sri[1] );
  283.     var err = errev(nn, qpr, qpi, ms, mp);
  284.     //console.log('compare mp=',mp,' to err=',err);
  285.     if( mp < 20 * err ) {
  286.       // Polynomial value is smaller in value than a bound on the error in evaluating P,
  287.       // terminate the iteration
  288.       conv = true;
  289.       zri[0] = sri[0];
  290.       zri[1] = sri[1];
  291.       //console.log('converged to',zri[0],'+',zri[1]+'i');
  292.       return conv;
  293.     }
  294.  
  295.     if( i !== 0 ) {
  296.       if( ! ( b || mp < omp || relstp >= 0.5 ) ) {
  297.         // Iteration has stalled. Probably a cluster of zeros. Do 5 fixed
  298.         // shift steps into the cluster to force one zero to dominate.
  299.         tp = relstp;
  300.         b = true;
  301.         if( relstp < eta ) {
  302.           tp = eta;
  303.         }
  304.         r1 = Math.sqrt(tp);
  305.         r2 = sri[0] * (1+r1) - sri[1] * r1;
  306.         sri[1] = sri[0] * r1 + sri[1] * (1+r1);
  307.         sri[0] = r2;
  308.         polyev( nn, sri[0], sri[1], pr, pi, qpr, qpi, pvri );
  309.         for(j=1; j<5; j++) {
  310.           calct( nn, sri[0], sri[1], hr, hi, qhr, qhi, pvri, tri );
  311.           nexth( nn, bool, qhr, qhi, qpr, qpi, hr, hi, tri );
  312.         }
  313.         omp = Infinity;
  314.       }
  315.       if( mp * 0.1 > omp ) {
  316.         return conv;
  317.       }
  318.     }
  319.     omp = mp;
  320.  
  321.     bool = calct( nn, sri[0], sri[1], hr, hi, qhr, qhi, pvri, tri );
  322.     bool = nexth( nn, bool, qhr, qhi, qpr, qpi, hr, hi, tri );
  323.     bool = calct( nn, sri[0], sri[1], hr, hi, qhr, qhi, pvri, tri );
  324.     if( ! bool ) {
  325.       relstp = Math.sqrt(tri[0]*tri[0]+tri[1]*tri[1]) / Math.sqrt(sri[0]*sri[0]+sri[1]*sri[1]);
  326.       sri[0] += tri[0];
  327.       sri[1] += tri[1];
  328.     }
  329.   }
  330.   return conv;
  331. }
  332.  
  333. function fxshft (l2, nn, zri, sri, hr, hi, pr, pi, qpr, qpi, qhr, qhi, shr, shi, pvri, tri) {
  334.   var i, j, n, conv;
  335.   var otr, oti, bool, svsr, svsi;
  336.  
  337.   n = nn - 1;
  338.  
  339.   //console.log('begin stage 2');
  340.  
  341.   //console.log('np.polyval(',showPolyPython(nn,pr,pi)+', '+sr+'+'+si+'j)');
  342.   //console.log('np.polyval(',showPolyPython(n,hr,hi)+', '+sr+'+'+si+'j)');
  343.  
  344.   // Evaluate P at s:
  345.   polyev( nn, sri[0], sri[1], pr, pi, qpr, qpi, pvri );
  346.   var test = true;
  347.   var pasd = false;
  348.  
  349.   // Calculate first t = -P(s)/H(s):
  350.   bool = calct( nn, sri[0], sri[1], hr, hi, qhr, qhi, pvri, tri );
  351.   zri[0] = sri[0] + tri[0];
  352.   zri[1] = sri[1] + tri[1];
  353.  
  354.   // Main loop for one second stage step
  355.   for(j=0; j<l2; j++) {
  356.     //console.log('fixed shift iteration #',j+1,' H =',showPolyPython(n,hr,hi));
  357.     otr = tri[0];
  358.     oti = tri[1];
  359.  
  360.     // Compute next h polynomial and new t:
  361.     nexth( nn, bool, qhr, qhi, qpr, qpi, hr, hi, tri );
  362.     bool = calct( nn, sri[0], sri[1], hr, hi, qhr, qhi, pvri, tri );
  363.     zri[0] = sri[0] + tri[0];
  364.     zri[1] = sri[1] + tri[1];
  365.  
  366.  
  367.     // Test for convergence unless stage 3 has failed once or this is the last H polynomial
  368.     if( ! ( bool || !test || j === l2-1) ) {
  369.       var tmpi = tri[0] - otr;
  370.       var tmpr = tri[1] - oti;
  371.       var m1 = Math.sqrt(tmpr*tmpr + tmpi*tmpi);
  372.       var m2 = Math.sqrt(zri[0]*zri[0] + zri[1]*zri[1]);
  373.       if( m1 < 0.5 * m2 ) {
  374.         if( pasd ) {
  375.           //console.log('weak convergence passed twice');
  376.           // The weak convergence test has been passed twice, start the third stage
  377.           // iteration, after saving the current H polynomial and shift
  378.           for(i=0; i<n; i++) {
  379.             shr[i] = hr[i];
  380.             shi[i] = hi[i];
  381.           }
  382.           svsr = sri[0];
  383.           svsi = sri[1];
  384.  
  385.           conv = vrshft( 10, nn, zri, sri, hr, hi, pr, pi, qpr, qpi, qhr, qhi, shr, shi, pvri, tri );
  386.           if( conv ) {
  387.             return conv;
  388.           }
  389.  
  390.           // The iteration failed to converge. Turn off testing and restore H, S, PV, and T.
  391.           //console.log('iteration failed to converge. Turn off testing & restore H, S, PV, T');
  392.           test = false;
  393.           for(i=0; i<n; i++) {
  394.             hr[i] = shr[i];
  395.             hi[i] = shi[i];
  396.           }
  397.           sri[0] = svsr;
  398.           sri[1] = svsi;
  399.           polyev( nn, sri[0], sri[1], pr, pi, qpr, qpi, pvri );
  400.           bool = calct( nn, sri[0], sri[1], hr, hi, qhr, qhi, pvri, tri );
  401.           continue;
  402.         }
  403.         pasd = true;
  404.       } else {
  405.         pasd = false;
  406.       }
  407.     }
  408.   }
  409.  
  410.   conv = vrshft( 10, nn, zri, sri, hr, hi, pr, pi, qpr, qpi, qhr, qhi, shr, shi, pvri, tri );
  411.  
  412.   return conv;
  413. }
  414.  
  415. var cpoly = function cpoly ( opr, opi ) {
  416.   var i, bound, xxx, conv, cnt1, cnt2, tmp, idnn2;
  417.  
  418.   if( opi === undefined ) {
  419.     opi = new Float64Array( opr.length );
  420.   }
  421.  
  422.   if( opr.length !== opi.length ) {
  423.     throw new TypeError('cpoly: error: real/complex polynomial coefficient input length mismatch');
  424.   }
  425.  
  426.   var degree = opr.length - 1;
  427.  
  428.   // Initialization of constants
  429.   var xx = 0.7071067811865475,
  430.       yy = -xx,
  431.       cosr = -0.06975647374412534,
  432.       sinr = 0.9975640502598242,
  433.       nn = degree + 1;
  434.  
  435.   // Output:
  436.   var zeror = new Float64Array(degree),
  437.       zeroi = new Float64Array(degree);
  438.  
  439.   // Allocate arrays
  440.   var pr = new Float64Array(nn),
  441.       pi = new Float64Array(nn),
  442.       hr = new Float64Array(degree),
  443.       hi = new Float64Array(degree),
  444.       qpr = new Float64Array(nn),
  445.       qpi = new Float64Array(nn),
  446.       qhr = new Float64Array(degree),
  447.       qhi = new Float64Array(degree),
  448.       shr = new Float64Array(nn),
  449.       shi = new Float64Array(nn),
  450.       zri = new Float64Array(2), // for pasing around zeros since js doesn't do by reference;
  451.       pvri = new Float64Array(2),  // for passing around the polynomial value, reason = ditto
  452.       tri = new Float64Array(2),   // for passing around T = -P(S)/H(S)
  453.       sri = new Float64Array(2);   // for passing around current s
  454.  
  455.   //console.log('degree =',degree);
  456.   //console.log('nn =',nn);
  457.   // Remove the zeros at the origin if any
  458.   while( opr[nn] === 0 && opi[nn] === 0 ) {
  459.     idnn2 = degree - nn + 1;
  460.     zeror[idnn2] = 0;
  461.     zeroi[idnn2] = 0;
  462.     nn--;
  463.   }
  464.  
  465.   // Make a copy of the coefficients
  466.   for(i=0; i<nn; i++) {
  467.     pr[i] = opr[i];
  468.     pi[i] = opi[i];
  469.     shr[i] = Math.sqrt(pr[i]*pr[i]+pi[i]*pi[i]);
  470.   }
  471.  
  472.   // Skip scaling the polynomial for unusually large
  473.   // or small coefficients. Caveat emptor.
  474.  
  475.   // Start the algorithm for one zero
  476.   while( nn > 2 ) {
  477.  
  478.     // Calculate bound, a lower bound on the modulus of the zeros:
  479.     for(i=0; i<nn; i++) {
  480.       shr[i] = Math.sqrt(pr[i]*pr[i] + pi[i]*pi[i]);
  481.     }
  482.     bound = cauchy( nn, shr, shi);
  483.  
  484.     // Outer loop to control two major passes with different sequences of shifts
  485.     for(cnt1=0; cnt1<2; cnt1++) {
  486.       //console.log("BEGIN OUTER LOOP");
  487.  
  488.       noshft(5, nn, hr, hi, pr, pi);
  489.  
  490.       // Inner loop to select a shift
  491.       for(cnt2=0; cnt2<9; cnt2++) {
  492.         //console.log("BEGIN INNER LOOP");
  493.  
  494.         // rotate shift angle xx and yy:
  495.         xxx = cosr * xx - sinr * yy;
  496.         yy = sinr * xx + cosr * yy;
  497.         xx = xxx;
  498.  
  499.         // Calculate the new shift:
  500.         sri[0] = bound * xx;
  501.         sri[1] = bound * yy;
  502.  
  503.         conv = fxshft( 10*cnt2, nn, zri, sri, hr, hi, pr, pi, qpr, qpi, qhr, qhi, shr, shi, pvri, tri );
  504.  
  505.         if( conv ) {
  506.           //console.log('MAIN LOOP CONVERGENCE. STORE ZERO');
  507.           idnn2 = degree - nn + 1;
  508.           zeror[idnn2] = zri[0];
  509.           zeroi[idnn2] = zri[1];
  510.           nn--;
  511.           //console.log('nn-- = ',nn);
  512.  
  513.           for(i=0; i<nn; i++) {
  514.             pr[i] = qpr[i];
  515.             pi[i] = qpi[i];
  516.           }
  517.           cnt1 = 3; // exit from the cnt2 loop also
  518.           break;
  519.         }
  520.       }
  521.     }
  522.   }
  523.  
  524.   if( nn <= 2 ) {
  525.     //console.log('END LOOP CONVERGENCE. STORE ZERO');
  526.     // Calculate the final zero and return ( - p[1] / p[0] ):
  527.     tmp = pr[0]*pr[0] + pi[0]*pi[0];
  528.     zeror[degree-1] = - ( pr[1]*pr[0] + pi[1]*pi[0] ) / tmp;
  529.     zeroi[degree-1] = - ( pi[1]*pr[0] - pr[1]*pi[0] ) / tmp;
  530.   }
  531.  
  532.   return [zeror, zeroi];
  533. };
RAW Paste Data
Top