* Another grand 2003 update.

[lilypond.git] / flower / polynomial.cc

diff --git a/flower/polynomial.cc b/flower/polynomial.cc
index df8847fe04ebb52f599d2f689d705de19a1e04cd..b9f93092d5996605b2c47edb665cc1723d8f48c9 100644 (file)
--- a/flower/polynomial.cc
+++ b/flower/polynomial.cc
@@ -1,7 +1,7 @@
 /*
    poly.cc -- routines for manipulation of polynomials in one var
 
-  (c) 1993--2000 Han-Wen Nienhuys <hanwen@cs.uu.nl>
+  (c)  1993--2003 Han-Wen Nienhuys <hanwen@cs.uu.nl>
  */
 
 #include <math.h>
@@ -94,23 +94,29 @@ Polynomial::clean ()
 }
 
 
-Polynomial 
-Polynomial::add (const Polynomial & p1, const Polynomial & p2)
+
+void
+Polynomial::operator += (Polynomial const &p)
 {
-  Polynomial dest;
-  int tempord =  p2.degree () >? p1.degree ();
-  for (int i = 0; i <= tempord; i++)
-    {
-      Real temp = 0.0;
-      if (i <= p1.degree ())
-       temp += p1.coefs_[i];
-      if (i <= p2.degree ())
-       temp += p2.coefs_[i];
-      dest.coefs_.push (temp);
-    }
-  return dest;
+  while (degree () < p.degree())
+    coefs_.push (0.0);
+
+  for (int i = 0; i <= p.degree(); i++)
+    coefs_[i] += p.coefs_[i];
 }
 
+
+void
+Polynomial::operator -= (Polynomial const &p)
+{
+  while (degree () < p.degree())
+    coefs_.push (0.0);
+
+  for (int i = 0; i <= p.degree(); i++)
+    coefs_[i] -= p.coefs_[i];
+}
+
+
 void
 Polynomial::scalarmultiply (Real fact)
 {
@@ -118,24 +124,7 @@ Polynomial::scalarmultiply (Real fact)
     coefs_[i] *= fact;
 }
 
-Polynomial
-Polynomial::subtract (const Polynomial & p1, const Polynomial & p2)
-{
-  Polynomial dest;
-  int tempord =  p2.degree () >? p1.degree ();
-  
-  for (int i = 0; i <= tempord; i++)
-    {
-      Real temp = 0.0; // can't store result directly.. a=a-b
-      if (i <= p1.degree ())
-       temp += p1.coefs_[i];
-      if (i <= p2.degree ())
-       temp -= p2.coefs_[i];
-      dest.coefs_.push (temp);
-    }
-  return dest;
-  
-}
+
 
 void
 Polynomial::set_negate (const Polynomial & src)
@@ -237,8 +226,8 @@ Polynomial::solve_cubic ()const
 
   /* use Cardano's formula */
 
-  Real cb_p = p * p * p;
-  Real D = q * q + cb_p;
+  Real cb = p * p * p;
+  Real D = q * q + cb;
 
   if (iszero (D)) {
     if (iszero (q)) {  /* one triple solution */
@@ -252,7 +241,7 @@ Polynomial::solve_cubic ()const
       sol.push (-u);
     }
   } else if (D < 0) {          /* Casus irreducibilis: three real solutions */
-    Real phi = 1.0 / 3 * acos (-q / sqrt (-cb_p));
+    Real phi = 1.0 / 3 * acos (-q / sqrt (-cb));
     Real t = 2 * sqrt (-p);
 
     sol.push (t * cos (phi));
@@ -356,14 +345,4 @@ Polynomial:: operator *= (Polynomial const &p2)
   *this = multiply (*this,p2);
 }  
 
-void
-Polynomial::operator += (Polynomial const &p)
-{
-  *this = add ( *this, p);
-}
 
-void
-Polynomial::operator -= (Polynomial const &p)
-{
-  *this = subtract (*this, p);
-}