Run `make grand-replace'.

[lilypond.git] / flower / polynomial.cc

diff --git a/flower/polynomial.cc b/flower/polynomial.cc
index 92447b63001046d6ac8cf4541eb4f380fd063388..56deefb714f739b111a150df40b69bd9b541e899 100644 (file)
--- 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 <hanwen@cs.uu.nl>
+  (c) 1993--2008 Han-Wen Nienhuys <hanwen@xs4all.nl>
 */
 
 #include "polynomial.hh"
 
+#include "warn.hh"
+
 #include <cmath>
 
+
+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;
@@ -32,13 +37,13 @@ Polynomial::multiply (const Polynomial &p1, const Polynomial &p2)
 {
   Polynomial dest;
 
-  int deg= p1.degree () + p2.degree ();
+  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;
@@ -47,11 +52,9 @@ Polynomial::multiply (const Polynomial &p1, const Polynomial &p2)
 void
 Polynomial::differentiate ()
 {
-  for (int i = 1; i<= degree (); i++)
-    {
-      coefs_[i - 1] = coefs_[i] * i;
-    }
-  coefs_.pop ();
+  for (int i = 1; i <= degree (); i++)
+    coefs_[i - 1] = coefs_[i] * i;
+  coefs_.pop_back ();
 }
 
 Polynomial
@@ -89,26 +92,26 @@ 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)
+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];
 }
 
 void
-Polynomial::operator-= (Polynomial const &p)
+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 ();
 }
 
@@ -167,24 +170,22 @@ Polynomial::check_sol (Real x) const
   p.differentiate ();
   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);        */
+  if (abs (f) > abs (d) * FUDGE)
+    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. */
@@ -194,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
@@ -204,10 +204,10 @@ iszero (Real r)
   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];
@@ -230,34 +230,41 @@ Polynomial::solve_cubic ()const
   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 */
+       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);
+    }
+  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_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 (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;
 
@@ -272,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
@@ -289,10 +296,10 @@ Polynomial::degree ()const
 /*
   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];
@@ -303,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<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;
@@ -334,12 +341,12 @@ Polynomial::solve () const
     case 3:
       return solve_cubic ();
     }
-  Array<Real> s;
+  vector<Real> s;
   return s;
 }
 
 void
-Polynomial:: operator*= (Polynomial const &p2)
+Polynomial::operator *= (Polynomial const &p2)
 {
   *this = multiply (*this, p2);
 }