flower/polynomial.cc

   1 /*
   2    poly.cc -- routines for manipulation of polynomials in one var
   3
   4   (c)  1993--2003 Han-Wen Nienhuys <hanwen@cs.uu.nl>
   5  */
   6
   7 #include <math.h>
   8
   9 #include "polynomial.hh"
  10
  11 /*
  12    Een beter milieu begint bij uzelf. Hergebruik!
  13
  14
  15    This was ripped from Rayce, a raytracer I once wrote.
  16 */
  17
  18 Real
  19 Polynomial::eval (Real x)const
  20 {
  21   Real p = 0.0;
  22
  23   // horner's scheme
  24   for (int i = coefs_.size (); i--; )
  25     p = x * p + coefs_[i];
  26
  27   return p;
  28 }
  29
  30
  31 Polynomial
  32 Polynomial::multiply (const Polynomial & p1, const Polynomial & p2)
  33 {
  34   Polynomial dest;
  35
  36   int deg= p1.degree () + p2.degree ();
  37   for (int i = 0; i <= deg; i++)
  38     {
  39       dest.coefs_.push (0);
  40       for (int j = 0; j <= i; j++)
  41         if (i - j <= p2.degree () && j <= p1.degree ())
  42           dest.coefs_.top () += p1.coefs_[j] * p2.coefs_[i - j];
  43     }
  44
  45   return dest;
  46 }
  47
  48 void
  49 Polynomial::differentiate ()
  50 {
  51   for (int i = 1; i<= degree (); i++)
  52     {
  53       coefs_[i-1] = coefs_[i] * i;
  54     }
  55   coefs_.pop ();
  56 }
  57
  58 Polynomial
  59 Polynomial::power (int exponent, const Polynomial & src)
  60 {
  61   int e = exponent;
  62   Polynomial dest (1), base (src);
  63
  64   /*
  65     classic int power. invariant: src^exponent = dest * src ^ e
  66     greetings go out to Lex Bijlsma & Jaap vd Woude */
  67   while (e > 0)
  68     {
  69       if (e % 2)
  70         {
  71           dest = multiply (dest, base);
  72           e--;
  73         } else
  74           {
  75             base = multiply (base, base);
  76             e /= 2;
  77           }
  78     }
  79   return  dest;
  80 }
  81
  82 static Real const FUDGE = 1e-8;
  83
  84 void
  85 Polynomial::clean ()
  86 {
  87 /*
  88   We only do relative comparisons. Absolute comparisons break down in
  89   degenerate cases.  */
  90   while (degree () > 0 &&
  91  (fabs (coefs_.top ()) < FUDGE * fabs (coefs_.top (1))
  92           || !coefs_.top ()))
  93     coefs_.pop ();
  94 }
  95
  96
  97
  98 void
  99 Polynomial::operator += (Polynomial const &p)
 100 {
 101   while (degree () < p.degree())
 102     coefs_.push (0.0);
 103
 104   for (int i = 0; i <= p.degree(); i++)
 105     coefs_[i] += p.coefs_[i];
 106 }
 107
 108
 109 void
 110 Polynomial::operator -= (Polynomial const &p)
 111 {
 112   while (degree () < p.degree())
 113     coefs_.push (0.0);
 114
 115   for (int i = 0; i <= p.degree(); i++)
 116     coefs_[i] -= p.coefs_[i];
 117 }
 118
 119
 120 void
 121 Polynomial::scalarmultiply (Real fact)
 122 {
 123   for (int i = 0; i <= degree (); i++)
 124     coefs_[i] *= fact;
 125 }
 126
 127
 128
 129 void
 130 Polynomial::set_negate (const Polynomial & src)
 131 {
 132   for (int i = 0; i <= src.degree (); i++)
 133     coefs_[i] = -src.coefs_[i];
 134 }
 135
 136 /// mod of #u/v#
 137 int
 138 Polynomial::set_mod (const Polynomial &u, const Polynomial &v)
 139 {
 140  (*this) = u;
 141
 142   if (v.lc () < 0.0) {
 143     for (int k = u.degree () - v.degree () - 1; k >= 0; k -= 2)
 144       coefs_[k] = -coefs_[k];
 145
 146     for (int k = u.degree () - v.degree (); k >= 0; k--)
 147       for (int j = v.degree () + k - 1; j >= k; j--)
 148         coefs_[j] = -coefs_[j] - coefs_[v.degree () + k] * v.coefs_[j - k];
 149   } else {
 150     for (int k = u.degree () - v.degree (); k >= 0; k--)
 151       for (int j = v.degree () + k - 1; j >= k; j--)
 152         coefs_[j] -= coefs_[v.degree () + k] * v.coefs_[j - k];
 153   }
 154
 155   int k = v.degree () - 1;
 156   while (k >= 0 && coefs_[k] == 0.0)
 157     k--;
 158
 159   coefs_.set_size (1+ ((k < 0) ? 0 : k));
 160   return degree ();
 161 }
 162
 163 void
 164 Polynomial::check_sol (Real x) const
 165 {
 166   Real f=eval (x);
 167   Polynomial p (*this);
 168   p.differentiate ();
 169   Real d = p.eval (x);
 170
 171   if ( abs (f) > abs (d) * FUDGE)
 172     ;
 173   /*
 174     warning ("x=%f is not a root of polynomial\n"
 175     "f (x)=%f, f' (x)=%f \n", x, f, d); */
 176 }
 177
 178 void
 179 Polynomial::check_sols (Array<Real> roots) const
 180 {
 181   for (int i=0; i< roots.size (); i++)
 182     check_sol (roots[i]);
 183 }
 184
 185 Polynomial::Polynomial (Real a, Real b)
 186 {
 187   coefs_.push (a);
 188   if (b)
 189     coefs_.push (b);
 190 }
 191
 192 /* cubic root. */
 193 inline Real cubic_root (Real x)
 194 {
 195   if (x > 0.0)
 196     return pow (x, 1.0/3.0) ;
 197   else if (x < 0.0)
 198     return -pow (-x, 1.0/3.0);
 199   else
 200     return 0.0;
 201 }
 202
 203 static bool
 204 iszero (Real r)
 205 {
 206   return !r;
 207 }
 208
 209 Array<Real>
 210 Polynomial::solve_cubic ()const
 211 {
 212   Array<Real> sol;
 213
 214   /* normal form: x^3 + Ax^2 + Bx + C = 0 */
 215   Real A = coefs_[2] / coefs_[3];
 216   Real B = coefs_[1] / coefs_[3];
 217   Real C = coefs_[0] / coefs_[3];
 218
 219   /*
 220    * substitute x = y - A/3 to eliminate quadric term: x^3 +px + q = 0
 221    */
 222
 223   Real sq_A = A * A;
 224   Real p = 1.0 / 3 * (-1.0 / 3 * sq_A + B);
 225   Real q = 1.0 / 2 * (2.0 / 27 * A * sq_A - 1.0 / 3 * A * B + C);
 226
 227   /* use Cardano's formula */
 228
 229   Real cb = p * p * p;
 230   Real D = q * q + cb;
 231
 232   if (iszero (D)) {
 233     if (iszero (q)) {   /* one triple solution */
 234       sol.push (0);
 235       sol.push (0);
 236       sol.push (0);
 237     } else {            /* one single and one double solution */
 238       Real u = cubic_root (-q);
 239
 240       sol.push (2 * u);
 241       sol.push (-u);
 242     }
 243   } else if (D < 0) {           /* Casus irreducibilis: three real solutions */
 244     Real phi = 1.0 / 3 * acos (-q / sqrt (-cb));
 245     Real t = 2 * sqrt (-p);
 246
 247     sol.push (t * cos (phi));
 248     sol.push (-t * cos (phi + M_PI / 3));
 249     sol.push ( -t * cos (phi - M_PI / 3));
 250   } else {                      /* one real solution */
 251     Real sqrt_D = sqrt (D);
 252     Real u = cubic_root (sqrt_D - q);
 253     Real v = -cubic_root (sqrt_D + q);
 254
 255     sol.push ( u + v);
 256   }
 257
 258   /* resubstitute */
 259   Real sub = 1.0 / 3 * A;
 260
 261   for (int i = sol.size (); i--;)
 262     {
 263       sol[i] -= sub;
 264
 265 #ifdef PARANOID
 266       assert (fabs (eval (sol[i]) ) < 1e-8);
 267 #endif
 268     }
 269
 270   return sol;
 271 }
 272
 273 Real
 274 Polynomial::lc () const
 275 {
 276   return coefs_.top ();
 277 }
 278
 279 Real&
 280 Polynomial::lc ()
 281 {
 282   return coefs_.top ();
 283 }
 284
 285 int
 286 Polynomial::degree ()const
 287 {
 288   return coefs_.size () -1;
 289 }
 290 /*
 291    all roots of quadratic eqn.
 292  */
 293 Array<Real>
 294 Polynomial::solve_quadric ()const
 295 {
 296   Array<Real> sol;
 297   /* normal form: x^2 + px + q = 0 */
 298   Real p = coefs_[1] / (2 * coefs_[2]);
 299   Real q = coefs_[0] / coefs_[2];
 300
 301   Real D = p * p - q;
 302
 303   if (D>0) {
 304     D = sqrt (D);
 305
 306     sol.push ( D - p);
 307     sol.push ( -D - p);
 308   }
 309   return sol;
 310 }
 311
 312 /* solve linear equation */
 313 Array<Real>
 314 Polynomial::solve_linear ()const
 315 {
 316   Array<Real> s;
 317   if (coefs_[1])
 318     s.push ( -coefs_[0] / coefs_[1]);
 319   return s;
 320 }
 321
 322
 323 Array<Real>
 324 Polynomial::solve () const
 325 {
 326   Polynomial * me = (Polynomial*) this;
 327   me->clean ();
 328
 329   switch (degree ())
 330     {
 331     case 1:
 332       return solve_linear ();
 333     case 2:
 334       return solve_quadric ();
 335     case 3:
 336       return solve_cubic ();
 337     }
 338   Array<Real> s;
 339   return s;
 340 }
 341
 342 void
 343 Polynomial:: operator *= (Polynomial const &p2)
 344 {
 345   *this = multiply (*this,p2);
 346 }
 347
 348