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