flower/polynomial.cc

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