2 poly.cc -- routines for manipulation of polynomials in one var
4 (c) 1993--2006 Han-Wen Nienhuys <hanwen@xs4all.nl>
7 #include "polynomial.hh"
13 Een beter milieu begint bij uzelf. Hergebruik!
16 This was ripped from Rayce, a raytracer I once wrote.
20 Polynomial::eval (Real x
) const
25 for (vsize i
= coefs_
.size (); i
--;)
26 p
= x
* p
+ coefs_
[i
];
32 Polynomial::multiply (const Polynomial
&p1
, const Polynomial
&p2
)
36 int deg
= p1
.degree () + p2
.degree ();
37 for (int i
= 0; i
<= deg
; i
++)
39 dest
.coefs_
.push_back (0);
40 for (int j
= 0; j
<= i
; j
++)
41 if (i
- j
<= p2
.degree () && j
<= p1
.degree ())
42 dest
.coefs_
.back () += p1
.coefs_
[j
] * p2
.coefs_
[i
- j
];
49 Polynomial::differentiate ()
51 for (int i
= 1; i
<= degree (); i
++)
52 coefs_
[i
- 1] = coefs_
[i
] * i
;
57 Polynomial::power (int exponent
, const Polynomial
&src
)
60 Polynomial
dest (1), base (src
);
63 classic int power. invariant: src^exponent = dest * src ^ e
64 greetings go out to Lex Bijlsma & Jaap vd Woude */
69 dest
= multiply (dest
, base
);
75 base
= multiply (base
, base
);
82 static Real
const FUDGE
= 1e-8;
88 We only do relative comparisons. Absolute comparisons break down in
91 && (fabs (coefs_
.back ()) < FUDGE
* fabs (back (coefs_
, 1))
97 Polynomial::operator += (Polynomial
const &p
)
99 while (degree () < p
.degree ())
100 coefs_
.push_back (0.0);
102 for (int i
= 0; i
<= p
.degree (); i
++)
103 coefs_
[i
] += p
.coefs_
[i
];
107 Polynomial::operator -= (Polynomial
const &p
)
109 while (degree () < p
.degree ())
110 coefs_
.push_back (0.0);
112 for (int i
= 0; i
<= p
.degree (); i
++)
113 coefs_
[i
] -= p
.coefs_
[i
];
117 Polynomial::scalarmultiply (Real fact
)
119 for (int i
= 0; i
<= degree (); i
++)
124 Polynomial::set_negate (const Polynomial
&src
)
126 for (int i
= 0; i
<= src
.degree (); i
++)
127 coefs_
[i
] = -src
.coefs_
[i
];
132 Polynomial::set_mod (const Polynomial
&u
, const Polynomial
&v
)
138 for (int k
= u
.degree () - v
.degree () - 1; k
>= 0; k
-= 2)
139 coefs_
[k
] = -coefs_
[k
];
141 for (int k
= u
.degree () - v
.degree (); k
>= 0; k
--)
142 for (int j
= v
.degree () + k
- 1; j
>= k
; j
--)
143 coefs_
[j
] = -coefs_
[j
] - coefs_
[v
.degree () + k
] * v
.coefs_
[j
- k
];
148 for (int k
= u
.degree () - v
.degree (); k
>= 0; k
--)
149 for (int j
= v
.degree () + k
- 1; j
>= k
; j
--)
150 coefs_
[j
] -= coefs_
[v
.degree () + k
] * v
.coefs_
[j
- k
];
153 int k
= v
.degree () - 1;
154 while (k
>= 0 && coefs_
[k
] == 0.0)
157 coefs_
.resize (1+ ((k
< 0) ? 0 : k
));
162 Polynomial::check_sol (Real x
) const
165 Polynomial
p (*this);
169 if (abs (f
) > abs (d
) * FUDGE
)
172 warning ("x=%f is not a root of polynomial\n"
173 "f (x)=%f, f' (x)=%f \n", x, f, d); */
177 Polynomial::check_sols (std::vector
<Real
> roots
) const
179 for (vsize i
= 0; i
< roots
.size (); i
++)
180 check_sol (roots
[i
]);
183 Polynomial::Polynomial (Real a
, Real b
)
185 coefs_
.push_back (a
);
187 coefs_
.push_back (b
);
191 inline Real
cubic_root (Real x
)
194 return pow (x
, 1.0 / 3.0);
196 return -pow (-x
, 1.0 / 3.0);
207 Polynomial::solve_cubic ()const
209 std::vector
<Real
> sol
;
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];
217 * substitute x = y - A/3 to eliminate quadric term: x^3 +px + q = 0
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
);
224 /* use Cardano's formula */
231 if (iszero (q
)) { /* one triple solution */
236 else { /* one single and one double solution */
237 Real u
= cubic_root (-q
);
239 sol
.push_back (2 * u
);
245 /* Casus irreducibilis: three real solutions */
246 Real phi
= 1.0 / 3 * acos (-q
/ sqrt (-cb
));
247 Real t
= 2 * sqrt (-p
);
249 sol
.push_back (t
* cos (phi
));
250 sol
.push_back (-t
* cos (phi
+ M_PI
/ 3));
251 sol
.push_back (-t
* cos (phi
- M_PI
/ 3));
255 /* one real solution */
256 Real sqrt_D
= sqrt (D
);
257 Real u
= cubic_root (sqrt_D
- q
);
258 Real v
= -cubic_root (sqrt_D
+ q
);
260 sol
.push_back (u
+ v
);
264 Real sub
= 1.0 / 3 * A
;
266 for (vsize i
= sol
.size (); i
--;)
271 assert (fabs (eval (sol
[i
])) < 1e-8);
279 Polynomial::lc () const
281 return coefs_
.back ();
287 return coefs_
.back ();
291 Polynomial::degree ()const
293 return coefs_
.size () -1;
296 all roots of quadratic eqn.
299 Polynomial::solve_quadric ()const
301 std::vector
<Real
> sol
;
302 /* normal form: x^2 + px + q = 0 */
303 Real p
= coefs_
[1] / (2 * coefs_
[2]);
304 Real q
= coefs_
[0] / coefs_
[2];
312 sol
.push_back (D
- p
);
313 sol
.push_back (-D
- p
);
318 /* solve linear equation */
320 Polynomial::solve_linear ()const
324 s
.push_back (-coefs_
[0] / coefs_
[1]);
329 Polynomial::solve () const
331 Polynomial
*me
= (Polynomial
*) this;
337 return solve_linear ();
339 return solve_quadric ();
341 return solve_cubic ();
348 Polynomial::operator *= (Polynomial
const &p2
)
350 *this = multiply (*this, p2
);