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