#include <iostream.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int precision=8;
const int zmaxsize=21;
const int cmaxsize=21;
int qqq;

#include "genprec.h"
#include "cgenprec.h"
#include "apps.h"

extern int setexample();
extern int wilkinson();

int npoly;
int npolymax;           
int ijkl;

realmp apoly[zmaxsizepad];
realmp bpoly[zmaxsizepad];
complexmp rootsofpoly[cmaxsizepad];
realmp zrootsofpoly[zmaxsizepad];
complexmp dpoly[cmaxsizepad];
realmp drealpoly[cmaxsizepad];

int main()
  {int i;
   zsupe();
   zsuplogsp();
cout<<"  "<<'\n';
cout<<"17 example polynomials are solved, including the notorious Wilkinson"<<'\n';
cout<<"polynomial (example 12). An expert will want to vary the coefficient"<<'\n'; 
cout<<"of x^19 in the 10th significant figure and observe how the roots vary"<<'\n';
cout<<"In each example, the coefficients are listed first, followed by the"<<'\n';
cout<<"roots,followed by the coefficients calculated by multiplying out"<<'\n';
cout<<"(x-r1)*(x-r2)*...*(x-rn) which is a check on the roots."<<'\n';   
cout<<"  "<<'\n';
                    
   for (ijkl=1;ijkl<20;ijkl++)  
     {if (ijkl==12)
        {wilkinson();
         goto aroundsetex; 
        } 
      setexample(); 
aroundsetex:
      rootsofpoly<solvpoly(npoly,apoly);
      cout<<'\n';
      cout<<"for example"<<'\t'<<ijkl<<'\t'<<"the coefficients a[npoly],...,a[0] are"<<'\n';
      cout<<'\n';
      for (i=npoly;i>=0;i--)
        {printmp(apoly[i]);}      //changed
      cout<<'\n';
      cout<<"the complexified roots are:"<<'\n';
      for (i=1;i<npoly+1;i++)
        {printmp(rootsofpoly[i]);}      //these are ok

      dpoly<cconstructpoly(npoly,rootsofpoly);  //where is cconstructpoly
//      zmessage="the coefficients calculated from the roots are";
//      zbreak();                //could print like this but the
//      for (i=npoly;i>=0;i--)   //coefficients are real
//        {printmp(dpoly[i]);}  
//      zmessage="next example:";
//      zbreak();

      zmessage="the coefficients calculated from the roots are";
      zbreak();
      drealpoly<realify(npoly,dpoly);
      for (i=npoly;i>=0;i--)
        {printmp(drealpoly[i]);}
      zmessage="next example:";
      zbreak();
     }   
   return 0;
  }

int setexample()   //a number of examples
 { int i;
   if (ijkl==1)
     {npoly=2;
      apoly[2]=1;     //x^2+2x+2
      apoly[1]=2;
      apoly[0]=2;
     }
   if (ijkl==2)
     {npoly=8;
      apoly[8]=1;
      apoly[7]=0;
      apoly[6]=0;
      apoly[5]=0;
      apoly[4]=68;
      apoly[3]=0;
      apoly[2]=0;
      apoly[1]=0;
      apoly[0]=256;
     }
   if (ijkl==3)
     {npoly=3;
      apoly[3]=1;
      apoly[2]=-1;
      apoly[1]=1;
      apoly[0]=-1;
     }
   if (ijkl==4)
     {npoly=5;
      apoly[5]=1;
      apoly[4]=1;
      apoly[3]=0;
      apoly[2]=-10;
      apoly[1]=24;
      apoly[0]=-16;
     }
   if (ijkl==5)
     {npoly=6;
      apoly[6]=1;
      apoly[5]=0;
      apoly[4]=-1;
      apoly[3]=-10;
      apoly[2]=34;
      apoly[1]=-40;
      apoly[0]=16;
     }
   if (ijkl==6)
     {npoly=6;
      apoly[6]=1;
      apoly[5]=2;
      apoly[4]=1;
      apoly[3]=-10;
      apoly[2]=14;
      apoly[1]=8;
      apoly[0]=-16;
     }
   if (ijkl==7)
     {npoly=7;
      zconinttomp(1);
      apoly[7]=1;
      apoly[6]=-1;
      apoly[5]=-1;
      apoly[4]=-9;
      apoly[3]=44;
      apoly[2]=-74;
      apoly[1]=56;
      apoly[0]=-16;
     }
   if (ijkl==8)
     {npoly=2;
      apoly[2]=1;
      apoly[1]=-3;
      apoly[0]=2;
     }
   if (ijkl==9)
     {npoly=2;
      apoly[2]=1;
      apoly[1]=-2;
      apoly[0]=1;
     }
   if (ijkl==10)
     {npoly=3;
      apoly[3]=1;
      apoly[2]=6;
      apoly[1]=12;
      apoly[0]=8;
     }
   if (ijkl==11)
     {npoly=3;
      apoly[3]=1;
      apoly[2]=-6;
      apoly[1]=11;
      apoly[0]=-6;
     } 
   if (ijkl==13)
     {npoly=4;
                               //x^4 + 5 * x^2 + 4 = 0
      apoly[4]=1;              //x = +/- i , +/- 2i
      apoly[3]=0;
      apoly[2]=5;
      apoly[1]=0;
      apoly[0]=4;
     }
   if (ijkl==14)    
     {npoly=8;
      apoly[8]=1;
      apoly[7]=0;
      apoly[6]=30;
      apoly[5]=0;
      apoly[4]=273;
      apoly[3]=0;
      apoly[2]=820;
      apoly[1]=0;
      apoly[0]=576;
     }
   if (ijkl==15)
     {npoly=4;
      apoly[4]=1;
      apoly[3]=-4;
      apoly[2]=11;
      apoly[1]=-14;
      apoly[0]=10;
     }
   if (ijkl==16)
     {npoly=4;
      apoly[4]=1;
      apoly[3]=-4;
      apoly[2]=8;
      apoly[1]=-8;
      apoly[0]=4;
     }
   if (ijkl==17)
     {npoly=8;
      apoly[8]=1;
      apoly[7]=0;
      apoly[6]=10;
      apoly[5]=0;
      apoly[4]=33;
      apoly[3]=0;
      apoly[2]=40;
      apoly[1]=0;
      apoly[0]=16;
     }
   if (ijkl==18)
     {npoly=4;
      apoly[4]=1;
      apoly[3]=-1.43855E00;
      apoly[2]=7.956803E-01;
      apoly[1]=-2.0091914E-01;
      apoly[0]=1.9366E-02;
     }

   if (ijkl==19)
     {npoly=8;
      apoly[8]=1;
      apoly[7]=0;
      apoly[6]=6;
      apoly[5]=0;
      apoly[4]=13;
      apoly[3]=0;
      apoly[2]=12;
      apoly[1]=0;
      apoly[0]=4;
     }


   return 0;
 }

int wilkinson()
  {int iw;
   npoly=20;
   for (iw=1;iw<npoly+1;iw++)
     {zrootsofpoly[iw]=iw;    //the zeros are 1,2,....,20
      ymess[iw]=2;           //first construct the poly and
     }                       //then try to rediscover the roots
   ynpasses=20;
   apoly<constructpoly(npoly,zrootsofpoly);
   for (iw=0;iw<npoly+1;iw++)
     {printmp(apoly[iw]);}
   zmessage="coefficients in wilkinson polynomial";
   zbreak();
   return 0;  
  }