flower/polynomial.cc

   1 /*
   2   This file is part of LilyPond, the GNU music typesetter.
   3
   4   Copyright (C) 1993--2009 Han-Wen Nienhuys <hanwen@xs4all.nl>
   5
   6   LilyPond is free software: you can redistribute it and/or modify
   7   it under the terms of the GNU General Public License as published by
   8   the Free Software Foundation, either version 3 of the License, or
   9   (at your option) any later version.
  10
  11   LilyPond is distributed in the hope that it will be useful,
  12   but WITHOUT ANY WARRANTY; without even the implied warranty of
  13   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14   GNU General Public License for more details.
  15
  16   You should have received a copy of the GNU General Public License
  17   along with LilyPond.  If not, see <http://www.gnu.org/licenses/>.
  18 */
  19
  20 #include "polynomial.hh"
  21
  22 #include "warn.hh"
  23
  24 #include <cmath>
  25
  26
  27 using namespace std;
  28
  29 /*
  30   Een beter milieu begint bij uzelf. Hergebruik!
  31
  32
  33   This was ripped from Rayce, a raytracer I once wrote.
  34 */
  35
  36 Real
  37 Polynomial::eval (Real x) const
  38 {
  39   Real p = 0.0;
  40
  41   // horner's scheme
  42   for (vsize i = coefs_.size (); i--;)
  43     p = x * p + coefs_[i];
  44
  45   return p;
  46 }
  47
  48 Polynomial
  49 Polynomial::multiply (const Polynomial &p1, const Polynomial &p2)
  50 {
  51   Polynomial dest;
  52
  53   int deg = p1.degree () + p2.degree ();
  54   for (int i = 0; i <= deg; i++)
  55     {
  56       dest.coefs_.push_back (0);
  57       for (int j = 0; j <= i; j++)
  58         if (i - j <= p2.degree () && j <= p1.degree ())
  59           dest.coefs_.back () += p1.coefs_[j] * p2.coefs_[i - j];
  60     }
  61
  62   return dest;
  63 }
  64
  65 void
  66 Polynomial::differentiate ()
  67 {
  68   for (int i = 1; i <= degree (); i++)
  69     coefs_[i - 1] = coefs_[i] * i;
  70   coefs_.pop_back ();
  71 }
  72
  73 Polynomial
  74 Polynomial::power (int exponent, const Polynomial &src)
  75 {
  76   int e = exponent;
  77   Polynomial dest (1), base (src);
  78
  79   /*
  80     classic int power. invariant: src^exponent = dest * src ^ e
  81     greetings go out to Lex Bijlsma & Jaap vd Woude */
  82   while (e > 0)
  83     {
  84       if (e % 2)
  85         {
  86           dest = multiply (dest, base);
  87           e--;
  88         }
  89       else
  90
  91         {
  92           base = multiply (base, base);
  93           e /= 2;
  94         }
  95     }
  96   return dest;
  97 }
  98
  99 static Real const FUDGE = 1e-8;
 100
 101 void
 102 Polynomial::clean ()
 103 {
 104   /*
 105     We only do relative comparisons. Absolute comparisons break down in
 106     degenerate cases.  */
 107   while (degree () > 0
 108          && (fabs (coefs_.back ()) < FUDGE * fabs (back (coefs_, 1))
 109              || !coefs_.back ()))
 110     coefs_.pop_back ();
 111 }
 112
 113 void
 114 Polynomial::operator += (Polynomial const &p)
 115 {
 116   while (degree () < p.degree ())
 117     coefs_.push_back (0.0);
 118
 119   for (int i = 0; i <= p.degree (); i++)
 120     coefs_[i] += p.coefs_[i];
 121 }
 122
 123 void
 124 Polynomial::operator -= (Polynomial const &p)
 125 {
 126   while (degree () < p.degree ())
 127     coefs_.push_back (0.0);
 128
 129   for (int i = 0; i <= p.degree (); i++)
 130     coefs_[i] -= p.coefs_[i];
 131 }
 132
 133 void
 134 Polynomial::scalarmultiply (Real fact)
 135 {
 136   for (int i = 0; i <= degree (); i++)
 137     coefs_[i] *= fact;
 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     {
 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     }
 162   else
 163
 164     {
 165       for (int k = u.degree () - v.degree (); k >= 0; k--)
 166         for (int j = v.degree () + k - 1; j >= k; j--)
 167           coefs_[j] -= coefs_[v.degree () + k] * v.coefs_[j - k];
 168     }
 169
 170   int k = v.degree () - 1;
 171   while (k >= 0 && coefs_[k] == 0.0)
 172     k--;
 173
 174   coefs_.resize (1+ ((k < 0) ? 0 : k));
 175   return degree ();
 176 }
 177
 178 void
 179 Polynomial::check_sol (Real x) const
 180 {
 181   Real f = eval (x);
 182   Polynomial p (*this);
 183   p.differentiate ();
 184   Real d = p.eval (x);
 185
 186   if (abs (f) > abs (d) * FUDGE)
 187     programming_error ("not a root of polynomial\n");
 188 }
 189
 190 void
 191 Polynomial::check_sols (vector<Real> roots) const
 192 {
 193   for (vsize i = 0; i < roots.size (); i++)
 194     check_sol (roots[i]);
 195 }
 196
 197 Polynomial::Polynomial (Real a, Real b)
 198 {
 199   coefs_.push_back (a);
 200   if (b)
 201     coefs_.push_back (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   return 0.0;
 212 }
 213
 214 static bool
 215 iszero (Real r)
 216 {
 217   return !r;
 218 }
 219
 220 vector<Real>
 221 Polynomial::solve_cubic ()const
 222 {
 223   vector<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     {
 245       if (iszero (q)) { /* one triple solution */
 246         sol.push_back (0);
 247         sol.push_back (0);
 248         sol.push_back (0);
 249       }
 250       else {            /* one single and one double solution */
 251         Real u = cubic_root (-q);
 252
 253         sol.push_back (2 * u);
 254         sol.push_back (-u);
 255       }
 256     }
 257   else if (D < 0)
 258     {
 259     /* Casus irreducibilis: three real solutions */
 260       Real phi = 1.0 / 3 * acos (-q / sqrt (-cb));
 261       Real t = 2 * sqrt (-p);
 262
 263       sol.push_back (t * cos (phi));
 264       sol.push_back (-t * cos (phi + M_PI / 3));
 265       sol.push_back (-t * cos (phi - M_PI / 3));
 266     }
 267   else
 268     {
 269       /* one real solution */
 270       Real sqrt_D = sqrt (D);
 271       Real u = cubic_root (sqrt_D - q);
 272       Real v = -cubic_root (sqrt_D + q);
 273
 274       sol.push_back (u + v);
 275     }
 276
 277   /* resubstitute */
 278   Real sub = 1.0 / 3 * A;
 279
 280   for (vsize i = sol.size (); i--;)
 281     {
 282       sol[i] -= sub;
 283
 284 #ifdef PARANOID
 285       assert (fabs (eval (sol[i])) < 1e-8);
 286 #endif
 287     }
 288
 289   return sol;
 290 }
 291
 292 Real
 293 Polynomial::lc () const
 294 {
 295   return coefs_.back ();
 296 }
 297
 298 Real &
 299 Polynomial::lc ()
 300 {
 301   return coefs_.back ();
 302 }
 303
 304 int
 305 Polynomial::degree ()const
 306 {
 307   return coefs_.size () -1;
 308 }
 309 /*
 310   all roots of quadratic eqn.
 311 */
 312 vector<Real>
 313 Polynomial::solve_quadric ()const
 314 {
 315   vector<Real> sol;
 316   /* normal form: x^2 + px + q = 0 */
 317   Real p = coefs_[1] / (2 * coefs_[2]);
 318   Real q = coefs_[0] / coefs_[2];
 319
 320   Real D = p * p - q;
 321
 322   if (D > 0)
 323     {
 324       D = sqrt (D);
 325
 326       sol.push_back (D - p);
 327       sol.push_back (-D - p);
 328     }
 329   return sol;
 330 }
 331
 332 /* solve linear equation */
 333 vector<Real>
 334 Polynomial::solve_linear ()const
 335 {
 336   vector<Real> s;
 337   if (coefs_[1])
 338     s.push_back (-coefs_[0] / coefs_[1]);
 339   return s;
 340 }
 341
 342 vector<Real>
 343 Polynomial::solve () const
 344 {
 345   Polynomial *me = (Polynomial *) this;
 346   me->clean ();
 347
 348   switch (degree ())
 349     {
 350     case 1:
 351       return solve_linear ();
 352     case 2:
 353       return solve_quadric ();
 354     case 3:
 355       return solve_cubic ();
 356     }
 357   vector<Real> s;
 358   return s;
 359 }
 360
 361 void
 362 Polynomial::operator *= (Polynomial const &p2)
 363 {
 364   *this = multiply (*this, p2);
 365 }
 366