flower/polynomial.cc

   1 /*
   2   poly.cc -- routines for manipulation of polynomials in one var
   3
   4   (c) 1993--2005 Han-Wen Nienhuys <hanwen@xs4all.nl>
   5 */
   6
   7 #include "polynomial.hh"
   8
   9 #include <cmath>
  10 using namespace std;
  11
  12 /*
  13   Een beter milieu begint bij uzelf. Hergebruik!
  14
  15
  16   This was ripped from Rayce, a raytracer I once wrote.
  17 */
  18
  19 Real
  20 Polynomial::eval (Real x) const
  21 {
  22   Real p = 0.0;
  23
  24   // horner's scheme
  25   for (int i = coefs_.size (); i--;)
  26     p = x * p + coefs_[i];
  27
  28   return p;
  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     coefs_[i - 1] = coefs_[i] * i;
  53   coefs_.pop ();
  54 }
  55
  56 Polynomial
  57 Polynomial::power (int exponent, const Polynomial &src)
  58 {
  59   int e = exponent;
  60   Polynomial dest (1), base (src);
  61
  62   /*
  63     classic int power. invariant: src^exponent = dest * src ^ e
  64     greetings go out to Lex Bijlsma & Jaap vd Woude */
  65   while (e > 0)
  66     {
  67       if (e % 2)
  68         {
  69           dest = multiply (dest, base);
  70           e--;
  71         }
  72       else
  73
  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 void
  97 Polynomial::operator += (Polynomial const &p)
  98 {
  99   while (degree () < p.degree ())
 100     coefs_.push (0.0);
 101
 102   for (int i = 0; i <= p.degree (); i++)
 103     coefs_[i] += p.coefs_[i];
 104 }
 105
 106 void
 107 Polynomial::operator -= (Polynomial const &p)
 108 {
 109   while (degree () < p.degree ())
 110     coefs_.push (0.0);
 111
 112   for (int i = 0; i <= p.degree (); i++)
 113     coefs_[i] -= p.coefs_[i];
 114 }
 115
 116 void
 117 Polynomial::scalarmultiply (Real fact)
 118 {
 119   for (int i = 0; i <= degree (); i++)
 120     coefs_[i] *= fact;
 121 }
 122
 123 void
 124 Polynomial::set_negate (const Polynomial &src)
 125 {
 126   for (int i = 0; i <= src.degree (); i++)
 127     coefs_[i] = -src.coefs_[i];
 128 }
 129
 130 /// mod of #u/v#
 131 int
 132 Polynomial::set_mod (const Polynomial &u, const Polynomial &v)
 133 {
 134   (*this) = u;
 135
 136   if (v.lc () < 0.0)
 137     {
 138       for (int k = u.degree () - v.degree () - 1; k >= 0; k -= 2)
 139         coefs_[k] = -coefs_[k];
 140
 141       for (int k = u.degree () - v.degree (); k >= 0; k--)
 142         for (int j = v.degree () + k - 1; j >= k; j--)
 143           coefs_[j] = -coefs_[j] - coefs_[v.degree () + k] * v.coefs_[j - k];
 144     }
 145   else
 146
 147     {
 148       for (int k = u.degree () - v.degree (); k >= 0; k--)
 149         for (int j = v.degree () + k - 1; j >= k; j--)
 150           coefs_[j] -= coefs_[v.degree () + k] * v.coefs_[j - k];
 151     }
 152
 153   int k = v.degree () - 1;
 154   while (k >= 0 && coefs_[k] == 0.0)
 155     k--;
 156
 157   coefs_.set_size (1+ ((k < 0) ? 0 : k));
 158   return degree ();
 159 }
 160
 161 void
 162 Polynomial::check_sol (Real x) const
 163 {
 164   Real f = eval (x);
 165   Polynomial p (*this);
 166   p.differentiate ();
 167   Real d = p.eval (x);
 168
 169   if (abs (f) > abs (d) * FUDGE)
 170     ;
 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   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)
 244     {
 245     /* 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
 254     {
 255       /* one real solution */
 256       Real sqrt_D = sqrt (D);
 257       Real u = cubic_root (sqrt_D - q);
 258       Real v = -cubic_root (sqrt_D + q);
 259
 260       sol.push (u + v);
 261     }
 262
 263   /* resubstitute */
 264   Real sub = 1.0 / 3 * A;
 265
 266   for (int i = sol.size (); i--;)
 267     {
 268       sol[i] -= sub;
 269
 270 #ifdef PARANOID
 271       assert (fabs (eval (sol[i])) < 1e-8);
 272 #endif
 273     }
 274
 275   return sol;
 276 }
 277
 278 Real
 279 Polynomial::lc () const
 280 {
 281   return coefs_.top ();
 282 }
 283
 284 Real &
 285 Polynomial::lc ()
 286 {
 287   return coefs_.top ();
 288 }
 289
 290 int
 291 Polynomial::degree ()const
 292 {
 293   return coefs_.size () -1;
 294 }
 295 /*
 296   all roots of quadratic eqn.
 297 */
 298 Array<Real>
 299 Polynomial::solve_quadric ()const
 300 {
 301   Array<Real> sol;
 302   /* normal form: x^2 + px + q = 0 */
 303   Real p = coefs_[1] / (2 * coefs_[2]);
 304   Real q = coefs_[0] / coefs_[2];
 305
 306   Real D = p * p - q;
 307
 308   if (D > 0)
 309     {
 310       D = sqrt (D);
 311
 312       sol.push (D - p);
 313       sol.push (-D - p);
 314     }
 315   return sol;
 316 }
 317
 318 /* solve linear equation */
 319 Array<Real>
 320 Polynomial::solve_linear ()const
 321 {
 322   Array<Real> s;
 323   if (coefs_[1])
 324     s.push (-coefs_[0] / coefs_[1]);
 325   return s;
 326 }
 327
 328 Array<Real>
 329 Polynomial::solve () const
 330 {
 331   Polynomial *me = (Polynomial *) this;
 332   me->clean ();
 333
 334   switch (degree ())
 335     {
 336     case 1:
 337       return solve_linear ();
 338     case 2:
 339       return solve_quadric ();
 340     case 3:
 341       return solve_cubic ();
 342     }
 343   Array<Real> s;
 344   return s;
 345 }
 346
 347 void
 348 Polynomial::operator *= (Polynomial const &p2)
 349 {
 350   *this = multiply (*this, p2);
 351 }
 352