/**************************************************************************
*               Test program for Bauhuber's method                        *
* ----------------------------------------------------------------------- *
*  This program uses Bauhuber's Method to find all real or complex        *
*  roots of a polynomial of degree n:                                     *
*                                               n-1           n           *
*      P(x) = a[0] + a[1] * x + ... + a[n-1] * x    + a[n] * x ,          *
*                                                                         *
*  where a[i], i=0, ..., n, can be complex.                               *
* ----------------------------------------------------------------------- *
* SAMPLE RUN:                                                             *
* (Find all roots of polynomial of degree 4:                              *
*  P(x) = 0.000089248 - 0.04368 X + 2.9948 X2 - 6.798 X3 + X4)            *
*                                                                         *
* --------------------------------------------------------------------    *
*   Complex and Real Roots of a Polynomial (Bauhuber's method)            *  
* --------------------------------------------------------------------    *
*                                                                         *
* Polynomial coefficients:                                                *
*                                                                         *
*  a[ 0] =         0.0000892480         0.0000000000                      *
*  a[ 1] =        -0.0436800000         0.0000000000                      *
*  a[ 2] =         2.9948000000         0.0000000000                      *
*  a[ 3] =        -6.7980000000         0.0000000000                      *
*  a[ 4] =         1.0000000000         0.0000000000                      *
*                                                                         *
* Roots (without scaling)                                                 *
*     No. Real part               Imaginary part         function value   *
*                                                                         *
*     0   2.4537212366448121e-003  0.0000000000000000e+000  1.355253e-020 *  
*     1   1.2573281756029101e-002  0.0000000000000000e+000  4.065758e-020 *
*     2   4.5731893744203350e-001  0.0000000000000000e+000  4.540097e-018 * 
*     3   6.3256540595652933e+000  0.0000000000000000e+000  1.618457e-013 * 
*                                                                         *
*                                                                         *
* Roots (with scaling)                                                    *
*     No. Real part               Imaginary part         function value   *
*                                                                         *
*     0   2.4537212366448130e-003  0.0000000000000000e+000  0.000000e+000 *  
*     1   1.2573281756029099e-002  0.0000000000000000e+000  2.710505e-020 * 
*     2   4.5731893744203356e-001  0.0000000000000000e+000  1.382358e-018 * 
*     3   6.3256540595652924e+000  0.0000000000000000e+000  5.143054e-014 *
*                                                                         *
* --------------------------------------------------------------------    *
* Ref.: "Numerical algorithms with C, By Gisela Engeln-Muellges and       *
*        Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11].                  *
*                                                                         *
*                                       C++ Release By J-P Moreau, Paris. *
*                                              (www.jpmoreau.fr)          *
***************************************************************************
Note: to link with basis_r.cpp, fbauhube.cpp, vmblock.cpp
---------------------------------------------------------                */
#include <basis.h>       // for general declarations 
#include <vmblock.h>     // for memory dynamic allocations 

int  WriteHead1( char *);
int  WriteEnd1(void);
void SetVec (int, REAL *, REAL);

int bauhub              /* Bauhuber's method for complex polynomials .*/
           (
            int    real,          /* Are the coefficients real ? .....*/
            int    scale,         /* Scaling ? .......................*/
            int     n,            /* degree of polynomial ............*/
            REAL  ar[],           /* Real parts of coefficients ......*/
            REAL  ai[],           /* Imaginary parts, coefficients ...*/
            REAL  rootr[],        /* Real parts of roots .............*/
            REAL  rooti[],        /* Imaginary parts of roots ........*/
            REAL  absf[]          /* Absolute value of function values*/
           ); 
  

void main ()  {

  REAL *ar, *ai, *rootr, *rooti, *val;
  int i, j, n, rc, skala;
  void *vmblock;
 
  WriteHead1("  Complex and Real Roots of a Polynomial (Bauhuber's method)");

  n = 4;  //order of polynomial

  vmblock = vminit();
  ar    = (REAL *) vmalloc(vmblock, VEKTOR, n + 1, 0);
  ai    = (REAL *) vmalloc(vmblock, VEKTOR, n + 1, 0);
  rootr = (REAL *) vmalloc(vmblock, VEKTOR, n + 1, 0);
  rooti = (REAL *) vmalloc(vmblock, VEKTOR, n + 1, 0);
  val   = (REAL *) vmalloc(vmblock, VEKTOR, n + 1, 0);

  if (! vmcomplete(vmblock))
      LogError ("No Memory", 0, __FILE__, __LINE__);

  //define ar vector
  ar[0] =  0.000089248;
  ar[1] = -0.04368;
  ar[2] =  2.9948;
  ar[3] = -6.798;
  ar[4] =  ONE;

  //ai vector is null
  SetVec (n + 1, ai, ZERO);

  printf ("Polynomial coefficients:\n");
  for (i = 0; i <= n; i++) {
    printf("\n a[%2d] = ", i);
    printf (FORMAT_2010LF, ar[i]);
    printf (FORMAT_2010LF, ai[i]);
  }

  for (j = 0; j < 2; j++) {
    skala = j;
    rc = bauhub (0, skala, n, ar, ai, rootr, rooti, val);
    if (rc == 0)
    {
      if (skala == 0)
        printf ("\n\nRoots (without scaling)\n");
      else
        printf ("\n\nRoots (with scaling)\n");

      printf ("    No. Real part               Imaginary part         "
              "function value\n");

      for (i = 0; i < n; i++)
      {
        printf (" \n %4d  ", i);
        printf (FORMAT_2016LE, rootr[i]);
        printf (FORMAT_2016LE, rooti[i]);
        printf (FORMAT_LE, val[i]);
      }
      printf ("\n");
    }
    else
      LogError ("bauhube", rc, __FILE__, __LINE__);
  }  
  WriteEnd1();
}

//end of file tbauhube.cpp