Issue 5167/6: Changes: show \markup xxx = ... \etc assignments

[lilypond.git] / flower / polynomial.cc

diff --git a/flower/polynomial.cc b/flower/polynomial.cc
index 56deefb714f739b111a150df40b69bd9b541e899..ec8605058c2bcdb0ce4c1ea5889d1463a882f434 100644 (file)
--- 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--2008 Han-Wen Nienhuys <hanwen@xs4all.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"
@@ -10,7 +23,6 @@
 
 #include <cmath>
 
-
 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<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 ()
 {
@@ -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.