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