2 This file is part of LilyPond, the GNU music typesetter.
4 Copyright (C) 1993--2015 Han-Wen Nienhuys <hanwen@xs4all.nl>
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.
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.
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/>.
20 #include "polynomial.hh"
30 Een beter milieu begint bij uzelf. Hergebruik!
33 This was ripped from Rayce, a raytracer I once wrote.
37 Polynomial::eval (Real x) const
42 for (vsize i = coefs_.size (); i--;)
43 p = x * p + coefs_[i];
49 Polynomial::multiply (const Polynomial &p1, const Polynomial &p2)
53 ssize_t deg = p1.degree () + p2.degree ();
54 for (ssize_t i = 0; i <= deg; i++)
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];
66 Polynomial::minmax (Real l, Real r, bool ret_max) const
71 programming_error ("left bound greater than right bound for polynomial minmax. flipping bounds.");
77 sols.push_back (eval (l));
78 sols.push_back (eval (r));
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> ());
88 return ret_max ? sols.back () : sols[0];
92 Polynomial::differentiate ()
94 for (int i = 1; i <= degree (); i++)
95 coefs_[i - 1] = coefs_[i] * i;
100 Polynomial::power (int exponent, const Polynomial &src)
103 Polynomial dest (1), base (src);
106 classic int power. invariant: src^exponent = dest * src ^ e
107 greetings go out to Lex Bijlsma & Jaap vd Woude */
112 dest = multiply (dest, base);
118 base = multiply (base, base);
125 static Real const FUDGE = 1e-8;
131 We only do relative comparisons. Absolute comparisons break down in
134 && (fabs (coefs_.back ()) < FUDGE * fabs (back (coefs_, 1))
140 Polynomial::operator += (Polynomial const &p)
142 while (degree () < p.degree ())
143 coefs_.push_back (0.0);
145 for (int i = 0; i <= p.degree (); i++)
146 coefs_[i] += p.coefs_[i];
150 Polynomial::operator -= (Polynomial const &p)
152 while (degree () < p.degree ())
153 coefs_.push_back (0.0);
155 for (int i = 0; i <= p.degree (); i++)
156 coefs_[i] -= p.coefs_[i];
160 Polynomial::scalarmultiply (Real fact)
162 for (int i = 0; i <= degree (); i++)
167 Polynomial::set_negate (const Polynomial &src)
169 for (int i = 0; i <= src.degree (); i++)
170 coefs_[i] = -src.coefs_[i];
175 Polynomial::set_mod (const Polynomial &u, const Polynomial &v)
181 for (ssize_t k = u.degree () - v.degree () - 1; k >= 0; k -= 2)
182 coefs_[k] = -coefs_[k];
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];
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];
196 ssize_t k = v.degree () - 1;
197 while (k >= 0 && coefs_[k] == 0.0)
200 coefs_.resize (1 + ((k < 0) ? 0 : k));
205 Polynomial::check_sol (Real x) const
208 Polynomial p (*this);
212 if (abs (f) > abs (d) * FUDGE)
213 programming_error ("not a root of polynomial\n");
217 Polynomial::check_sols (vector<Real> roots) const
219 for (vsize i = 0; i < roots.size (); i++)
220 check_sol (roots[i]);
223 Polynomial::Polynomial (Real a, Real b)
225 coefs_.push_back (a);
227 coefs_.push_back (b);
231 inline Real cubic_root (Real x)
234 return pow (x, 1.0 / 3.0);
236 return -pow (-x, 1.0 / 3.0);
247 Polynomial::solve_cubic ()const
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];
257 * substitute x = y - A/3 to eliminate quadric term: x^3 +px + q = 0
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);
264 /* use Cardano's formula */
271 if (iszero (q)) /* one triple solution */
277 else /* one single and one double solution */
279 Real u = cubic_root (-q);
281 sol.push_back (2 * u);
287 /* Casus irreducibilis: three real solutions */
288 Real phi = 1.0 / 3 * acos (-q / sqrt (-cb));
289 Real t = 2 * sqrt (-p);
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));
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);
302 sol.push_back (u + v);
306 Real sub = 1.0 / 3 * A;
308 for (vsize i = sol.size (); i--;)
313 assert (fabs (eval (sol[i])) < 1e-8);
321 Polynomial::lc () const
323 return coefs_.back ();
329 return coefs_.back ();
333 Polynomial::degree ()const
335 return coefs_.size () - 1;
338 all roots of quadratic eqn.
341 Polynomial::solve_quadric ()const
344 /* normal form: x^2 + px + q = 0 */
345 Real p = coefs_[1] / (2 * coefs_[2]);
346 Real q = coefs_[0] / coefs_[2];
354 sol.push_back (D - p);
355 sol.push_back (-D - p);
360 /* solve linear equation */
362 Polynomial::solve_linear ()const
366 s.push_back (-coefs_[0] / coefs_[1]);
371 Polynomial::solve () const
373 Polynomial *me = (Polynomial *) this;
379 return solve_linear ();
381 return solve_quadric ();
383 return solve_cubic ();
390 Polynomial::operator *= (Polynomial const &p2)
392 *this = multiply (*this, p2);