flower/polynomial.cc

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