2 poly.cc -- routines for manipulation of polynomials in one var
4 (c) 1993--2004 Han-Wen Nienhuys <hanwen@cs.uu.nl>
9 #include "polynomial.hh"
12 Een beter milieu begint bij uzelf. Hergebruik!
15 This was ripped from Rayce, a raytracer I once wrote.
19 Polynomial::eval (Real x
)const
24 for (int i
= coefs_
.size (); i
--; )
25 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
++)
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
];
49 Polynomial::differentiate ()
51 for (int i
= 1; i
<= degree (); i
++)
53 coefs_
[i
-1] = coefs_
[i
] * i
;
59 Polynomial::power (int exponent
, const Polynomial
& src
)
62 Polynomial
dest (1), base (src
);
65 classic int power. invariant: src^exponent = dest * src ^ e
66 greetings go out to Lex Bijlsma & Jaap vd Woude */
71 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
90 while (degree () > 0 &&
91 (fabs (coefs_
.top ()) < FUDGE
* fabs (coefs_
.top (1))
99 Polynomial::operator += (Polynomial
const &p
)
101 while (degree () < p
.degree())
104 for (int i
= 0; i
<= p
.degree(); i
++)
105 coefs_
[i
] += p
.coefs_
[i
];
110 Polynomial::operator -= (Polynomial
const &p
)
112 while (degree () < p
.degree())
115 for (int i
= 0; i
<= p
.degree(); i
++)
116 coefs_
[i
] -= p
.coefs_
[i
];
121 Polynomial::scalarmultiply (Real fact
)
123 for (int i
= 0; i
<= degree (); i
++)
130 Polynomial::set_negate (const Polynomial
& src
)
132 for (int i
= 0; i
<= src
.degree (); i
++)
133 coefs_
[i
] = -src
.coefs_
[i
];
138 Polynomial::set_mod (const Polynomial
&u
, const Polynomial
&v
)
143 for (int k
= u
.degree () - v
.degree () - 1; k
>= 0; k
-= 2)
144 coefs_
[k
] = -coefs_
[k
];
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
];
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
];
155 int k
= v
.degree () - 1;
156 while (k
>= 0 && coefs_
[k
] == 0.0)
159 coefs_
.set_size (1+ ((k
< 0) ? 0 : k
));
164 Polynomial::check_sol (Real x
) const
167 Polynomial
p (*this);
171 if ( abs (f
) > abs (d
) * FUDGE
)
174 warning ("x=%f is not a root of polynomial\n"
175 "f (x)=%f, f' (x)=%f \n", x, f, d); */
179 Polynomial::check_sols (Array
<Real
> roots
) const
181 for (int i
=0; i
< roots
.size (); i
++)
182 check_sol (roots
[i
]);
185 Polynomial::Polynomial (Real a
, Real b
)
193 inline Real
cubic_root (Real x
)
196 return pow (x
, 1.0/3.0) ;
198 return -pow (-x
, 1.0/3.0);
210 Polynomial::solve_cubic ()const
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];
220 * substitute x = y - A/3 to eliminate quadric term: x^3 +px + q = 0
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
);
227 /* use Cardano's formula */
233 if (iszero (q
)) { /* one triple solution */
237 } else { /* one single and one double solution */
238 Real u
= cubic_root (-q
);
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
);
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
);
259 Real sub
= 1.0 / 3 * A
;
261 for (int i
= sol
.size (); i
--;)
266 assert (fabs (eval (sol
[i
]) ) < 1e-8);
274 Polynomial::lc () const
276 return coefs_
.top ();
282 return coefs_
.top ();
286 Polynomial::degree ()const
288 return coefs_
.size () -1;
291 all roots of quadratic eqn.
294 Polynomial::solve_quadric ()const
297 /* normal form: x^2 + px + q = 0 */
298 Real p
= coefs_
[1] / (2 * coefs_
[2]);
299 Real q
= coefs_
[0] / coefs_
[2];
312 /* solve linear equation */
314 Polynomial::solve_linear ()const
318 s
.push ( -coefs_
[0] / coefs_
[1]);
324 Polynomial::solve () const
326 Polynomial
* me
= (Polynomial
*) this;
332 return solve_linear ();
334 return solve_quadric ();
336 return solve_cubic ();
343 Polynomial:: operator *= (Polynomial
const &p2
)
345 *this = multiply (*this,p2
);