X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=flower%2Fpolynomial.cc;h=ec8605058c2bcdb0ce4c1ea5889d1463a882f434;hb=b872748c6aa8bb721ced458691b38ac2fac5dfc8;hp=d9784c277069e363866077b3302c1d5b7887a425;hpb=f9214bac21e9926dc3248416f58190c98c4167a9;p=lilypond.git diff --git a/flower/polynomial.cc b/flower/polynomial.cc index d9784c2770..ec8605058c 100644 --- a/flower/polynomial.cc +++ b/flower/polynomial.cc @@ -1,7 +1,20 @@ /* - poly.cc -- routines for manipulation of polynomials in one var + This file is part of LilyPond, the GNU music typesetter. - (c) 1993--2007 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" @@ -10,7 +23,6 @@ #include - using namespace std; /* @@ -37,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 () { @@ -69,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; } @@ -92,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 (); } @@ -139,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 (); } @@ -218,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 */ @@ -229,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); @@ -288,10 +328,10 @@ Polynomial::lc () return coefs_.back (); } -int +ssize_t Polynomial::degree ()const { - return coefs_.size () -1; + return coefs_.size () - 1; } /* all roots of quadratic eqn.