0

I'm trying to implement robust / stable Laguerre's method. My code works for most of polynomials but for few it "fails". Respectively I don't know how to properly handle "corner" cases. Following code tries to find single root (F64 - 64bit float, C64 - 64bit complex):

    private static C64 GetSingle( C64 guess, int degree, C64 [] coeffs )
    {
        var i = 0;
        var x = guess;
        var n = (F64) degree;
        
        while( true )
        {
            ++Iters;
            
            var v0 = PolyUtils.Evaluate( x, coeffs, degree );
            
            if( v0.Abs() <= ACCURACY )
                break;
            
            var v1 = PolyUtils.EvaluateDeriv1( x, coeffs, degree );
            var v2 = PolyUtils.EvaluateDeriv2( x, coeffs, degree );
            
            var g  = v1 / v0;
            var gg = g * g;
            
            var h  = gg - ( v2 / v0 );
            var f  = C64.Sqrt(( n - 1.0 ) * ( n * h - gg ));
            
            var d0 = g - f;
            var d1 = g + f;
            
            var dx = d0.Abs() >= d1.Abs() ? ( n / d0 ) : ( n / d1 );
            
                x -= dx;
            
    //      if( dx.Abs() <= ACCURACY )
    //          break;
            
            if( ++i == ITERS_PER_ROOT ) // even after trying all guesses we didn't converted to the root (within given accuracy)
                break;
            
            if(( i & ( ITERS_PER_GUESS - 1 )) == 0 ) // didn't converge yet --> restart with different guess
            {
                x = GUESSES[ i / ITERS_PER_GUESS ];
            }
        }
        
        return x;
    }

At the end if it didn't found root it tries different guess, first quess (if not specified) is always 'zero'.

For example for 'x^4 + x^3 + x + 1' it founds 1st root '-1'. Deflates (divides) original poly by 'x + 1' so 2nd root search continues with polynomial 'x^3 + 1'. Again it starts with 'zero' as initial guess... but now both 1st and 2nd derivates are 'zero' which leads to 'zero' in 'd0' and 'd1'... ending by division-by-zero (and NaNs in root).

Another such example is 'x^5 - 1' - while searching for 1st root we again ends with zero derivates.

Can someone tell me how to handle such situations?

Should I just try another guess if derivates are 'zero'? I saw many implementation on net but none had such conditions so I don't know if I'm missing something.

Thank you

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51
Martin Capoušek
  • 3
  • 3
  • The usual recommendation is to guess the inner root radius and place the initial guess on the circle of that radius. It is well-known that the polynomials `x^n-1` or `x^n+1`, `n>3`, are the exceptions to the good behavior of the Laguerre method, they have a region around 0 where the iteration diverges or oscillates, this is (geometrically) reflected on the unit circle, so that convergence only happens from some annulus around the unit circle. – Lutz Lehmann Feb 01 '21 at 08:03
  • Thank you. For now I restart root search with different guess. – Martin Capoušek Feb 12 '21 at 12:05

0 Answers0