/*********************************************************************
*   Test program for pegasus (6 functions defined in tfunc1.cpp)     *
* ------------------------------------------------------------------ *
* SAMPLE RUN:                                                        *
*                                                                    *
* File tpegasus.lst contains:                                        *
*                                                                    *
*--------------------------------------------------------------------*
* Pegasus method for searching real roots of real valued, nonlinear  *
* functions.                                                         *
*--------------------------------------------------------------------*
*                                                                    *
* f = 5 * x - exp (x)                                                *
* Starting value x0 =  0.000000                                      *
* Return code       =  1                                             *
* Root              =   0.2591711018190737                           *
* Function value    = -2.220446e-016                                 *
* Iterations        =  6                                             *
*                                                                    *
* f = (((((x-6)*x+15)*x-20)*x+15)*x-6)*x+1                           *
* Starting value x0 =  0.500000                                      *
* Return code       = -1                                             *
* Root              =   1.5000000000000000                           *
* Function value    =  1.562500e-002                                 *
* Iterations        =  0                                             *
*                                                                    *
* f = sin (x)                                                        *
* Starting value x0 =  3.000000                                      *
* Return code       =  1                                             *
* Root              =   3.1415926535897931                           *
* Function value    =  1.224606e-016                                 *
* Iterations        =  6                                             *
*                                                                    *
* f = 1 + sin (x)                                                    *
* Starting value x0 = -2.000000                                      *
* Return code       = -1                                             *
* Root              =  -1.5000000000000000                           * 
* Function value    =  2.505013e-003                                 *
* Iterations        =  0                                             *
*                                                                    *
* f = exp(x) -(1.0 + x + x*x*0.5)                                    *
* Starting value x0 =  2.000000                                      *
* Return code       = -1                                             *
* Root              =   3.0000000000000000                           * 
* Function value    =  1.158554e+001                                 *
* Iterations        =  0                                             *
*                                                                    *
* f = (x-1.0)*(x-1.0)*( sin(PI*x)  - log(fabs(2.0*x/(x+1.0)))        *
* Starting value x0 =  2.000000                                      *
* Return code       = -1                                             *
* Root              =   3.0000000000000000                           *
* Function value    = -1.621860e+000                                 *
* Iterations        =  0                                             *
*                                                                    *
*--------------------------------------------------------------------*
* 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)         *
*********************************************************************/
#include <basis.h>
#include "tfunc1.h"   


int pegasus             /* Pegasus Method    .........................*/
            (
             REALFCT   fct,       /* Function ........................*/
             REAL *    x1,        /* Starting value 1 ................*/
             REAL *    x2,        /* Starting value 2 / solution .....*/
             REAL *    f2,        /* Function value at x2 ............*/
             int *     iter       /* Number of iterations ............*/
            );
  
           
void main()  {

  int i, rc, iter;
  REAL x1,x2, f;
  REAL (*fct) (REAL);
  char *text;
  FILE *fp;
  
  fp=fopen("tpegasus.lst","w"); // output file
  WriteHead(fp," Pegasus method for searching real roots of real valued, nonlinear functions"); 

  for (i = 1; i <= 6; i++)
  {
    switch (i)
    {
      case 1: x1 = ZERO;
		      x2 = ONE;
              fct  = f1;
              text = "f = 5 * x - exp (x)";
              break;
      case 2: x1 = 0.5;
		      x2 = 1.5;
              fct  = f2;
              text = "f = (((((x-6)*x+15)*x-20)*x+15)*x-6)*x+1";
              break;
      case 3: x1 = 3.0;
		      x2 = 3.5;
              fct  = f3;
              text = "f = sin (x)";
              break;
      case 4: x1 = -2.0;
		      x2 = -1.5;
              fct  = f4;
              text = "f = 1 + sin (x)";
              break;
      case 5: x1 = TWO;
		      x2 = 3.0; 
              fct  = f5;
              text = "f = exp(x) -(1.0 + x + x*x*0.5)";
              break;
      case 6: x1 = TWO;
		      x2 = 3.0;
              fct  = f6;
              text = "f = (x-1.0)*(x-1.0)*( sin(PI*x)  - "
                     "log(fabs(2.0*x/(x+1.0)))";
              break;
      default:fprintf(fp," *** Example not registered\n\n"); 
		      return;
    }

    fprintf (fp," %s\n", text);
    fprintf (fp," Starting value x0 = ");
    fprintf (fp,FORMAT_LF, x1);
    fprintf (fp,"\n");

    rc = pegasus (fct, &x1, &x2, &f, &iter);

    fprintf (fp," Return code       = % d\n", rc);
    fprintf (fp," Root              = "); fprintf (fp,FORMAT_2016LF, x2);
    fprintf (fp,"\n Function value    = "); fprintf (fp,FORMAT_LE, f);
    fprintf (fp,"\n Iterations        = % d\n\n", iter);
  } 
  WriteEnd(fp);
  fclose(fp);

  printf("\nResults in tpegasus.lst.\n\n");

}

//end of file tpegasus.cpp