]> git.donarmstrong.com Git - lilypond.git/blob - flower/polynomial.cc
Merge branch 'master' of git://git.sv.gnu.org/lilypond
[lilypond.git] / flower / polynomial.cc
1 /*
2   poly.cc -- routines for manipulation of polynomials in one var
3
4   (c) 1993--2007 Han-Wen Nienhuys <hanwen@xs4all.nl>
5 */
6
7 #include "polynomial.hh"
8
9 #include "warn.hh"
10
11 #include <cmath>
12
13
14 using namespace std;
15
16 /*
17   Een beter milieu begint bij uzelf. Hergebruik!
18
19
20   This was ripped from Rayce, a raytracer I once wrote.
21 */
22
23 Real
24 Polynomial::eval (Real x) const
25 {
26   Real p = 0.0;
27
28   // horner's scheme
29   for (vsize i = coefs_.size (); i--;)
30     p = x * p + coefs_[i];
31
32   return p;
33 }
34
35 Polynomial
36 Polynomial::multiply (const Polynomial &p1, const Polynomial &p2)
37 {
38   Polynomial dest;
39
40   int deg = p1.degree () + p2.degree ();
41   for (int i = 0; i <= deg; i++)
42     {
43       dest.coefs_.push_back (0);
44       for (int j = 0; j <= i; j++)
45         if (i - j <= p2.degree () && j <= p1.degree ())
46           dest.coefs_.back () += p1.coefs_[j] * p2.coefs_[i - j];
47     }
48
49   return dest;
50 }
51
52 void
53 Polynomial::differentiate ()
54 {
55   for (int i = 1; i <= degree (); i++)
56     coefs_[i - 1] = coefs_[i] * i;
57   coefs_.pop_back ();
58 }
59
60 Polynomial
61 Polynomial::power (int exponent, const Polynomial &src)
62 {
63   int e = exponent;
64   Polynomial dest (1), base (src);
65
66   /*
67     classic int power. invariant: src^exponent = dest * src ^ e
68     greetings go out to Lex Bijlsma & Jaap vd Woude */
69   while (e > 0)
70     {
71       if (e % 2)
72         {
73           dest = multiply (dest, base);
74           e--;
75         }
76       else
77
78         {
79           base = multiply (base, base);
80           e /= 2;
81         }
82     }
83   return dest;
84 }
85
86 static Real const FUDGE = 1e-8;
87
88 void
89 Polynomial::clean ()
90 {
91   /*
92     We only do relative comparisons. Absolute comparisons break down in
93     degenerate cases.  */
94   while (degree () > 0
95          && (fabs (coefs_.back ()) < FUDGE * fabs (back (coefs_, 1))
96              || !coefs_.back ()))
97     coefs_.pop_back ();
98 }
99
100 void
101 Polynomial::operator += (Polynomial const &p)
102 {
103   while (degree () < p.degree ())
104     coefs_.push_back (0.0);
105
106   for (int i = 0; i <= p.degree (); i++)
107     coefs_[i] += p.coefs_[i];
108 }
109
110 void
111 Polynomial::operator -= (Polynomial const &p)
112 {
113   while (degree () < p.degree ())
114     coefs_.push_back (0.0);
115
116   for (int i = 0; i <= p.degree (); i++)
117     coefs_[i] -= p.coefs_[i];
118 }
119
120 void
121 Polynomial::scalarmultiply (Real fact)
122 {
123   for (int i = 0; i <= degree (); i++)
124     coefs_[i] *= fact;
125 }
126
127 void
128 Polynomial::set_negate (const Polynomial &src)
129 {
130   for (int i = 0; i <= src.degree (); i++)
131     coefs_[i] = -src.coefs_[i];
132 }
133
134 /// mod of #u/v#
135 int
136 Polynomial::set_mod (const Polynomial &u, const Polynomial &v)
137 {
138   (*this) = u;
139
140   if (v.lc () < 0.0)
141     {
142       for (int k = u.degree () - v.degree () - 1; k >= 0; k -= 2)
143         coefs_[k] = -coefs_[k];
144
145       for (int k = u.degree () - v.degree (); k >= 0; k--)
146         for (int j = v.degree () + k - 1; j >= k; j--)
147           coefs_[j] = -coefs_[j] - coefs_[v.degree () + k] * v.coefs_[j - k];
148     }
149   else
150
151     {
152       for (int k = u.degree () - v.degree (); k >= 0; k--)
153         for (int j = v.degree () + k - 1; j >= k; j--)
154           coefs_[j] -= coefs_[v.degree () + k] * v.coefs_[j - k];
155     }
156
157   int k = v.degree () - 1;
158   while (k >= 0 && coefs_[k] == 0.0)
159     k--;
160
161   coefs_.resize (1+ ((k < 0) ? 0 : k));
162   return degree ();
163 }
164
165 void
166 Polynomial::check_sol (Real x) const
167 {
168   Real f = eval (x);
169   Polynomial p (*this);
170   p.differentiate ();
171   Real d = p.eval (x);
172
173   if (abs (f) > abs (d) * FUDGE)
174     programming_error ("not a root of polynomial\n");
175 }
176
177 void
178 Polynomial::check_sols (vector<Real> roots) const
179 {
180   for (vsize i = 0; i < roots.size (); i++)
181     check_sol (roots[i]);
182 }
183
184 Polynomial::Polynomial (Real a, Real b)
185 {
186   coefs_.push_back (a);
187   if (b)
188     coefs_.push_back (b);
189 }
190
191 /* cubic root. */
192 inline Real cubic_root (Real x)
193 {
194   if (x > 0.0)
195     return pow (x, 1.0 / 3.0);
196   else if (x < 0.0)
197     return -pow (-x, 1.0 / 3.0);
198   return 0.0;
199 }
200
201 static bool
202 iszero (Real r)
203 {
204   return !r;
205 }
206
207 vector<Real>
208 Polynomial::solve_cubic ()const
209 {
210   vector<Real> sol;
211
212   /* normal form: x^3 + Ax^2 + Bx + C = 0 */
213   Real A = coefs_[2] / coefs_[3];
214   Real B = coefs_[1] / coefs_[3];
215   Real C = coefs_[0] / coefs_[3];
216
217   /*
218    * substitute x = y - A/3 to eliminate quadric term: x^3 +px + q = 0
219    */
220
221   Real sq_A = A *A;
222   Real p = 1.0 / 3 * (-1.0 / 3 * sq_A + B);
223   Real q = 1.0 / 2 * (2.0 / 27 * A *sq_A - 1.0 / 3 * A *B + C);
224
225   /* use Cardano's formula */
226
227   Real cb = p * p * p;
228   Real D = q * q + cb;
229
230   if (iszero (D))
231     {
232       if (iszero (q)) { /* one triple solution */
233         sol.push_back (0);
234         sol.push_back (0);
235         sol.push_back (0);
236       }
237       else {            /* one single and one double solution */
238         Real u = cubic_root (-q);
239
240         sol.push_back (2 * u);
241         sol.push_back (-u);
242       }
243     }
244   else if (D < 0)
245     {
246     /* Casus irreducibilis: three real solutions */
247       Real phi = 1.0 / 3 * acos (-q / sqrt (-cb));
248       Real t = 2 * sqrt (-p);
249
250       sol.push_back (t * cos (phi));
251       sol.push_back (-t * cos (phi + M_PI / 3));
252       sol.push_back (-t * cos (phi - M_PI / 3));
253     }
254   else
255     {
256       /* one real solution */
257       Real sqrt_D = sqrt (D);
258       Real u = cubic_root (sqrt_D - q);
259       Real v = -cubic_root (sqrt_D + q);
260
261       sol.push_back (u + v);
262     }
263
264   /* resubstitute */
265   Real sub = 1.0 / 3 * A;
266
267   for (vsize i = sol.size (); i--;)
268     {
269       sol[i] -= sub;
270
271 #ifdef PARANOID
272       assert (fabs (eval (sol[i])) < 1e-8);
273 #endif
274     }
275
276   return sol;
277 }
278
279 Real
280 Polynomial::lc () const
281 {
282   return coefs_.back ();
283 }
284
285 Real &
286 Polynomial::lc ()
287 {
288   return coefs_.back ();
289 }
290
291 int
292 Polynomial::degree ()const
293 {
294   return coefs_.size () -1;
295 }
296 /*
297   all roots of quadratic eqn.
298 */
299 vector<Real>
300 Polynomial::solve_quadric ()const
301 {
302   vector<Real> sol;
303   /* normal form: x^2 + px + q = 0 */
304   Real p = coefs_[1] / (2 * coefs_[2]);
305   Real q = coefs_[0] / coefs_[2];
306
307   Real D = p * p - q;
308
309   if (D > 0)
310     {
311       D = sqrt (D);
312
313       sol.push_back (D - p);
314       sol.push_back (-D - p);
315     }
316   return sol;
317 }
318
319 /* solve linear equation */
320 vector<Real>
321 Polynomial::solve_linear ()const
322 {
323   vector<Real> s;
324   if (coefs_[1])
325     s.push_back (-coefs_[0] / coefs_[1]);
326   return s;
327 }
328
329 vector<Real>
330 Polynomial::solve () const
331 {
332   Polynomial *me = (Polynomial *) this;
333   me->clean ();
334
335   switch (degree ())
336     {
337     case 1:
338       return solve_linear ();
339     case 2:
340       return solve_quadric ();
341     case 3:
342       return solve_cubic ();
343     }
344   vector<Real> s;
345   return s;
346 }
347
348 void
349 Polynomial::operator *= (Polynomial const &p2)
350 {
351   *this = multiply (*this, p2);
352 }
353