home *** CD-ROM | disk | FTP | other *** search
/ Stars of Shareware: Raytrace & Morphing / SOS-RAYTRACE.ISO / programm / source / rayce27s / poly.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-02-02  |  3.5 KB  |  192 lines
  1. /*
  2.  * poly.c -- routines for manipulation of polynomials in one var
  3.  * 
  4.  * (c) 1993, 1994 Han-Wen Nienhuys
  5.  * 
  6.  * This program is free software; you can redistribute it and/or modify it
  7.  * under the terms of the GNU General Public License as published by the
  8.  * Free Software Foundation;
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License along
  16.  * with this program; if not, write to the Free Software Foundation, Inc.,
  17.  * 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include "ray.h"
  21. #include "proto.h"
  22.  
  23. /* fudge factor for small coeffs */
  24. #define    FUDGE        1.0e-12
  25.  
  26.  
  27. PUBLIC void
  28. poly_copy(spoly *dst, spoly *src)
  29. {
  30.     int             i;
  31.  
  32.     dst->ord = src->ord;
  33.  
  34.     for (i = 0; i <= src->ord; i++)
  35.     dst->coef[i] = src->coef[i];
  36. }
  37.  
  38.  
  39. PUBLIC void
  40. poly_multiply(spoly *dest, spoly *p1, spoly *p2)
  41. {
  42.     int             i,
  43.                     j;
  44.  
  45.     dest->ord = p1->ord + p2->ord;
  46.     if (dest->ord > MAX_ORDER)
  47.     errormsg("MAX_ORDER exceeded");
  48.  
  49.     for (i = 0; i <= dest->ord; i++) {
  50.     dest->coef[i] = 0;
  51.     for (j = 0; j <= i; j++)
  52.         if (i - j <= p2->ord && j <= p1->ord)
  53.         dest->coef[i] += p1->coef[j] * p2->coef[i - j];
  54.     }
  55.  
  56. }
  57.  
  58.  
  59. PUBLIC void
  60. poly_power(spoly *dst, int exponent, spoly *src)
  61. {
  62.     spoly           temp;
  63.     int             e;
  64.  
  65.     dst->ord = 0;
  66.     dst->coef[0] = 1;
  67.     e = exponent;
  68.  
  69.     /* classic int power. invariant: dst^exponent = dst*src ^ e */
  70.     while (e > 0) {
  71.     if (e % 2) {
  72.         poly_multiply(&temp, dst, src);
  73.         poly_copy(dst, &temp);
  74.         e--;
  75.     } else {
  76.         poly_multiply(&temp, src, src);
  77.         poly_copy(src, &temp);
  78.         e /= 2;
  79.     }
  80.     }
  81.     return;
  82. }
  83.  
  84.  
  85. PUBLIC void
  86. poly_clean(spoly *p)
  87. {
  88.     int             i;
  89.  
  90.     for (i = 0; i <= p->ord; i++) {
  91.     if (ABS(p->coef[i]) < FUDGE)
  92.         p->coef[i] = 0.0;
  93.     }
  94.  
  95.     for (i = p->ord + 1; i <= MAX_ORDER; i++)
  96.     p->coef[i] = 0;
  97.  
  98.     while (p->ord > 0 && p->coef[p->ord] == 0.0)
  99.     p->ord--;
  100. }
  101.  
  102.  
  103. PUBLIC void
  104. poly_add(spoly *dst, spoly *p1, spoly *p2)
  105. {
  106.     int             i,
  107.                     tempord;
  108.     double          temp;
  109.  
  110.  
  111.     tempord = MAX(p1->ord, p2->ord);
  112.  
  113.     for (i = 0; i <= tempord; i++) {
  114.     temp = 0.0;
  115.     if (i <= p1->ord)
  116.         temp += p1->coef[i];
  117.     if (i <= p2->ord)
  118.         temp += p2->coef[i];
  119.     dst->coef[i] = temp;
  120.     }
  121.     dst->ord = tempord;
  122. }
  123.  
  124.  
  125. PUBLIC void
  126. poly_scalmult(spoly *d, double fact)
  127. {
  128.     int             i;
  129.  
  130.     for (i = 0; i <= d->ord; i++)
  131.     d->coef[i] *= fact;
  132. }
  133.  
  134.  
  135. PUBLIC void
  136. poly_sub(spoly *dst, spoly *p1, spoly *p2)
  137. {
  138.     int             i,
  139.                     tempord;
  140.     double          temp;
  141.  
  142.  
  143.     tempord = MAX(p1->ord, p2->ord);
  144.  
  145.     for (i = 0; i <= tempord; i++) {
  146.     temp = 0.0;
  147.     if (i <= p1->ord)
  148.         temp += p1->coef[i];
  149.     if (i <= p2->ord)
  150.         temp -= p2->coef[i];
  151.     dst->coef[i] = temp;
  152.     }
  153.     dst->ord = tempord;
  154.  
  155. }
  156.  
  157.  
  158. PUBLIC void
  159. poly_neg(spoly *dst, spoly *src)
  160. {
  161.     int             i;
  162.  
  163.     dst->ord = src->ord;
  164.  
  165.     for (i = 0; i <= src->ord; i++)
  166.     dst->coef[i] = -src->coef[i];
  167. }
  168.  
  169. /* print an spoly */
  170. PUBLIC void
  171. print_poly(spoly *p)
  172. {
  173. #ifdef DEBUG
  174.     int             d;
  175.  
  176.     d = p->ord;
  177.     printf("poly ord %d: ", d);
  178.     do {
  179.     if (p->coef[d]) {
  180.         printf("%lf", p->coef[d]);
  181.         if (d == 1)
  182.         printf(" t + ");
  183.         else if (d > 1) {
  184.         printf(" t^%d + ", d);
  185.         }
  186.     }
  187.     }
  188.     while (d-- > 0);
  189.     printf("\n");
  190. #endif
  191. }
  192.