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