home *** CD-ROM | disk | FTP | other *** search
/ Education Sampler 1992 [NeXTSTEP] / Education_1992_Sampler.iso / SoundAndMusic / cmix / lpc / stabilization / factor.c < prev    next >
C/C++ Source or Header  |  1990-02-03  |  14KB  |  545 lines
  1. /* factor.f -- translated by f2c (version of 26 January 1990  18:57:16).
  2.    You must link the resulting object file with the libraries:
  3.     -lF77 -lI77 -lm -lc   (in that order)
  4. */
  5.  
  6. #include "f2c.h"
  7.  
  8. /* Subroutine */ int factor_(b, k4, rootr, rooti, kinsid, kprint, eps)
  9. doublereal *b;
  10. integer *k4;
  11. doublereal *rootr, *rooti;
  12. integer *kinsid, *kprint;
  13. doublereal *eps;
  14. {
  15.     /* System generated locals */
  16.     integer i_1, i_2, i_3;
  17.     doublereal d_1, d_2;
  18.     doublecomplex z_1, z_2;
  19.  
  20.     /* Builtin functions */
  21.     double atan2();
  22.     /* Subroutine */ int s_stop();
  23.     double pow_dd(), sqrt();
  24.  
  25.     /* Local variables */
  26.     static doublereal amax;
  27.     static integer kerr;
  28.     static doublereal dist, rmin, rmax;
  29.     static integer i, j, k;
  30.     static doublereal r;
  31.     static doublecomplex z;
  32.     static doublereal parti, distr, partr, r2, pi, resmag, resmax;
  33.     extern /* Subroutine */ int dproot_();
  34.     static integer k4m;
  35.     static doublereal coe[37];
  36.     static doublecomplex jay, res;
  37.  
  38. /*
  39. printf("in factornpoles = %d\n",*k4);
  40. for(i=0; i<36; i++) printf(" %g ",b[i]);
  41. */
  42.     /* Parameter adjustments */
  43.     --b;
  44.     --rootr;
  45.     --rooti;
  46.     /* Function Body */
  47. /*       sets up problem, calls dproot, */
  48. /*       and checks residual values at roots */
  49.     jay.r = 0., jay.i = 1.;
  50.     pi = atan2(1., 1.) * 4.;
  51.     i_1 = *k4;
  52.     for (i = 1; i <= i_1; ++i) {
  53. /* L550: */
  54.     coe[i - 1] = b[i];
  55.     }
  56.     k4m = *k4 - 1;
  57.     dproot_(&k4m, coe, &rootr[1], &rooti[1], &kerr, kprint, eps);
  58. /*      write(0,600)kerr */
  59. /* L600: */
  60.     if (kerr > 0) {
  61.     s_stop("", 0L);
  62.     }
  63.     *kinsid = 0;
  64.     resmax = 0.;
  65.     rmax = 0.;
  66.     rmin = 4294967296.;
  67.     dist = 4294967296.;
  68.     amax = 4294967296.;
  69.     d_1 = 1. / *k4;
  70.     r2 = pow_dd(&amax, &d_1);
  71.     i_1 = k4m;
  72.     for (j = 1; j <= i_1; ++j) {
  73.     i_2 = j;
  74.     i_3 = j;
  75.     z_2.r = rooti[i_3] * jay.r, z_2.i = rooti[i_3] * jay.i;
  76.     z_1.r = rootr[i_2] + z_2.r, z_1.i = z_2.i;
  77.     z.r = z_1.r, z.i = z_1.i;
  78. /* Computing 2nd power */
  79.     d_1 = rootr[j];
  80. /* Computing 2nd power */
  81.     d_2 = rooti[j];
  82.     r = sqrt(d_1 * d_1 + d_2 * d_2);
  83. /*         skip residue calculation if root is too big */
  84.     if (r > r2) {
  85.         goto L713;
  86.     }
  87.     i_2 = *k4;
  88.     res.r = b[i_2], res.i = 0.;
  89.     i_2 = *k4;
  90.     for (k = 2; k <= i_2; ++k) {
  91. /* L705: */
  92.         z_2.r = res.r * z.r - res.i * z.i, z_2.i = res.r * z.i + res.i * 
  93.             z.r;
  94.         i_3 = *k4 - k + 1;
  95.         z_1.r = z_2.r + b[i_3], z_1.i = z_2.i;
  96.         res.r = z_1.r, res.i = z_1.i;
  97.     }
  98.     partr = res.r;
  99.     z_2.r = -jay.r, z_2.i = -jay.i;
  100.     z_1.r = z_2.r * res.r - z_2.i * res.i, z_1.i = z_2.r * res.i + z_2.i *
  101.          res.r;
  102.     parti = z_1.r;
  103. /* Computing 2nd power */
  104.     d_1 = partr;
  105. /* Computing 2nd power */
  106.     d_2 = parti;
  107.     resmag = sqrt(d_1 * d_1 + d_2 * d_2);
  108.     if (resmax <= resmag) {
  109.         resmax = resmag;
  110.     }
  111. L713:
  112.     if (rmax < r) {
  113.         rmax = r;
  114.     }
  115.     if (rmin > r) {
  116.         rmin = r;
  117.     }
  118.     if (r < 1.) {
  119.         ++(*kinsid);
  120.     }
  121.     distr = (d_1 = r - 1., abs(d_1));
  122.     if (dist > distr) {
  123.         dist = distr;
  124.     }
  125. /* L701: */
  126.     }
  127. /*     write(0,703)resmax */
  128. /*     write(0,704)rmax,rmin,dist */
  129. /* 703   format('resmax= ',d20.10) */
  130. /* 704   format('rmax= ',d20.10/'rmin= ',d20.10/'dist=',d20.10) */
  131.     return 0;
  132. } /* factor_ */
  133.  
  134. /* Subroutine */ int dproot_(mm, a, rootr, rooti, kerr, kprint, eps)
  135. integer *mm;
  136. doublereal *a, *rootr, *rooti;
  137. integer *kerr, *kprint;
  138. doublereal *eps;
  139. {
  140.     /* System generated locals */
  141.     integer i_1, i_2, i_3, i_4;
  142.     doublereal d_1, d_2;
  143.     doublecomplex z_1, z_2, z_3;
  144.  
  145.     /* Builtin functions */
  146.     double pow_dd(), sqrt(), cos(), sin();
  147.     void z_div();
  148.  
  149.     /* Local variables */
  150.     static integer mdec;
  151.     static doublereal amin, amax, save[37];
  152.     static integer kmax, kpol;
  153.     static doublereal temp, size;
  154.     static integer ktry;
  155.     static doublereal real1, real2;
  156.     static doublecomplex b[37], c[37];
  157.     static integer i, j, k, m;
  158.     static doublecomplex p, w, z;
  159.     static doublereal parti;
  160.     static integer kpolm;
  161.     static doublereal partr;
  162.     static integer kount, newst, nroot;
  163.     static doublereal r1, r2;
  164.     static integer ktrym;
  165.     static doublecomplex bb[37], cc[37];
  166.     static integer mp;
  167.     static doublecomplex pp;
  168.     static doublereal sqteps, rkount, rr1, rr2;
  169.     static doublecomplex jay;
  170.     static integer mmp;
  171.  
  172.     /* Parameter adjustments */
  173.     --a;
  174.     --rootr;
  175.     --rooti;
  176.  
  177.     /* Function Body */
  178. /*        mm=degree of polynomial */
  179. /*        a=coefficient array, lowest to highest degree */
  180. /*        kprint=1 for full printing */
  181. /*        kerr=0 is normal return */
  182.     jay.r = 0., jay.i = 1.;
  183.     mmp = *mm + 1;
  184.     m = *mm;
  185.     mp = mmp;
  186.     i_1 = mp;
  187.     for (i = 1; i <= i_1; ++i) {
  188. /* L700: */
  189.     save[i - 1] = a[i];
  190.     }
  191. /*       kount is number of iterations so far */
  192.     kount = 0;
  193. /*       kmax is maximum total number of iterations allowed */
  194.     kmax = m * 20;
  195. /*       newst is number of re-starts */
  196.     newst = 0;
  197. /*       ktrym is number of attempted iterations before re-starting */
  198.     ktrym = 20;
  199. /*      kpolm is number of attempted iterations before polishing is 
  200. stopped*/
  201.     kpolm = 20;
  202. /*       amax is the largest number we allow */
  203.     amax = 4294967296.;
  204.     amin = 1. / amax;
  205. /*       rr1 and rr2 are radii within which we work for polishing */
  206.     d_1 = 1. / m;
  207.     rr1 = pow_dd(&amin, &d_1);
  208.     d_1 = 1. / m;
  209.     rr2 = pow_dd(&amax, &d_1);
  210. /*     eps is the tolerance for convergence */
  211.     sqteps = sqrt(*eps);
  212. /*        main loop; m is current degree */
  213. L10:
  214.     if (m <= 0) {
  215.     goto L200;
  216.     }
  217. /*        new z, a point on the unit circle */
  218.     rkount = (doublereal) kount;
  219.     d_1 = cos(rkount);
  220.     d_2 = sin(rkount);
  221.     z_2.r = d_2 * jay.r, z_2.i = d_2 * jay.i;
  222.     z_1.r = d_1 + z_2.r, z_1.i = z_2.i;
  223.     z.r = z_1.r, z.i = z_1.i;
  224.     ktry = 0;
  225. /*      r1 and r2 are boundaries of an expanding annulus within which we 
  226. work*/
  227.     d_1 = 1. / m;
  228.     r1 = pow_dd(&amin, &d_1);
  229.     d_1 = 1. / m;
  230.     r2 = pow_dd(&amax, &d_1);
  231. /*        inside loop */
  232. L20:
  233.     partr = z.r;
  234.     z_2.r = -jay.r, z_2.i = -jay.i;
  235.     z_1.r = z_2.r * z.r - z_2.i * z.i, z_1.i = z_2.r * z.i + z_2.i * z.r;
  236.     parti = z_1.r;
  237. /* Computing 2nd power */
  238.     d_1 = partr;
  239. /* Computing 2nd power */
  240.     d_2 = parti;
  241.     size = sqrt(d_1 * d_1 + d_2 * d_2);
  242.     if (size < r1 || size > r2) {
  243.     goto L300;
  244.     }
  245.     if (ktry >= ktrym) {
  246.     goto L300;
  247.     }
  248.     ++ktry;
  249.     if (kount >= kmax) {
  250.     goto L400;
  251.     }
  252.     ++kount;
  253. /*        get value of polynomial at z, synthetic division */
  254.     i_1 = mp - 1;
  255.     i_2 = mp;
  256.     b[i_1].r = a[i_2], b[i_1].i = 0.;
  257.     i_1 = m;
  258.     for (j = 1; j <= i_1; ++j) {
  259.     k = m - j + 1;
  260. /* L30: */
  261.     i_2 = k - 1;
  262.     i_3 = k;
  263.     z_2.r = z.r * b[i_3].r - z.i * b[i_3].i, z_2.i = z.r * b[i_3].i + z.i 
  264.         * b[i_3].r;
  265.     i_4 = k;
  266.     z_1.r = z_2.r + a[i_4], z_1.i = z_2.i;
  267.     b[i_2].r = z_1.r, b[i_2].i = z_1.i;
  268.     }
  269.     p.r = b[0].r, p.i = b[0].i;
  270.     partr = p.r;
  271.     z_2.r = -jay.r, z_2.i = -jay.i;
  272.     z_1.r = z_2.r * p.r - z_2.i * p.i, z_1.i = z_2.r * p.i + z_2.i * p.r;
  273.     parti = z_1.r;
  274. /* Computing 2nd power */
  275.     d_1 = partr;
  276. /* Computing 2nd power */
  277.     d_2 = parti;
  278.     if (sqrt(d_1 * d_1 + d_2 * d_2) > amax) {
  279.     goto L300;
  280.     }
  281. /*        get value of derivative at z, synthetic division */
  282.     i_2 = mp - 1;
  283.     i_3 = mp - 1;
  284.     c[i_2].r = b[i_3].r, c[i_2].i = b[i_3].i;
  285.     mdec = m - 1;
  286.     i_2 = mdec;
  287.     for (j = 1; j <= i_2; ++j) {
  288.     k = m - j + 1;
  289. /* L60: */
  290.     i_3 = k - 1;
  291.     i_4 = k;
  292.     z_2.r = z.r * c[i_4].r - z.i * c[i_4].i, z_2.i = z.r * c[i_4].i + z.i 
  293.         * c[i_4].r;
  294.     i_1 = k - 1;
  295.     z_1.r = z_2.r + b[i_1].r, z_1.i = z_2.i + b[i_1].i;
  296.     c[i_3].r = z_1.r, c[i_3].i = z_1.i;
  297.     }
  298.     pp.r = c[1].r, pp.i = c[1].i;
  299.     partr = pp.r;
  300.     z_2.r = -jay.r, z_2.i = -jay.i;
  301.     z_1.r = z_2.r * pp.r - z_2.i * pp.i, z_1.i = z_2.r * pp.i + z_2.i * pp.r;
  302.     parti = z_1.r;
  303. /* Computing 2nd power */
  304.     d_1 = partr;
  305. /* Computing 2nd power */
  306.     d_2 = parti;
  307.     if (sqrt(d_1 * d_1 + d_2 * d_2) < amin) {
  308.     goto L300;
  309.     }
  310. /*        test for convergence */
  311.     partr = p.r;
  312.     z_2.r = -jay.r, z_2.i = -jay.i;
  313.     z_1.r = z_2.r * p.r - z_2.i * p.i, z_1.i = z_2.r * p.i + z_2.i * p.r;
  314.     parti = z_1.r;
  315. /* Computing 2nd power */
  316.     d_1 = partr;
  317. /* Computing 2nd power */
  318.     d_2 = parti;
  319.     size = sqrt(d_1 * d_1 + d_2 * d_2);
  320.     if (size > *eps) {
  321.     goto L775;
  322.     }
  323.     nroot = *mm - m + 1;
  324.     goto L500;
  325. L775:
  326.     z_div(&z_2, &p, &pp);
  327.     z_1.r = z.r - z_2.r, z_1.i = z.i - z_2.i;
  328.     z.r = z_1.r, z.i = z_1.i;
  329.     goto L20;
  330. /*        end of main loop */
  331. /*        normal return */
  332. L200:
  333.     *kerr = 0;
  334.     goto L600;
  335. /*        new start */
  336. L300:
  337.     rkount = (doublereal) kount;
  338.     d_1 = cos(rkount);
  339.     d_2 = sin(rkount);
  340.     z_2.r = d_2 * jay.r, z_2.i = d_2 * jay.i;
  341.     z_1.r = d_1 + z_2.r, z_1.i = z_2.i;
  342.     z.r = z_1.r, z.i = z_1.i;
  343.     ktry = 0;
  344.     ++newst;
  345.     goto L20;
  346. /*        too many iterations */
  347. L400:
  348.     *kerr = 400;
  349.     goto L600;
  350. /*        root z located */
  351. /*        polish z to get w */
  352. L500:
  353.     w.r = z.r, w.i = z.i;
  354.     kpol = 0;
  355. L510:
  356.     partr = w.r;
  357.     z_2.r = -jay.r, z_2.i = -jay.i;
  358.     z_1.r = z_2.r * w.r - z_2.i * w.i, z_1.i = z_2.r * w.i + z_2.i * w.r;
  359.     parti = z_1.r;
  360. /* Computing 2nd power */
  361.     d_1 = partr;
  362. /* Computing 2nd power */
  363.     d_2 = parti;
  364.     size = sqrt(d_1 * d_1 + d_2 * d_2);
  365. /*       give up polishing if w is outside annulus */
  366.     if (size < rr1 || size > rr2) {
  367.     goto L501;
  368.     }
  369. /*       give up polishing if kpol>=kpolm */
  370.     if (kpol >= kpolm) {
  371.     goto L501;
  372.     }
  373.     ++kpol;
  374.     if (kount >= kmax) {
  375.     goto L400;
  376.     }
  377.     ++kount;
  378.     i_3 = mmp - 1;
  379.     i_4 = mmp - 1;
  380.     bb[i_3].r = save[i_4], bb[i_3].i = 0.;
  381.     i_3 = *mm;
  382.     for (j = 1; j <= i_3; ++j) {
  383.     k = *mm - j + 1;
  384. /* L530: */
  385.     i_4 = k - 1;
  386.     i_1 = k;
  387.     z_2.r = w.r * bb[i_1].r - w.i * bb[i_1].i, z_2.i = w.r * bb[i_1].i + 
  388.         w.i * bb[i_1].r;
  389.     i_2 = k - 1;
  390.     z_1.r = z_2.r + save[i_2], z_1.i = z_2.i;
  391.     bb[i_4].r = z_1.r, bb[i_4].i = z_1.i;
  392.     }
  393.     p.r = bb[0].r, p.i = bb[0].i;
  394.     partr = p.r;
  395.     z_2.r = -jay.r, z_2.i = -jay.i;
  396.     z_1.r = z_2.r * p.r - z_2.i * p.i, z_1.i = z_2.r * p.i + z_2.i * p.r;
  397.     parti = z_1.r;
  398. /* Computing 2nd power */
  399.     d_1 = partr;
  400. /* Computing 2nd power */
  401.     d_2 = parti;
  402.     if (sqrt(d_1 * d_1 + d_2 * d_2) > amax) {
  403.     goto L300;
  404.     }
  405.     i_4 = mmp - 1;
  406.     i_1 = mmp - 1;
  407.     cc[i_4].r = bb[i_1].r, cc[i_4].i = bb[i_1].i;
  408.     mdec = *mm - 1;
  409.     i_4 = mdec;
  410.     for (j = 1; j <= i_4; ++j) {
  411.     k = *mm - j + 1;
  412. /* L560: */
  413.     i_1 = k - 1;
  414.     i_2 = k;
  415.     z_2.r = w.r * cc[i_2].r - w.i * cc[i_2].i, z_2.i = w.r * cc[i_2].i + 
  416.         w.i * cc[i_2].r;
  417.     i_3 = k - 1;
  418.     z_1.r = z_2.r + bb[i_3].r, z_1.i = z_2.i + bb[i_3].i;
  419.     cc[i_1].r = z_1.r, cc[i_1].i = z_1.i;
  420.     }
  421.     pp.r = cc[1].r, pp.i = cc[1].i;
  422.     partr = pp.r;
  423.     z_2.r = -jay.r, z_2.i = -jay.i;
  424.     z_1.r = z_2.r * pp.r - z_2.i * pp.i, z_1.i = z_2.r * pp.i + z_2.i * pp.r;
  425.     parti = z_1.r;
  426. /* Computing 2nd power */
  427.     d_1 = partr;
  428. /* Computing 2nd power */
  429.     d_2 = parti;
  430.     if (sqrt(d_1 * d_1 + d_2 * d_2) < amin) {
  431.     goto L300;
  432.     }
  433.     partr = p.r;
  434.     z_2.r = -jay.r, z_2.i = -jay.i;
  435.     z_1.r = z_2.r * p.r - z_2.i * p.i, z_1.i = z_2.r * p.i + z_2.i * p.r;
  436.     parti = z_1.r;
  437. /* Computing 2nd power */
  438.     d_1 = partr;
  439. /* Computing 2nd power */
  440.     d_2 = parti;
  441.     size = sqrt(d_1 * d_1 + d_2 * d_2);
  442. /*       test for convergence of polishing */
  443.     if (size <= *eps) {
  444.     goto L501;
  445.     }
  446.     z_div(&z_2, &p, &pp);
  447.     z_1.r = w.r - z_2.r, z_1.i = w.i - z_2.i;
  448.     w.r = z_1.r, w.i = z_1.i;
  449.     goto L510;
  450. /*        deflate */
  451. L501:
  452.     i_1 = mp - 1;
  453.     i_2 = mp;
  454.     b[i_1].r = a[i_2], b[i_1].i = 0.;
  455.     i_1 = m;
  456.     for (j = 1; j <= i_1; ++j) {
  457.     k = m - j + 1;
  458. /* L830: */
  459.     i_2 = k - 1;
  460.     i_3 = k;
  461.     z_2.r = z.r * b[i_3].r - z.i * b[i_3].i, z_2.i = z.r * b[i_3].i + z.i 
  462.         * b[i_3].r;
  463.     i_4 = k;
  464.     z_1.r = z_2.r + a[i_4], z_1.i = z_2.i;
  465.     b[i_2].r = z_1.r, b[i_2].i = z_1.i;
  466.     }
  467.     p.r = b[0].r, p.i = b[0].i;
  468.     i_2 = m;
  469.     rootr[i_2] = w.r;
  470.     i_2 = m;
  471.     z_2.r = -jay.r, z_2.i = -jay.i;
  472.     z_1.r = z_2.r * w.r - z_2.i * w.i, z_1.i = z_2.r * w.i + z_2.i * w.r;
  473.     rooti[i_2] = z_1.r;
  474.     --m;
  475.     --mp;
  476.     z_2.r = -jay.r, z_2.i = -jay.i;
  477.     z_1.r = z_2.r * w.r - z_2.i * w.i, z_1.i = z_2.r * w.i + z_2.i * w.r;
  478.     parti = z_1.r;
  479.     if (abs(parti) > sqteps) {
  480.     goto L140;
  481.     }
  482. /*        real root */
  483.     rooti[m + 1] = 0.;
  484.     i_2 = mp;
  485.     for (j = 1; j <= i_2; ++j) {
  486. /* L100: */
  487.     i_3 = j;
  488.     i_4 = j;
  489.     a[i_3] = b[i_4].r;
  490.     }
  491.     goto L10;
  492. /*        complex root */
  493. L140:
  494.     partr = z.r;
  495.     z_2.r = -jay.r, z_2.i = -jay.i;
  496.     z_1.r = z_2.r * z.r - z_2.i * z.i, z_1.i = z_2.r * z.i + z_2.i * z.r;
  497.     parti = z_1.r;
  498.     z_2.r = parti * jay.r, z_2.i = parti * jay.i;
  499.     z_1.r = partr - z_2.r, z_1.i = -z_2.i;
  500.     z.r = z_1.r, z.i = z_1.i;
  501.     i_3 = mp - 1;
  502.     i_4 = mp;
  503.     c[i_3].r = b[i_4].r, c[i_3].i = b[i_4].i;
  504.     i_3 = m;
  505.     for (j = 1; j <= i_3; ++j) {
  506.     k = m - j + 1;
  507. /* L110: */
  508.     i_4 = k - 1;
  509.     i_2 = k;
  510.     z_2.r = z.r * c[i_2].r - z.i * c[i_2].i, z_2.i = z.r * c[i_2].i + z.i 
  511.         * c[i_2].r;
  512.     i_1 = k;
  513.     z_1.r = z_2.r + b[i_1].r, z_1.i = z_2.i + b[i_1].i;
  514.     c[i_4].r = z_1.r, c[i_4].i = z_1.i;
  515.     }
  516.     i_4 = m;
  517.     rootr[i_4] = w.r;
  518.     i_4 = m;
  519.     z_3.r = -jay.r, z_3.i = -jay.i;
  520.     z_2.r = z_3.r * w.r - z_3.i * w.i, z_2.i = z_3.r * w.i + z_3.i * w.r;
  521.     z_1.r = -z_2.r, z_1.i = -z_2.i;
  522.     rooti[i_4] = z_1.r;
  523.     --m;
  524.     --mp;
  525.     i_4 = mp;
  526.     for (j = 1; j <= i_4; ++j) {
  527. /* L130: */
  528.     i_2 = j;
  529.     i_1 = j;
  530.     a[i_2] = c[i_1].r;
  531.     }
  532.     goto L10;
  533. /*        report and return */
  534. L600:
  535.     real1 = (doublereal) kount;
  536.     real2 = (doublereal) (*mm);
  537.     temp = real1 / real2;
  538. /*     write(0,150)kount,temp */
  539. /* L150: */
  540. /*     write(0,151)newst,kerr */
  541. /* L151: */
  542.     return 0;
  543. } /* dproot_ */
  544.  
  545.