flower/polynomial.cc

   1 /*
   2    poly.cc -- routines for manipulation of polynomials in one var
   3
   4    (c) 1993--1999 Han-Wen Nienhuys <hanwen@cs.uu.nl>
   5  */
   6
   7
   8 #include "polynomial.hh"
   9
  10 /*
  11    Een beter milieu begint bij uzelf. Hergebruik!
  12
  13
  14    This was ripped from Rayce, a raytracer I once wrote.
  15 */
  16
  17 Real
  18 Polynomial::eval (Real x)const
  19 {
  20   Real p = 0.0;
  21
  22   // horner's scheme
  23   for (int i = coefs_.size (); i--; )
  24     p = x * p + coefs_[i];
  25
  26   return p;
  27 }
  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   // classicint 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         } else
  72           {
  73             base = multiply(base, base);
  74             e /= 2;
  75           }
  76     }
  77   return  dest;
  78 }
  79
  80 static Real const FUDGE = 1e-8;
  81
  82 void
  83 Polynomial::clean()
  84 {
  85   int i;
  86   for (i = 0; i <= degree (); i++)
  87     {
  88       if (abs(coefs_[i]) < FUDGE)
  89         coefs_[i] = 0.0;
  90     }
  91
  92   while (degree () > 0 && fabs (coefs_.top ()) < FUDGE * fabs (coefs_.top (1)))
  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 * p;
 241   Real D = q * q + cb_p;
 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_p));
 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       assert (fabs (eval (sol[i]) ) < 1e-8);
 277     }
 278
 279   return sol;
 280 }
 281
 282 Real
 283 Polynomial::lc () const
 284 {
 285   return coefs_.top();
 286 }
 287
 288 Real&
 289 Polynomial::lc ()
 290 {
 291   return coefs_.top ();
 292 }
 293
 294 int
 295 Polynomial::degree ()const
 296 {
 297   return coefs_.size () -1;
 298 }
 299 /*
 300    all roots of quadratic eqn.
 301  */
 302 Array<Real>
 303 Polynomial::solve_quadric()const
 304 {
 305   Array<Real> sol;
 306   /* normal form: x^2 + px + q = 0 */
 307   Real p = coefs_[1] / (2 * coefs_[2]);
 308   Real q = coefs_[0] / coefs_[2];
 309
 310   Real D = p * p - q;
 311
 312   if (D>0) {
 313     D = sqrt(D);
 314
 315     sol.push ( D - p);
 316     sol.push ( -D - p);
 317   }
 318   return sol;
 319 }
 320
 321 /* solve linear equation */
 322 Array<Real>
 323 Polynomial::solve_linear()const
 324 {
 325   Array<Real> s;
 326   if (coefs_[1])
 327     s.push ( -coefs_[0] / coefs_[1]);
 328   return s;
 329 }
 330
 331
 332 Array<Real>
 333 Polynomial::solve () const
 334 {
 335   Polynomial * me = (Polynomial*) this;
 336   me->clean ();
 337
 338   switch (degree ())
 339     {
 340     case 1:
 341       return solve_linear ();
 342     case 2:
 343       return solve_quadric ();
 344     case 3:
 345       return solve_cubic ();
 346     }
 347   assert (false);
 348   Array<Real> s;
 349   return s;
 350 }
 351
 352 void
 353 Polynomial:: operator *= (Polynomial const &p2)
 354 {
 355   *this = multiply (*this,p2);
 356 }
 357
 358 void
 359 Polynomial::operator += (Polynomial const &p)
 360 {
 361   *this = add( *this, p);
 362 }
 363
 364 void
 365 Polynomial::operator -= (Polynomial const &p)
 366 {
 367   *this = subtract(*this, p);
 368 }