• API
• FAQ
• Tools
• Trends
• Archive
daily pastebin goal
6%
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