>>void laguer(Complex a[], int m, Complex *x, int *its) >>{ >> int iter, j; >> double abx, abp, abm, err; >> Complex dx,x1, b, d, f, g, h, sq, gp, gm, g2; >> static double frac[MR+1] = {0.0, 0.5, 0.25, 0.75, 0.13, >> 0.38, 0.62, 0.88, 1.0}; >> Complex xx, yy; >> >> for (iter = 1; iter <= MAXIT; iter++) { >> *its = iter; >> b = a[m]; >> err = abs(b); >> d = f = Complex(0.0, 0.0); >> abx = abs(*x); >> for (j=m-1;j>=0;j--) { >> f = f = (*x) * f + d; >> d = (*x) * d + b; >> b = (*x) * b + a[j]; >> err = abs(b) + abx * err; >> } >> err *= EPSS; >> >> if ( abs(b) <= err ) return; >> g = d/b; >> g2 = g*g; >> h = g2 - f/b * 2.0; >> sq = sqrt( (h*(m*1.0) - g2) * ((m-1)*1.0) ); >> gp = g + sq; >> gm = g - sq; >> abp = abs(gp); >> abm = abs(gm); >> if (abp < abm) gp=gm; >> >> dx = (((abp > 0.0) || (abm>0)) >> ? Complex(m,0)/gp >> : Complex( cos(iter), sin(iter)) * exp( log(1+abx) ) ); >> >> // : Complex( cos(iter), sin(iter) ) * exp( log(1+abx) ) ); >> >> x1 = *x - dx; >> if (real(*x) == real(x1) && imag(*x) == imag(x1)) return; >> if (iter % MT) *x = x1; >> else *x = *x - frac[iter/MT] * dx; >> } >> assert( 0 ); >> // exit(1); >>} >>>