flower/polynomial.cc

   1 /*
   2   This file is part of LilyPond, the GNU music typesetter.
   3
   4   Copyright (C) 1993--2011 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 using namespace std;
  27
  28 /*
  29   Een beter milieu begint bij uzelf. Hergebruik!
  30
  31
  32   This was ripped from Rayce, a raytracer I once wrote.
  33 */
  34
  35 Real
  36 Polynomial::eval (Real x) const
  37 {
  38   Real p = 0.0;
  39
  40   // horner's scheme
  41   for (vsize i = coefs_.size (); i--;)
  42     p = x * p + coefs_[i];
  43
  44   return p;
  45 }
  46
  47 Polynomial
  48 Polynomial::multiply (const Polynomial &p1, const Polynomial &p2)
  49 {
  50   Polynomial dest;
  51
  52   int deg = p1.degree () + p2.degree ();
  53   for (int i = 0; i <= deg; i++)
  54     {
  55       dest.coefs_.push_back (0);
  56       for (int j = 0; j <= i; j++)
  57         if (i - j <= p2.degree () && j <= p1.degree ())
  58           dest.coefs_.back () += p1.coefs_[j] * p2.coefs_[i - j];
  59     }
  60
  61   return dest;
  62 }
  63
  64 void
  65 Polynomial::differentiate ()
  66 {
  67   for (int i = 1; i <= degree (); i++)
  68     coefs_[i - 1] = coefs_[i] * i;
  69   coefs_.pop_back ();
  70 }
  71
  72 Polynomial
  73 Polynomial::power (int exponent, const Polynomial &src)
  74 {
  75   int e = exponent;
  76   Polynomial dest (1), base (src);
  77
  78   /*
  79     classic int power. invariant: src^exponent = dest * src ^ e
  80     greetings go out to Lex Bijlsma & Jaap vd Woude */
  81   while (e > 0)
  82     {
  83       if (e % 2)
  84         {
  85           dest = multiply (dest, base);
  86           e--;
  87         }
  88       else
  89
  90         {
  91           base = multiply (base, base);
  92           e /= 2;
  93         }
  94     }
  95   return dest;
  96 }
  97
  98 static Real const FUDGE = 1e-8;
  99
 100 void
 101 Polynomial::clean ()
 102 {
 103   /*
 104     We only do relative comparisons. Absolute comparisons break down in
 105     degenerate cases.  */
 106   while (degree () > 0
 107          && (fabs (coefs_.back ()) < FUDGE * fabs (back (coefs_, 1))
 108              || !coefs_.back ()))
 109     coefs_.pop_back ();
 110 }
 111
 112 void
 113 Polynomial::operator += (Polynomial const &p)
 114 {
 115   while (degree () < p.degree ())
 116     coefs_.push_back (0.0);
 117
 118   for (int i = 0; i <= p.degree (); i++)
 119     coefs_[i] += p.coefs_[i];
 120 }
 121
 122 void
 123 Polynomial::operator -= (Polynomial const &p)
 124 {
 125   while (degree () < p.degree ())
 126     coefs_.push_back (0.0);
 127
 128   for (int i = 0; i <= p.degree (); i++)
 129     coefs_[i] -= p.coefs_[i];
 130 }
 131
 132 void
 133 Polynomial::scalarmultiply (Real fact)
 134 {
 135   for (int i = 0; i <= degree (); i++)
 136     coefs_[i] *= fact;
 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     {
 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     }
 161   else
 162
 163     {
 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_.resize (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     programming_error ("not a root of polynomial\n");
 187 }
 188
 189 void
 190 Polynomial::check_sols (vector<Real> roots) const
 191 {
 192   for (vsize i = 0; i < roots.size (); i++)
 193     check_sol (roots[i]);
 194 }
 195
 196 Polynomial::Polynomial (Real a, Real b)
 197 {
 198   coefs_.push_back (a);
 199   if (b)
 200     coefs_.push_back (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   return 0.0;
 211 }
 212
 213 static bool
 214 iszero (Real r)
 215 {
 216   return !r;
 217 }
 218
 219 vector<Real>
 220 Polynomial::solve_cubic ()const
 221 {
 222   vector<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;
 240   Real D = q * q + cb;
 241
 242   if (iszero (D))
 243     {
 244       if (iszero (q))   /* one triple solution */
 245         {
 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         {
 252           Real u = cubic_root (-q);
 253
 254           sol.push_back (2 * u);
 255           sol.push_back (-u);
 256         }
 257     }
 258   else if (D < 0)
 259     {
 260       /* Casus irreducibilis: three real solutions */
 261       Real phi = 1.0 / 3 * acos (-q / sqrt (-cb));
 262       Real t = 2 * sqrt (-p);
 263
 264       sol.push_back (t * cos (phi));
 265       sol.push_back (-t * cos (phi + M_PI / 3));
 266       sol.push_back (-t * cos (phi - M_PI / 3));
 267     }
 268   else
 269     {
 270       /* one real solution */
 271       Real sqrt_D = sqrt (D);
 272       Real u = cubic_root (sqrt_D - q);
 273       Real v = -cubic_root (sqrt_D + q);
 274
 275       sol.push_back (u + v);
 276     }
 277
 278   /* resubstitute */
 279   Real sub = 1.0 / 3 * A;
 280
 281   for (vsize i = sol.size (); i--;)
 282     {
 283       sol[i] -= sub;
 284
 285 #ifdef PARANOID
 286       assert (fabs (eval (sol[i])) < 1e-8);
 287 #endif
 288     }
 289
 290   return sol;
 291 }
 292
 293 Real
 294 Polynomial::lc () const
 295 {
 296   return coefs_.back ();
 297 }
 298
 299 Real &
 300 Polynomial::lc ()
 301 {
 302   return coefs_.back ();
 303 }
 304
 305 int
 306 Polynomial::degree ()const
 307 {
 308   return coefs_.size () - 1;
 309 }
 310 /*
 311   all roots of quadratic eqn.
 312 */
 313 vector<Real>
 314 Polynomial::solve_quadric ()const
 315 {
 316   vector<Real> sol;
 317   /* normal form: x^2 + px + q = 0 */
 318   Real p = coefs_[1] / (2 * coefs_[2]);
 319   Real q = coefs_[0] / coefs_[2];
 320
 321   Real D = p * p - q;
 322
 323   if (D > 0)
 324     {
 325       D = sqrt (D);
 326
 327       sol.push_back (D - p);
 328       sol.push_back (-D - p);
 329     }
 330   return sol;
 331 }
 332
 333 /* solve linear equation */
 334 vector<Real>
 335 Polynomial::solve_linear ()const
 336 {
 337   vector<Real> s;
 338   if (coefs_[1])
 339     s.push_back (-coefs_[0] / coefs_[1]);
 340   return s;
 341 }
 342
 343 vector<Real>
 344 Polynomial::solve () const
 345 {
 346   Polynomial *me = (Polynomial *) this;
 347   me->clean ();
 348
 349   switch (degree ())
 350     {
 351     case 1:
 352       return solve_linear ();
 353     case 2:
 354       return solve_quadric ();
 355     case 3:
 356       return solve_cubic ();
 357     }
 358   vector<Real> s;
 359   return s;
 360 }
 361
 362 void
 363 Polynomial::operator *= (Polynomial const &p2)
 364 {
 365   *this = multiply (*this, p2);
 366 }
 367