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