Função para Cálculo pelo Método Bairstow para Gretl

function [raiz_neg,raiz_pos,k,info]=bairstow(a,r,kmax,eps)

  // [raiz_neg,raiz_pos,k,info]=bairstow(a,r,kmax,eps)

  // Calcula uma aproxima��o para um par complexo de ra�zes do polin�mio p(x) cujos coeficientes

  // encontram-se no vetor "a"=[a_0,a_1,...,a_n], usando o m�todo de Bairstow, com estimativa inicial "r",

  // sujeita a |p(x)|<eps, desde que n�o exceda a "kmax" itera��es.

  // Os coeficientes devem ser armazenados em "a" na ordem a_0, a_1, ...   

  // Se info=0, localizou uma raiz; se info=-1, a derivada "df" � numericamente nula na itera��o "k";

  // se info=-2, |x_k-x_(k-1)|<del, mas |f(x)|>=eps; se info=-3, o n�mero m�ximo de itera��es foi excedido.

  // (C) Rudnei Dias da Cunha 2009

  info = 0;

  quase_zero = sqrt(number_properties("eps"));

  x = 2*real(r)

  y = -(real(r)^2+imag(r)^2)

  for k=1:kmax

    [matriz,termo_indep]=horner_quadratico(a,x,y)

        // b0 e b1 s�o -termo_indep(1), -termo_indep(2)

    if ((abs(termo_indep(1))<eps)&(abs(termo_indep(2))<eps))

               // resolve z^2-alpha*z-beta         

               [raiz_neg,raiz_pos] = quadratica_1(-x,-y)

               break

    else

        [L,U,P] = lu(matriz);

        delta = U\(L\P*termo_indep);

        x = x+delta(1)

        y = y+delta(2)

    end   

  end 

  if (k>=kmax)

    info = -3

  end

endfunction