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   /*
 173     warning ("x=%f is not a root of polynomial\n"
 174     "f (x)=%f, f' (x)=%f \n", x, f, d); */
 175 }
 176
 177 void
 178 Polynomial::check_sols (Array<Real> roots) const
 179 {
 180   for (int i = 0; i < roots.size (); i++)
 181     check_sol (roots[i]);
 182 }
 183
 184 Polynomial::Polynomial (Real a, Real b)
 185 {
 186   coefs_.push (a);
 187   if (b)
 188     coefs_.push (b);
 189 }
 190
 191 /* cubic root. */
 192 inline Real cubic_root (Real x)
 193 {
 194   if (x > 0.0)
 195     return pow (x, 1.0 / 3.0);
 196   else if (x < 0.0)
 197     return -pow (-x, 1.0 / 3.0);
 198   else
 199     return 0.0;
 200 }
 201
 202 static bool
 203 iszero (Real r)
 204 {
 205   return !r;
 206 }
 207
 208 Array<Real>
 209 Polynomial::solve_cubic ()const
 210 {
 211   Array<Real> sol;
 212
 213   /* normal form: x^3 + Ax^2 + Bx + C = 0 */
 214   Real A = coefs_[2] / coefs_[3];
 215   Real B = coefs_[1] / coefs_[3];
 216   Real C = coefs_[0] / coefs_[3];
 217
 218   /*
 219    * substitute x = y - A/3 to eliminate quadric term: x^3 +px + q = 0
 220    */
 221
 222   Real sq_A = A *A;
 223   Real p = 1.0 / 3 * (-1.0 / 3 * sq_A + B);
 224   Real q = 1.0 / 2 * (2.0 / 27 * A *sq_A - 1.0 / 3 * A *B + C);
 225
 226   /* use Cardano's formula */
 227
 228   Real cb = p * p * p;
 229   Real D = q * q + cb;
 230
 231   if (iszero (D))
 232     {
 233       if (iszero (q)) { /* one triple solution */
 234         sol.push (0);
 235         sol.push (0);
 236         sol.push (0);
 237       }
 238       else {            /* one single and one double solution */
 239         Real u = cubic_root (-q);
 240
 241         sol.push (2 * u);
 242         sol.push (-u);
 243       }
 244     }
 245   else if (D < 0) {             /* Casus irreducibilis: three real solutions */
 246     Real phi = 1.0 / 3 * acos (-q / sqrt (-cb));
 247     Real t = 2 * sqrt (-p);
 248
 249     sol.push (t * cos (phi));
 250     sol.push (-t * cos (phi + M_PI / 3));
 251     sol.push (-t * cos (phi - M_PI / 3));
 252   }
 253   else {                        /* one real solution */
 254     Real sqrt_D = sqrt (D);
 255     Real u = cubic_root (sqrt_D - q);
 256     Real v = -cubic_root (sqrt_D + q);
 257
 258     sol.push (u + v);
 259   }
 260
 261   /* resubstitute */
 262   Real sub = 1.0 / 3 * A;
 263
 264   for (int i = sol.size (); i--;)
 265     {
 266       sol[i] -= sub;
 267
 268 #ifdef PARANOID
 269       assert (fabs (eval (sol[i])) < 1e-8);
 270 #endif
 271     }
 272
 273   return sol;
 274 }
 275
 276 Real
 277 Polynomial::lc () const
 278 {
 279   return coefs_.top ();
 280 }
 281
 282 Real &
 283 Polynomial::lc ()
 284 {
 285   return coefs_.top ();
 286 }
 287
 288 int
 289 Polynomial::degree ()const
 290 {
 291   return coefs_.size () -1;
 292 }
 293 /*
 294   all roots of quadratic eqn.
 295 */
 296 Array<Real>
 297 Polynomial::solve_quadric ()const
 298 {
 299   Array<Real> sol;
 300   /* normal form: x^2 + px + q = 0 */
 301   Real p = coefs_[1] / (2 * coefs_[2]);
 302   Real q = coefs_[0] / coefs_[2];
 303
 304   Real D = p * p - q;
 305
 306   if (D > 0)
 307     {
 308       D = sqrt (D);
 309
 310       sol.push (D - p);
 311       sol.push (-D - p);
 312     }
 313   return sol;
 314 }
 315
 316 /* solve linear equation */
 317 Array<Real>
 318 Polynomial::solve_linear ()const
 319 {
 320   Array<Real> s;
 321   if (coefs_[1])
 322     s.push (-coefs_[0] / coefs_[1]);
 323   return s;
 324 }
 325
 326 Array<Real>
 327 Polynomial::solve () const
 328 {
 329   Polynomial *me = (Polynomial *) this;
 330   me->clean ();
 331
 332   switch (degree ())
 333     {
 334     case 1:
 335       return solve_linear ();
 336     case 2:
 337       return solve_quadric ();
 338     case 3:
 339       return solve_cubic ();
 340     }
 341   Array<Real> s;
 342   return s;
 343 }
 344
 345 void
 346 Polynomial::operator *= (Polynomial const &p2)
 347 {
 348   *this = multiply (*this, p2);
 349 }
 350