Re: Bairstow's Method

[Nicholas Jankowski]

On Sun, Sep 6, 2015 at 11:42 AM, superklutz <address@hidden> wrote:

Here is the code in case the link does not show up:

function [xre, xim, j, ear, eas] = Bairstow(a, deg, es, ri, si, maxit) %a =
matrix for coefficients
printf("Root evaluation using Bairstow's method.\n"); %Title
r = ri;
s = si;
b = zeros(1, deg+1);
c = zeros(1, deg);
xre = xim = zeros(deg, 1);
n = deg;
ier = 0;
ear = eas = 1;
j = 0;
ier = 0;
i = x = n-2;
do
  do
    j = j+1;
    b(1) = a(1);
    b(2) = a(2)+r.*b(n+1);
    c(1) = b(2);
    c(2) = b(n)+r.*c(n);
    do
      b(x) = a(x) + r.*b(x-1) + s.*b(x-2);
      c(i) = b(x) + r.*c(i-1) + s.*c(i-2);
      i=i+1;
      x=i+1;
    until(i == n && x == n+1)
    det =  c(n-1).^2 - c(n-2).*c(n);
    if (det != 0)
      dr = (-b(n).*c(n-1) + b(n+1).*c(n-2))./det;
      ds = (-b(n+1).*c(n-1) + b(n).*c(n))./det;
      r = r + dr;
      s = s + ds;
    else
      r = r+1;
      s = s+1;
      j = 0;
    endif
  until (ear<es && eas<es || j >= maxit) %Stopping criterion
[r1, i1, r2, i2] = quadform (r, s)
xre(n) = r1;
xim(n) = i1;
xre(n-1) = r2;
xim(n-1) = i2;
n = n-2;
  do
    a(x) = b(x-2);
    x = x-1;
  until (x == 1)
until (n < 3 || j >= maxit)
if (j < maxit)
  if (n == 2)
  r = -a(deg)./a(deg-1);
  s = -a(deg+1)./a(deg-1);
  [r1, i1, r2, i2] = quadform (r, s)
    xre(n) = r1;
    xim(n) = i1;
    xre(n-1) = r2;
    xim(n-1) = i2;
  else
  xre(n) = -a(deg+1)./a(deg);
  xim(n) = 0;
  endif
else
  ier = 1;
endif
printf("After %d iterations, the computed roots are as follows: \n", j);
%Prints the root and the error
xre
xim
clear
endfunction
function [r1, i1, r2, i2] = quadform (r, s)
dis = r^2 + 4*s;
if (dis>0)
  r1 = (r + sqrt(dis))/2;
  r2 = (r - sqrt(dis))/2;
  i1 = i2 = 0;
else
  r2 = r1 = r/2;
  i1 = sqrt(abs(dis))/2;
  i2 = -i1;
endif
endfunction

This is the error:

>> Bairstow([1, -3.5, 2.75, 2.125, -3.875, 1.25], 5, 0.00005, -1, -1, 100)
Root evaluation using Bairstow's method.
error: Bairstow: A(I): index out of bounds; value 7 out of bound 6
error: called from
    Bairstow at line 22 column 12
>>

in line 22, it throws an error when x is incremented to 7.  a(x) is then out of bounds as it only has 6 elements.

are you using Octave 4.0.0?  start trying to use the debugger. add breakpoints in the editor, and then step through the code testing values  leading up to the error. you can do the same without the gui if you're using an earlier version of Octave. add 'keyboard' commands to pause execution, and dbstep or dbcont commands to walk through the program.