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