test.c

Polynomial Root Finder is a reliable and fast C program (+ Matlab gateway) for finding all roots of

/*******************************************************************/
/*                                                                 */
/* test environment for the function null(); null() calculates     */
/* the roots of a polynomial with complex coefficients with a      */
/* combination of muller method and newton method                  */
/*                                                                 */
/*******************************************************************/
#define TEST
#include "header.h"










                      /* read coefficients stored in file FILENAME */
void datei_einlesen(int *grad)
/*int *grad;             delivers degree of polynomial             */
{
     FILE *fp;        /* filepointer                               */
     int  i;          /* counter                                   */

                      /* open file                                 */
     if ((fp = fopen(FILENAME,MODE)) == NULL) {
          printf("Can't open file %s!\n",FILENAME);
          return;
     } else {
                      /* read file: file format:                   */
                      /* degree | Flag | RE(a0) ... | RE(an)       */
                      /*                 IM(a0)...  | IM(an)       */
 
                      /* read degree                               */
          (void) fscanf(fp,"%d",grad);
                      /* read flag                                 */
          (void) fscanf(fp,"%d",&flag);
                      /* read coefficients                         */
          for (i=0; i<=(*grad); i++) {
               (void) fscanf(fp,"%le",&(p[i].r));
          if (flag) {
               for (i=0;i<=(*grad); i++)
                    (void) fscanf(fp,"%le",&(p[i].i));
               }
          }
                      /* close file                                */
          (void) fclose(fp);
     }
}

     


                      /* read coefficients from terminal           */
void tastatur_einlesen(int *grad)
/*int *grad;*/
{
     int i;           /* counter                                   */
     char ent;        /* read choice                               */ 

                      /* enter degree of polynomial                */
     printf("Please enter degree: ");
     (void) scanf("%d",grad);
     printf("Enter coefficients a0 ... an:\n");
                      /* Realteilabfrage                           */
     printf("real parts:\n");
     for (i=0; i<=(*grad); i++) {
          printf("a%d (real part): ",i);
          (void) scanf("%le",&(p[i].r));
     }
     printf("Exists also imaginary parts (y/n)? ");
     (void) scanf("%c",&ent);    /* clear keyboard                */
     (void) scanf("%c",&ent);    /* ans choice                    */
                                 /* complex coefficients          */
     if (ent == 'y' || ent == 'Y') {
          printf("imaginary parts:\n");
          for (i=0; i<=(*grad);i++) {
               printf("a%d (imaginary part): ",i);
               (void) scanf("%le",&(p[i].i));
          }
     } else flag = FALSE;        /* real coeff. => flag = FALSE   */
}





                      /* main loop                                 */
char main(void) {
     double maxerr;   /* max. error of all determined roots        */
     unsigned char   error;
     char   sig1,     /* signum of real part                       */
            sig2,     /* signum of imaginary part                  */
            ent;      /* read choice                               */ 
     int    degree,   /* degree of polynomial                      */
            i;        /* counter variable                          */
     FILE   *fp;      /* file pointer                              */

     printf("Input from file %s (y/n)?\n",FILENAME);
     (void) scanf("%c",&ent);
     if (ent == 'y' || ent == 'Y') 
          datei_einlesen(&degree);
     else
          tastatur_einlesen(&degree);

     error=null(p,pred,&degree,root,&maxerr,(unsigned char)flag);

     if (!error) {
                                /* generate matlab-file */
       fp=fopen("null.m","w"); 
       fprintf(fp,"re=[");
       for (i=1;i<=degree; i++) 
            fprintf(fp,"%.18e\n",root[i-1].r);
       fprintf(fp,"];\nim=[");
       for (i=1;i<=degree; i++)
            fprintf(fp,"%.18e\n",root[i-1].i);
       fprintf(fp,"];\ni=sqrt(-1);\nnull=re+i*im;\npn_plot(null,0)");
       (void) fclose(fp);
                               /* write stdoutput */ 
       for (i=1;i<=degree; i++) {
          sig1=(root[i-1].r>=0) ? ' ':'-';
          sig2=(root[i-1].i>=0) ? ' ':'-';
          printf("%3d. root:    %c%.18e     %c%.18e\n",i,
          sig1,fabs(root[i-1].r),sig2,fabs(root[i-1].i));
       }
       printf("max. error: %.18e\n",maxerr);
     } else {
          printf("Error %d occured!\n",(int)error);
          switch (error) {
          case 1:
               printf("Power of polynomial lower null!\n");
               break;
          case 2:
               printf("Polynomial is a null vector!\n");
               break;
          case 3:
               printf("Polynomial is a constant unequal null!\n");
               break;
          }
          return 1;
     }
     return 0;        /* no error occured       */
}