X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=flower%2Fpolynomial.cc;h=56deefb714f739b111a150df40b69bd9b541e899;hb=5b4b0d6e9a197e8f9eb085b7c2ad78b8be3e5cfc;hp=7b60880e7d901d40f4779fb4c47f9cc2fda72062;hpb=bdf4ab13203502e7ec7cf9cf5896527643a07c1f;p=lilypond.git diff --git a/flower/polynomial.cc b/flower/polynomial.cc index 7b60880e7d..56deefb714 100644 --- a/flower/polynomial.cc +++ b/flower/polynomial.cc @@ -1,13 +1,18 @@ /* poly.cc -- routines for manipulation of polynomials in one var - (c) 1993--2005 Han-Wen Nienhuys + (c) 1993--2008 Han-Wen Nienhuys */ #include "polynomial.hh" +#include "warn.hh" + #include + +using namespace std; + /* Een beter milieu begint bij uzelf. Hergebruik! @@ -16,12 +21,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 +40,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; @@ -48,10 +53,8 @@ void Polynomial::differentiate () { for (int i = 1; i <= degree (); i++) - { - coefs_[i - 1] = coefs_[i] * i; - } - coefs_.pop (); + coefs_[i - 1] = coefs_[i] * i; + coefs_.pop_back (); } Polynomial @@ -89,16 +92,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]; @@ -108,7 +111,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]; @@ -155,7 +158,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 (); } @@ -168,24 +171,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. */ @@ -195,8 +195,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 @@ -205,10 +204,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]; @@ -231,37 +230,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; @@ -276,13 +279,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 @@ -293,10 +296,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]; @@ -307,23 +310,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; @@ -338,7 +341,7 @@ Polynomial::solve () const case 3: return solve_cubic (); } - Array s; + vector s; return s; }