flower/polynomial.cc

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