2 This file is part of LilyPond, the GNU music typesetter.
4 Copyright (C) 1993--2011 Han-Wen Nienhuys <hanwen@xs4all.nl>
6 LilyPond is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 3 of the License, or
9 (at your option) any later version.
11 LilyPond is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
16 You should have received a copy of the GNU General Public License
17 along with LilyPond. If not, see <http://www.gnu.org/licenses/>.
20 #include "polynomial.hh"
30 Een beter milieu begint bij uzelf. Hergebruik!
33 This was ripped from Rayce, a raytracer I once wrote.
37 Polynomial::eval (Real x
) const
42 for (vsize i
= coefs_
.size (); i
--;)
43 p
= x
* p
+ coefs_
[i
];
49 Polynomial::multiply (const Polynomial
&p1
, const Polynomial
&p2
)
53 int deg
= p1
.degree () + p2
.degree ();
54 for (int i
= 0; i
<= deg
; i
++)
56 dest
.coefs_
.push_back (0);
57 for (int j
= 0; j
<= i
; j
++)
58 if (i
- j
<= p2
.degree () && j
<= p1
.degree ())
59 dest
.coefs_
.back () += p1
.coefs_
[j
] * p2
.coefs_
[i
- j
];
66 Polynomial::differentiate ()
68 for (int i
= 1; i
<= degree (); i
++)
69 coefs_
[i
- 1] = coefs_
[i
] * i
;
74 Polynomial::power (int exponent
, const Polynomial
&src
)
77 Polynomial
dest (1), base (src
);
80 classic int power. invariant: src^exponent = dest * src ^ e
81 greetings go out to Lex Bijlsma & Jaap vd Woude */
86 dest
= multiply (dest
, base
);
92 base
= multiply (base
, base
);
99 static Real
const FUDGE
= 1e-8;
105 We only do relative comparisons. Absolute comparisons break down in
108 && (fabs (coefs_
.back ()) < FUDGE
* fabs (back (coefs_
, 1))
114 Polynomial::operator += (Polynomial
const &p
)
116 while (degree () < p
.degree ())
117 coefs_
.push_back (0.0);
119 for (int i
= 0; i
<= p
.degree (); i
++)
120 coefs_
[i
] += p
.coefs_
[i
];
124 Polynomial::operator -= (Polynomial
const &p
)
126 while (degree () < p
.degree ())
127 coefs_
.push_back (0.0);
129 for (int i
= 0; i
<= p
.degree (); i
++)
130 coefs_
[i
] -= p
.coefs_
[i
];
134 Polynomial::scalarmultiply (Real fact
)
136 for (int i
= 0; i
<= degree (); i
++)
141 Polynomial::set_negate (const Polynomial
&src
)
143 for (int i
= 0; i
<= src
.degree (); i
++)
144 coefs_
[i
] = -src
.coefs_
[i
];
149 Polynomial::set_mod (const Polynomial
&u
, const Polynomial
&v
)
155 for (int k
= u
.degree () - v
.degree () - 1; k
>= 0; k
-= 2)
156 coefs_
[k
] = -coefs_
[k
];
158 for (int k
= u
.degree () - v
.degree (); k
>= 0; k
--)
159 for (int j
= v
.degree () + k
- 1; j
>= k
; j
--)
160 coefs_
[j
] = -coefs_
[j
] - coefs_
[v
.degree () + k
] * v
.coefs_
[j
- k
];
165 for (int k
= u
.degree () - v
.degree (); k
>= 0; k
--)
166 for (int j
= v
.degree () + k
- 1; j
>= k
; j
--)
167 coefs_
[j
] -= coefs_
[v
.degree () + k
] * v
.coefs_
[j
- k
];
170 int k
= v
.degree () - 1;
171 while (k
>= 0 && coefs_
[k
] == 0.0)
174 coefs_
.resize (1+ ((k
< 0) ? 0 : k
));
179 Polynomial::check_sol (Real x
) const
182 Polynomial
p (*this);
186 if (abs (f
) > abs (d
) * FUDGE
)
187 programming_error ("not a root of polynomial\n");
191 Polynomial::check_sols (vector
<Real
> roots
) const
193 for (vsize i
= 0; i
< roots
.size (); i
++)
194 check_sol (roots
[i
]);
197 Polynomial::Polynomial (Real a
, Real b
)
199 coefs_
.push_back (a
);
201 coefs_
.push_back (b
);
205 inline Real
cubic_root (Real x
)
208 return pow (x
, 1.0 / 3.0);
210 return -pow (-x
, 1.0 / 3.0);
221 Polynomial::solve_cubic ()const
225 /* normal form: x^3 + Ax^2 + Bx + C = 0 */
226 Real A
= coefs_
[2] / coefs_
[3];
227 Real B
= coefs_
[1] / coefs_
[3];
228 Real C
= coefs_
[0] / coefs_
[3];
231 * substitute x = y - A/3 to eliminate quadric term: x^3 +px + q = 0
235 Real p
= 1.0 / 3 * (-1.0 / 3 * sq_A
+ B
);
236 Real q
= 1.0 / 2 * (2.0 / 27 * A
*sq_A
- 1.0 / 3 * A
*B
+ C
);
238 /* use Cardano's formula */
245 if (iszero (q
)) { /* one triple solution */
250 else { /* one single and one double solution */
251 Real u
= cubic_root (-q
);
253 sol
.push_back (2 * u
);
259 /* Casus irreducibilis: three real solutions */
260 Real phi
= 1.0 / 3 * acos (-q
/ sqrt (-cb
));
261 Real t
= 2 * sqrt (-p
);
263 sol
.push_back (t
* cos (phi
));
264 sol
.push_back (-t
* cos (phi
+ M_PI
/ 3));
265 sol
.push_back (-t
* cos (phi
- M_PI
/ 3));
269 /* one real solution */
270 Real sqrt_D
= sqrt (D
);
271 Real u
= cubic_root (sqrt_D
- q
);
272 Real v
= -cubic_root (sqrt_D
+ q
);
274 sol
.push_back (u
+ v
);
278 Real sub
= 1.0 / 3 * A
;
280 for (vsize i
= sol
.size (); i
--;)
285 assert (fabs (eval (sol
[i
])) < 1e-8);
293 Polynomial::lc () const
295 return coefs_
.back ();
301 return coefs_
.back ();
305 Polynomial::degree ()const
307 return coefs_
.size () -1;
310 all roots of quadratic eqn.
313 Polynomial::solve_quadric ()const
316 /* normal form: x^2 + px + q = 0 */
317 Real p
= coefs_
[1] / (2 * coefs_
[2]);
318 Real q
= coefs_
[0] / coefs_
[2];
326 sol
.push_back (D
- p
);
327 sol
.push_back (-D
- p
);
332 /* solve linear equation */
334 Polynomial::solve_linear ()const
338 s
.push_back (-coefs_
[0] / coefs_
[1]);
343 Polynomial::solve () const
345 Polynomial
*me
= (Polynomial
*) this;
351 return solve_linear ();
353 return solve_quadric ();
355 return solve_cubic ();
362 Polynomial::operator *= (Polynomial
const &p2
)
364 *this = multiply (*this, p2
);