X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=flower%2Fpolynomial.cc;h=1530baee83e6c7ca00052394e600f9bbf3d04aad;hb=3bd8002b58b0408bf9d0f796744fb012b40ba47d;hp=906c948c9f905c0c725a0e9a81e45d4bb7f6cf46;hpb=108cf0e8c08c8e15e2a800feb161cfad9057daa8;p=lilypond.git diff --git a/flower/polynomial.cc b/flower/polynomial.cc index 906c948c9f..1530baee83 100644 --- a/flower/polynomial.cc +++ b/flower/polynomial.cc @@ -1,13 +1,31 @@ /* - poly.cc -- routines for manipulation of polynomials in one var + This file is part of LilyPond, the GNU music typesetter. - (c) 1993--2005 Han-Wen Nienhuys + Copyright (C) 1993--2011 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; + /* Een beter milieu begint bij uzelf. Hergebruik! @@ -16,12 +34,12 @@ */ Real -Polynomial::eval (Real x)const +Polynomial::eval (Real x) const { Real p = 0.0; // horner's scheme - for (int i = coefs_.size (); i--;) + for (vsize i = coefs_.size (); i--;) p = x * p + coefs_[i]; return p; @@ -35,10 +53,10 @@ Polynomial::multiply (const Polynomial &p1, const Polynomial &p2) int deg = p1.degree () + p2.degree (); for (int i = 0; i <= deg; i++) { - dest.coefs_.push (0); + dest.coefs_.push_back (0); for (int j = 0; j <= i; j++) if (i - j <= p2.degree () && j <= p1.degree ()) - dest.coefs_.top () += p1.coefs_[j] * p2.coefs_[i - j]; + dest.coefs_.back () += p1.coefs_[j] * p2.coefs_[i - j]; } return dest; @@ -49,7 +67,7 @@ Polynomial::differentiate () { for (int i = 1; i <= degree (); i++) coefs_[i - 1] = coefs_[i] * i; - coefs_.pop (); + coefs_.pop_back (); } Polynomial @@ -87,16 +105,16 @@ Polynomial::clean () We only do relative comparisons. Absolute comparisons break down in degenerate cases. */ while (degree () > 0 - && (fabs (coefs_.top ()) < FUDGE *fabs (coefs_.top (1)) - || !coefs_.top ())) - coefs_.pop (); + && (fabs (coefs_.back ()) < FUDGE * fabs (back (coefs_, 1)) + || !coefs_.back ())) + coefs_.pop_back (); } void Polynomial::operator += (Polynomial const &p) { while (degree () < p.degree ()) - coefs_.push (0.0); + coefs_.push_back (0.0); for (int i = 0; i <= p.degree (); i++) coefs_[i] += p.coefs_[i]; @@ -106,7 +124,7 @@ void Polynomial::operator -= (Polynomial const &p) { while (degree () < p.degree ()) - coefs_.push (0.0); + coefs_.push_back (0.0); for (int i = 0; i <= p.degree (); i++) coefs_[i] -= p.coefs_[i]; @@ -153,7 +171,7 @@ Polynomial::set_mod (const Polynomial &u, const Polynomial &v) while (k >= 0 && coefs_[k] == 0.0) k--; - coefs_.set_size (1+ ((k < 0) ? 0 : k)); + coefs_.resize (1+ ((k < 0) ? 0 : k)); return degree (); } @@ -166,24 +184,21 @@ 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 (Array roots) const +Polynomial::check_sols (vector roots) const { - for (int i = 0; i < roots.size (); i++) + for (vsize i = 0; i < roots.size (); i++) check_sol (roots[i]); } Polynomial::Polynomial (Real a, Real b) { - coefs_.push (a); + coefs_.push_back (a); if (b) - coefs_.push (b); + coefs_.push_back (b); } /* cubic root. */ @@ -193,8 +208,7 @@ inline Real cubic_root (Real x) return pow (x, 1.0 / 3.0); else if (x < 0.0) return -pow (-x, 1.0 / 3.0); - else - return 0.0; + return 0.0; } static bool @@ -203,10 +217,10 @@ iszero (Real r) return !r; } -Array +vector Polynomial::solve_cubic ()const { - Array sol; + vector sol; /* normal form: x^3 + Ax^2 + Bx + C = 0 */ Real A = coefs_[2] / coefs_[3]; @@ -229,37 +243,41 @@ Polynomial::solve_cubic ()const if (iszero (D)) { if (iszero (q)) { /* one triple solution */ - sol.push (0); - sol.push (0); - sol.push (0); + 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 (2 * u); - sol.push (-u); + sol.push_back (2 * u); + sol.push_back (-u); } } - else if (D < 0) { /* Casus irreducibilis: three real solutions */ - Real phi = 1.0 / 3 * acos (-q / sqrt (-cb)); - Real t = 2 * sqrt (-p); - - sol.push (t * cos (phi)); - sol.push (-t * cos (phi + M_PI / 3)); - sol.push (-t * cos (phi - M_PI / 3)); - } - else { /* one real solution */ - Real sqrt_D = sqrt (D); - Real u = cubic_root (sqrt_D - q); - Real v = -cubic_root (sqrt_D + q); - - sol.push (u + v); - } + else if (D < 0) + { + /* Casus irreducibilis: three real solutions */ + Real phi = 1.0 / 3 * acos (-q / sqrt (-cb)); + Real t = 2 * sqrt (-p); + + sol.push_back (t * cos (phi)); + sol.push_back (-t * cos (phi + M_PI / 3)); + sol.push_back (-t * cos (phi - M_PI / 3)); + } + else + { + /* one real solution */ + Real sqrt_D = sqrt (D); + Real u = cubic_root (sqrt_D - q); + Real v = -cubic_root (sqrt_D + q); + + sol.push_back (u + v); + } /* resubstitute */ Real sub = 1.0 / 3 * A; - for (int i = sol.size (); i--;) + for (vsize i = sol.size (); i--;) { sol[i] -= sub; @@ -274,13 +292,13 @@ Polynomial::solve_cubic ()const Real Polynomial::lc () const { - return coefs_.top (); + return coefs_.back (); } Real & Polynomial::lc () { - return coefs_.top (); + return coefs_.back (); } int @@ -291,10 +309,10 @@ Polynomial::degree ()const /* all roots of quadratic eqn. */ -Array +vector Polynomial::solve_quadric ()const { - Array sol; + vector sol; /* normal form: x^2 + px + q = 0 */ Real p = coefs_[1] / (2 * coefs_[2]); Real q = coefs_[0] / coefs_[2]; @@ -305,23 +323,23 @@ Polynomial::solve_quadric ()const { D = sqrt (D); - sol.push (D - p); - sol.push (-D - p); + sol.push_back (D - p); + sol.push_back (-D - p); } return sol; } /* solve linear equation */ -Array +vector Polynomial::solve_linear ()const { - Array s; + vector s; if (coefs_[1]) - s.push (-coefs_[0] / coefs_[1]); + s.push_back (-coefs_[0] / coefs_[1]); return s; } -Array +vector Polynomial::solve () const { Polynomial *me = (Polynomial *) this; @@ -336,7 +354,7 @@ Polynomial::solve () const case 3: return solve_cubic (); } - Array s; + vector s; return s; }