/*
- 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 <hanwen@cs.uu.nl>
+ Copyright (C) 1993--2015 Han-Wen Nienhuys <hanwen@xs4all.nl>
+
+ 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 <http://www.gnu.org/licenses/>.
*/
#include "polynomial.hh"
+#include "warn.hh"
+
#include <cmath>
+
using namespace std;
/*
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;
{
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 (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_.push_back (0);
+ 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<Real> 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<Real> 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<Real> ());
+
+ return ret_max ? sols.back () : sols[0];
+}
+
void
Polynomial::differentiate ()
{
for (int i = 1; i <= degree (); i++)
coefs_[i - 1] = coefs_[i] * i;
- coefs_.pop ();
+ coefs_.pop_back ();
}
Polynomial
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;
}
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];
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];
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_.set_size (1+ ((k < 0) ? 0 : k));
+ coefs_.resize (1 + ((k < 0) ? 0 : k));
return degree ();
}
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<Real> roots) const
+Polynomial::check_sols (vector<Real> 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. */
return !r;
}
-Array<Real>
+vector<Real>
Polynomial::solve_cubic ()const
{
- Array<Real> sol;
+ vector<Real> sol;
/* normal form: x^3 + Ax^2 + Bx + C = 0 */
Real A = coefs_[2] / coefs_[3];
* 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 */
if (iszero (D))
{
- if (iszero (q)) { /* one triple solution */
- sol.push (0);
- sol.push (0);
- sol.push (0);
- }
- else { /* one single and one double solution */
- Real u = cubic_root (-q);
-
- sol.push (2 * u);
- sol.push (-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);
- sol.push (t * cos (phi));
- sol.push (-t * cos (phi + M_PI / 3));
- sol.push (-t * cos (phi - M_PI / 3));
+ 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
{
Real u = cubic_root (sqrt_D - q);
Real v = -cubic_root (sqrt_D + q);
- sol.push (u + v);
+ 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;
Real
Polynomial::lc () const
{
- return coefs_.top ();
+ return coefs_.back ();
}
Real &
Polynomial::lc ()
{
- return coefs_.top ();
+ return coefs_.back ();
}
-int
+ssize_t
Polynomial::degree ()const
{
- return coefs_.size () -1;
+ return coefs_.size () - 1;
}
/*
all roots of quadratic eqn.
*/
-Array<Real>
+vector<Real>
Polynomial::solve_quadric ()const
{
- Array<Real> sol;
+ vector<Real> sol;
/* normal form: x^2 + px + q = 0 */
Real p = coefs_[1] / (2 * coefs_[2]);
Real q = coefs_[0] / coefs_[2];
{
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<Real>
+vector<Real>
Polynomial::solve_linear ()const
{
- Array<Real> s;
+ vector<Real> s;
if (coefs_[1])
- s.push (-coefs_[0] / coefs_[1]);
+ s.push_back (-coefs_[0] / coefs_[1]);
return s;
}
-Array<Real>
+vector<Real>
Polynomial::solve () const
{
Polynomial *me = (Polynomial *) this;
case 3:
return solve_cubic ();
}
- Array<Real> s;
+ vector<Real> s;
return s;
}
projects / lilypond.git / blobdiff