home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 3 / RISC_DISC_3.iso / resources / etexts / gems / gemsi / sturm / main.c next >
C/C++ Source or Header  |  1992-09-08  |  2KB  |  120 lines
  1. /*
  2. Using Sturm Sequences to Bracket Real Roots of Polynomial Equations
  3. by D.G. Hook and P.R. McAree
  4. from "Graphics Gems", Academic Press, 1990
  5. */
  6.  
  7. /*
  8.  * main.c
  9.  *
  10.  *    a sample driver program.
  11.  */
  12. #include <stdio.h>
  13. #include <math.h>
  14. #include "solve.h"
  15.  
  16. /*
  17.  * a driver program for a root solver.
  18.  */
  19. main()
  20. {
  21.     poly    sseq[MAX_ORDER];
  22.     double     min, max, roots[MAX_ORDER];
  23.     int        i, j, order, nroots, nchanges, np, atmin, atmax;
  24.  
  25.     /*
  26.      * get the details...
  27.      */
  28.  
  29.     printf("Please enter order of polynomial: ");
  30.     scanf("%d", &order);
  31.  
  32.     printf("\n");
  33.  
  34.     for (i = order; i >= 0; i--) {
  35.             printf("Please enter coefficient number %d: ", i);
  36.             scanf("%lf", &sseq[0].coef[i]);
  37.     }
  38.  
  39.     printf("\n");
  40.  
  41.     /*
  42.      * build the Sturm sequence
  43.      */
  44.     np = buildsturm(order, sseq);
  45.  
  46.     printf("Sturm sequence for:\n");
  47.  
  48.     for (i = order; i >= 0; i--)
  49.             printf("%lf ", sseq[0].coef[i]);
  50.  
  51.     printf("\n\n");
  52.  
  53.     for (i = 0; i <= np; i++) {
  54.             for (j = sseq[i].ord; j >= 0; j--)
  55.                 printf("%lf ", sseq[i].coef[j]);
  56.             printf("\n");
  57.     }
  58.  
  59.     printf("\n");
  60.  
  61.  
  62.     /* 
  63.      * get the number of real roots
  64.      */
  65.     nroots = numroots(np, sseq, &atmin, &atmax);
  66.  
  67.     if (nroots == 0) {
  68.             printf("solve: no real roots\n");
  69.             exit(0);
  70.     }
  71.  
  72.     /*
  73.      * calculate the bracket that the roots live in
  74.      */
  75.     min = -1.0;
  76.     nchanges = numchanges(np, sseq, min);
  77.     for (i = 0; nchanges != atmin && i != MAXPOW; i++) { 
  78.             min *= 10.0;
  79.             nchanges = numchanges(np, sseq, min);
  80.     }
  81.  
  82.     if (nchanges != atmin) {
  83.             printf("solve: unable to bracket all negative roots\n");
  84.             atmin = nchanges;
  85.     }
  86.  
  87.     max = 1.0;
  88.     nchanges = numchanges(np, sseq, max);
  89.     for (i = 0; nchanges != atmax && i != MAXPOW; i++) { 
  90.             max *= 10.0;
  91.             nchanges = numchanges(np, sseq, max);
  92.     }
  93.  
  94.     if (nchanges != atmax) {
  95.             printf("solve: unable to bracket all positive roots\n");
  96.             atmax = nchanges;
  97.     }
  98.  
  99.     nroots = atmin - atmax;
  100.  
  101.     /*
  102.      * perform the bisection.
  103.      */
  104.     sbisect(np, sseq, min, max, atmin, atmax, roots);
  105.  
  106.  
  107.     /*
  108.      * write out the roots...
  109.      */
  110.     if (nroots == 1) {
  111.             printf("\n1 distinct real root at x = %f\n", roots[0]);
  112.     } else {
  113.             printf("\n%d distinct real roots for x: ", nroots);
  114.  
  115.             for (i = 0; i != nroots; i++)
  116.                 printf("%f ", roots[i]);
  117.             printf("\n");
  118.     }
  119. }
  120.