X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=flower%2Fpolynomial.cc;h=ec8605058c2bcdb0ce4c1ea5889d1463a882f434;hb=90e4d7057f3857da049dfda3d130017d4719bd6b;hp=5f719147674860a3cca8d4035b4282e80660625e;hpb=ff90d34b47b91a7321ff9b21fe20133911376364;p=lilypond.git diff --git a/flower/polynomial.cc b/flower/polynomial.cc index 5f71914767..ec8605058c 100644 --- a/flower/polynomial.cc +++ b/flower/polynomial.cc @@ -1,12 +1,28 @@ /* - poly.cc -- routines for manipulation of polynomials in one var + This file is part of LilyPond, the GNU music typesetter. - (c) 1993--2006 Han-Wen Nienhuys + Copyright (C) 1993--2015 Han-Wen Nienhuys + + LilyPond is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + LilyPond is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with LilyPond. If not, see . */ #include "polynomial.hh" +#include "warn.hh" + #include + using namespace std; /* @@ -33,18 +49,44 @@ Polynomial::multiply (const Polynomial &p1, const Polynomial &p2) { Polynomial dest; - int deg = p1.degree () + p2.degree (); - for (int i = 0; i <= deg; i++) + ssize_t deg = p1.degree () + p2.degree (); + for (ssize_t i = 0; i <= deg; i++) { dest.coefs_.push_back (0); - for (int j = 0; j <= i; j++) - if (i - j <= p2.degree () && j <= p1.degree ()) - dest.coefs_.back () += p1.coefs_[j] * p2.coefs_[i - j]; + for (ssize_t j = 0; j <= i; j++) + if (i - j <= p2.degree () && j <= p1.degree ()) + dest.coefs_.back () += p1.coefs_[j] * p2.coefs_[i - j]; } return dest; } +Real +Polynomial::minmax (Real l, Real r, bool ret_max) const +{ + vector sols; + if (l > r) + { + programming_error ("left bound greater than right bound for polynomial minmax. flipping bounds."); + l = l + r; + r = l - r; + l = l - r; + } + + sols.push_back (eval (l)); + sols.push_back (eval (r)); + + Polynomial deriv (*this); + deriv.differentiate (); + vector maxmins = deriv.solve (); + for (vsize i = 0; i < maxmins.size (); i++) + if (maxmins[i] >= l && maxmins[i] <= r) + sols.push_back (eval (maxmins[i])); + vector_sort (sols, less ()); + + return ret_max ? sols.back () : sols[0]; +} + void Polynomial::differentiate () { @@ -65,16 +107,16 @@ Polynomial::power (int exponent, const Polynomial &src) while (e > 0) { if (e % 2) - { - dest = multiply (dest, base); - e--; - } + { + dest = multiply (dest, base); + e--; + } else - { - base = multiply (base, base); - e /= 2; - } + { + base = multiply (base, base); + e /= 2; + } } return dest; } @@ -88,8 +130,8 @@ Polynomial::clean () We only do relative comparisons. Absolute comparisons break down in degenerate cases. */ while (degree () > 0 - && (fabs (coefs_.back ()) < FUDGE * fabs (back (coefs_, 1)) - || !coefs_.back ())) + && (fabs (coefs_.back ()) < FUDGE * fabs (back (coefs_, 1)) + || !coefs_.back ())) coefs_.pop_back (); } @@ -135,26 +177,26 @@ Polynomial::set_mod (const Polynomial &u, const Polynomial &v) if (v.lc () < 0.0) { - for (int k = u.degree () - v.degree () - 1; k >= 0; k -= 2) - coefs_[k] = -coefs_[k]; + for (ssize_t k = u.degree () - v.degree () - 1; k >= 0; k -= 2) + coefs_[k] = -coefs_[k]; - for (int k = u.degree () - v.degree (); k >= 0; k--) - for (int j = v.degree () + k - 1; j >= k; j--) - coefs_[j] = -coefs_[j] - coefs_[v.degree () + k] * v.coefs_[j - k]; + for (ssize_t k = u.degree () - v.degree (); k >= 0; k--) + for (ssize_t j = v.degree () + k - 1; j >= k; j--) + coefs_[j] = -coefs_[j] - coefs_[v.degree () + k] * v.coefs_[j - k]; } else { - for (int k = u.degree () - v.degree (); k >= 0; k--) - for (int j = v.degree () + k - 1; j >= k; j--) - coefs_[j] -= coefs_[v.degree () + k] * v.coefs_[j - k]; + for (ssize_t k = u.degree () - v.degree (); k >= 0; k--) + for (ssize_t j = v.degree () + k - 1; j >= k; j--) + coefs_[j] -= coefs_[v.degree () + k] * v.coefs_[j - k]; } - int k = v.degree () - 1; + ssize_t k = v.degree () - 1; while (k >= 0 && coefs_[k] == 0.0) k--; - coefs_.resize (1+ ((k < 0) ? 0 : k)); + coefs_.resize (1 + ((k < 0) ? 0 : k)); return degree (); } @@ -167,14 +209,11 @@ Polynomial::check_sol (Real x) const Real d = p.eval (x); if (abs (f) > abs (d) * FUDGE) - ; - /* - warning ("x=%f is not a root of polynomial\n" - "f (x)=%f, f' (x)=%f \n", x, f, d); */ + programming_error ("not a root of polynomial\n"); } void -Polynomial::check_sols (std::vector roots) const +Polynomial::check_sols (vector roots) const { for (vsize i = 0; i < roots.size (); i++) check_sol (roots[i]); @@ -203,10 +242,10 @@ iszero (Real r) return !r; } -std::vector +vector Polynomial::solve_cubic ()const { - std::vector sol; + vector sol; /* normal form: x^3 + Ax^2 + Bx + C = 0 */ Real A = coefs_[2] / coefs_[3]; @@ -217,9 +256,9 @@ Polynomial::solve_cubic ()const * substitute x = y - A/3 to eliminate quadric term: x^3 +px + q = 0 */ - Real sq_A = A *A; + Real sq_A = A * A; Real p = 1.0 / 3 * (-1.0 / 3 * sq_A + B); - Real q = 1.0 / 2 * (2.0 / 27 * A *sq_A - 1.0 / 3 * A *B + C); + Real q = 1.0 / 2 * (2.0 / 27 * A * sq_A - 1.0 / 3 * A * B + C); /* use Cardano's formula */ @@ -228,21 +267,23 @@ Polynomial::solve_cubic ()const if (iszero (D)) { - if (iszero (q)) { /* one triple solution */ - sol.push_back (0); - sol.push_back (0); - sol.push_back (0); - } - else { /* one single and one double solution */ - Real u = cubic_root (-q); - - sol.push_back (2 * u); - sol.push_back (-u); - } + if (iszero (q)) /* one triple solution */ + { + sol.push_back (0); + sol.push_back (0); + sol.push_back (0); + } + else /* one single and one double solution */ + { + Real u = cubic_root (-q); + + sol.push_back (2 * u); + sol.push_back (-u); + } } else if (D < 0) { - /* Casus irreducibilis: three real solutions */ + /* Casus irreducibilis: three real solutions */ Real phi = 1.0 / 3 * acos (-q / sqrt (-cb)); Real t = 2 * sqrt (-p); @@ -287,18 +328,18 @@ Polynomial::lc () return coefs_.back (); } -int +ssize_t Polynomial::degree ()const { - return coefs_.size () -1; + return coefs_.size () - 1; } /* all roots of quadratic eqn. */ -std::vector +vector Polynomial::solve_quadric ()const { - std::vector sol; + vector sol; /* normal form: x^2 + px + q = 0 */ Real p = coefs_[1] / (2 * coefs_[2]); Real q = coefs_[0] / coefs_[2]; @@ -316,16 +357,16 @@ Polynomial::solve_quadric ()const } /* solve linear equation */ -std::vector +vector Polynomial::solve_linear ()const { - std::vector s; + vector s; if (coefs_[1]) s.push_back (-coefs_[0] / coefs_[1]); return s; } -std::vector +vector Polynomial::solve () const { Polynomial *me = (Polynomial *) this; @@ -340,7 +381,7 @@ Polynomial::solve () const case 3: return solve_cubic (); } - std::vector s; + vector s; return s; }