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