-/* This is one of those 'mathemeticians should not write code' kind of
- cases. Newton's method of polishing roots is straightforward
- enough... except in those cases where it just fails in the real
- world. In our case below, we're worried about a local mini/maxima
- shooting a root estimation off to infinity, or the new estimation
- chaotically oscillating about convergence (shouldn't actually be a
- problem in our usage.
-
- Maehly's modification (zero suppression, to prevent two tenative
- roots from collapsing to the same actual root) similarly can
- temporarily shoot a root off toward infinity. It would come
- back... if it were not for the fact that machine representation has
- limited dynamic range and resolution. This too is guarded by
- limiting delta.
-
- Last problem is convergence criteria; we don't know what a 'double'
- is on our hardware/compiler, and the convergence limit is bounded
- by roundoff noise. So, we hack convergence:
-
- Require at most 1e-6 mean squared error for all zeroes. When
- converging, start the clock ticking at 1e-6; limit our polishing to
- as many more iterations as took us to get this far, 100 max.
-
- Past max iters, quit when MSE is no longer decreasing *or* we go
- below ~1e-20 MSE, whichever happens first. */
-
-static void Newton_Raphson_Maehly(float *a,int ord,float *r){
- int i, k, count=0, maxiter=0;
- double error=1.,besterror=1.;
- double *root=alloca(ord*sizeof(double));
-
- for(i=0; i<ord;i++) root[i] = 2.0 * (i+0.5) / ord - 1.0;
-
- while(error>1.e-20){
+/* Newton-Raphson-Maehly actually functioned as a decent root finder,
+ but there are root sets for which it gets into limit cycles
+ (exacerbated by zero suppression) and fails. We can't afford to
+ fail, even if the failure is 1 in 100,000,000, so we now use
+ Laguerre and later polish with Newton-Raphson (which can then
+ afford to fail) */
+
+#define EPSILON 10e-7
+static int Laguerre_With_Deflation(float *a,int ord,float *r){
+ int i,m;
+ double *defl=alloca(sizeof(*defl)*(ord+1));
+ for(i=0;i<=ord;i++)defl[i]=a[i];
+
+ for(m=ord;m>0;m--){
+ double new=0.f,delta;
+
+ /* iterate a root */
+ while(1){
+ double p=defl[m],pp=0.f,ppp=0.f,denom;
+
+ /* eval the polynomial and its first two derivatives */
+ for(i=m;i>0;i--){
+ ppp = new*ppp + pp;
+ pp = new*pp + p;
+ p = new*p + defl[i-1];
+ }
+
+ /* Laguerre's method */
+ denom=(m-1) * ((m-1)*pp*pp - m*p*ppp);
+ if(denom<0)
+ return(-1); /* complex root! The LPC generator handed us a bad filter */
+
+ if(pp>0){
+ denom = pp + sqrt(denom);
+ if(denom<EPSILON)denom=EPSILON;
+ }else{
+ denom = pp - sqrt(denom);
+ if(denom>-(EPSILON))denom=-(EPSILON);
+ }
+
+ delta = m*p/denom;
+ new -= delta;
+
+ if(delta<0.f)delta*=-1;
+
+ if(fabs(delta/new)<10e-12)break;
+ }
+
+ r[m-1]=new;
+
+ /* forward deflation */
+
+ for(i=m;i>0;i--)
+ defl[i-1]+=new*defl[i];
+ defl++;
+
+ }
+ return(0);
+}
+
+
+/* for spit-and-polish only */
+static int Newton_Raphson(float *a,int ord,float *r){
+ int i, k, count=0;
+ double error=1.f;
+ double *root=alloca(ord*sizeof(*root));
+
+ for(i=0; i<ord;i++) root[i] = r[i];
+
+ while(error>1e-20){