'************************************************
'* Roots of a polynomial by Bernouilli's method *
'* -------------------------------------------- *
'* Reference: "Methodes de calcul numerique by  *
'*             Claude Nowakowski, PSI Editions, *
'*             France, 1981" [BIBLI 07].        *
'* -------------------------------------------- *
'* SAMPLE RUN:                                  *
'*                                              *
'* Roots of a polynomial by Bernouilli's method *
'*                                              *
'* Input order of polynomial: 4                 *
'*                                              *
'* Input coefficients of the polynomial:        *
'*                                              *
'*   A( 0 ) = ? 1                               *
'*   A( 1 ) = ? 2                               *
'*   A( 2 ) = ? -13                             *
'*   A( 3 ) = ? -14                             *
'*   A( 4 ) = ? 24                              *
'*                                              *
'* ===== R O O T S =====                        *
'*                                              *
'*   X1 =   -4.0000                             *
'*   X2 =    3.0000                             *
'*   X3 =   -2.0001                             *
'*   X4 =    1.0001                             *
'*                                              *
'************************************************
defint i-n
defdbl a-h,o-z

  cls
  print
  print " Roots of a polynomial by Bernouilli's method"
  print
  input " Input order of polynomial: ", nd
  print

  dim A(nd), Y(nd)

  print " Input coefficients of the polynomial:"
  print
  for k=0 to nd
    print "   A(";k;") = "; : input A(k)
  next k

  gosub 1000  'call Bernouilli subroutine
              'and print results
end

1000 'Bernouilli's subroutine
  print
  print " ===== R O O T S ====="
  print
  n=nd
  if ABS(A(0))<1e-12 then
    print " Error: A(0) must be <> 0 !"
    return
  end if
  for k=0 to n
    A(k)=A(k)/A(0)
  next k
  nroot=0
5 x1=0# : nroot=nroot+1
  for k=2 to n
    Y(k)=0
  next k
  Y(1)=1#
  for l=1 to 50
    s=0#
    for k=1 to n
      s=s+Y(k)*A(k)
    next k
    Y(0)=-s
    x=Y(0)/Y(1)
    er=ABS(x1-x)
    if l<n then goto 50
    if er<5e-5 then goto 100
50  x1=x
    for k=n to 1 step -1
      Y(k)=Y(k-1)
    next k
  next l
  print " Error: no convergence (error=";er;")"
  return
100 print using "   X# = ####.####"; nroot; x
  n=n-1
  if n=0 then return
  A(1)=A(1)+x
  for k=2 to n
    A(k)=A(k)+x*A(k-1)
  next k
  goto 5
return

'end of file bernou.bas