/***********************************************
* Roots of a polynomial by Bernouilli's method *
* -------------------------------------------- *
* Reference: "Methodes de calcul numerique By  *
*             Claude Nowakowski, PSI Editions, *
*             France, 1981" [BIBLI 04].        *
*                                              *
*                C++ version by J-P Moreau     *
*                    (www.jpmoreau.fr)         *
* -------------------------------------------- *
* 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                             *
*                                              *
***********************************************/
#include <stdio.h>
#include <math.h>

#define SIZE  25

int    k,nd;
double temp,A[SIZE],Y[SIZE];


// Bernouilli's subroutine
void Bernouilli() {
// Labels: e5,e50,e100
  int k,l,n,nroot;
  double er,s,x,x1;
  printf("\n ===== R O O T S =====\n\n");
  n=nd;
  if (fabs(A[0])<1e-12) {
    printf(" Error: A(0) must be <> 0 !\n");
    return;
  }
  for (k=0; k<n+1; k++)  A[k]=A[k]/A[0];
  nroot=0;
e5: x1=0.0; nroot++;
  for (k=2; k<n+1; k++)  Y[k]=0.0;
  Y[1]=1.0;
  for (l=1; l<51; l++) {
    s=0.0;
    for (k=1; k<n+1; k++)  s += Y[k]*A[k];
    Y[0]=-s;
    x=Y[0]/Y[1];
    er=fabs(x1-x);
    if (l<n) goto e50;
    if (er<5e-5) goto e100;
e50: x1=x;
    for (k=n; k>0; k--)  Y[k]=Y[k-1];
  } //l loop
  printf("\n Error: no convergence (error=%e)\n",er);
  return;
e100: printf("   X%d = %9.4f\n",nroot,x);
  n--;
  if (n==0) return;
  A[1] += x;
  for (k=2; k<n+1; k++)  A[k] += x*A[k-1];
  goto e5;
}


void main()  {
  
  printf("\n Roots of a polynomial by Bernouilli's method\n\n");
  printf(" Input order of polynomial: "); scanf("%d",&nd);
  
  printf("\n Input coefficients of the polynomial:\n\n");
  
  for (k=0; k<nd+1; k++) {
    printf("   A(%d) = ",k); scanf("%lf",&temp);
    A[k] = temp;
  }

  // call Bernouilli procedure and printf results
  Bernouilli();

  printf("\n");
  printf("\n");

}

// end of file bernou.cpp