flower/polynomial.cc

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