home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 3 / RISC_DISC_3.iso / resources / etexts / gems / gemsv / ch5_4 / poly.cpp next >
C/C++ Source or Header  |  1995-03-04  |  3KB  |  131 lines
  1. /************************************************
  2.  * POLY.CPP
  3.  * Andreas Leipelt, "Ray Tracing a Swept Sphere"
  4.  * from "Graphics Gems", Academic Press
  5.  *
  6.  * Implementation of the polynomial class. The code is
  7.  * not complete ! You need to insert a root solver in
  8.  * the method 'root_between' .
  9.  */
  10.  
  11. #include <math.h>
  12. #include "poly.h"
  13.  
  14. // constructor of the polynomial class 
  15. polynomial::polynomial()
  16. {
  17.   deg = 0;
  18.   for (double *fp = &coef[MAX_DEGREE]; fp >= coef; fp--) *fp = 0.0;
  19. }
  20.  
  21. // evaluates the polynomial with Horner's scheme.
  22. double polynomial::eval(double x)
  23. {
  24.   double *fp = &coef[deg], val;
  25.   for (val = *fp--; fp >= coef; fp--) val = val*x + *fp;
  26.   return val;
  27. }
  28.  
  29. // returns the first derivative of the polynomial.
  30. polynomial polynomial::derivative()
  31. {
  32.   polynomial ret;
  33.  
  34.   if (!deg) return ret;
  35.   ret.deg = deg-1;
  36.   for (int i=0; i <= ret.deg; i++) ret.coef[i] = (i+1)*coef[i+1];
  37.   return ret;
  38. }
  39.  
  40. // returns the absolute minimum of the given polynomial in the
  41. // interval [a , b]
  42. double polynomial::min(double a, double b)
  43. {
  44.   double roots[MAX_DEGREE], tmp, Min = eval(a);
  45.  
  46.   int n = derivative().roots_between(a, b, roots);
  47.   roots[n] = b;
  48.   for (int i=0; i <= n; i++) {
  49.     tmp = eval(roots[i]);
  50.     if (tmp < Min) Min = tmp;
  51.   }
  52.   return Min;
  53. }
  54.  
  55. // returns the absolute maximum of the given polynomial in the
  56. // interval [a ; b]
  57. double polynomial::max(double a, double b)
  58. {
  59.   double roots[MAX_DEGREE], tmp, Max = eval(a);
  60.  
  61.   int n = derivative().roots_between(a, b, roots);
  62.   roots[n] = b;
  63.   for (int i=0; i <= n; i++) {
  64.     tmp = eval(roots[i]);
  65.     if (tmp > Max) Max = tmp;
  66.   }
  67.   return Max;
  68. }
  69.  
  70. int polynomial::roots_between(double a, double b, double *roots)
  71. {
  72.   // This function should return the number of roots between
  73.   // a  and  b  and the array 'roots' should contain these roots.
  74.   // Refer to Hook and McAree, "Using Sturm Sequences to Bracket
  75.   // Real Roots of Polynomial Equations" in "Graphics Gems I"
  76.   return 0;
  77. }
  78.  
  79. polynomial operator+(polynomial& p, polynomial& q)
  80. {
  81.   polynomial sum;
  82.  
  83.   if (p.deg < q.deg) sum.deg = q.deg;
  84.   else sum.deg = p.deg;
  85.   for (int i=0; i <= sum.deg; i++)
  86.      sum.coef[i] = p.coef[i] + q.coef[i];
  87.   if (p.deg == q.deg) {
  88.     while (sum.deg > -1 && fabs(sum.coef[sum.deg]) < polyeps)
  89.        sum.coef[sum.deg--] = 0.0;
  90.     if (sum.deg < 0) sum.deg = 0;
  91.   }
  92.   return sum;
  93. }
  94.  
  95. polynomial operator-(polynomial& p, polynomial& q)
  96. {
  97.   polynomial dif;
  98.  
  99.   if (p.deg < q.deg) dif.deg = q.deg;
  100.   else dif.deg = p.deg;
  101.   for (int i=0; i <= dif.deg; i++)
  102.      dif.coef[i] = p.coef[i] - q.coef[i];
  103.   if (p.deg == q.deg) {
  104.     while (dif.deg > -1 && fabs(dif.coef[dif.deg]) < polyeps)
  105.        dif.coef[dif.deg--] = 0.0;
  106.     if (dif.deg < 0) dif.deg = 0;
  107.   }
  108.   return dif;
  109. }
  110.  
  111. polynomial operator*(polynomial& p, polynomial& q)
  112. {
  113.   polynomial prod;
  114.  
  115.   prod.deg = p.deg + q.deg;
  116.   for (int i=0; i <= p.deg; i++)
  117.     for (int j=0; j <= q.deg; j++)
  118.        prod.coef[i+j] += p.coef[i]*q.coef[j];
  119.   return prod;
  120. }
  121.  
  122. polynomial operator*(double s, polynomial& p)
  123. {
  124.   polynomial scale;
  125.  
  126.   if (s == 0.0) return scale;
  127.   scale.deg = p.deg;
  128.   for (int i=0; i <= p.deg; i++) scale.coef[i] = s*p.coef[i];
  129.   return scale;
  130. }
  131.