2 poly.cc -- routines for manipulation of polynomials in one var
4 (c) 1993--2000 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))
98 Polynomial::add (const Polynomial
& p1
, const Polynomial
& p2
)
101 int tempord
= p2
.degree () >? p1
.degree ();
102 for (int i
= 0; i
<= tempord
; i
++)
105 if (i
<= p1
.degree ())
106 temp
+= p1
.coefs_
[i
];
107 if (i
<= p2
.degree ())
108 temp
+= p2
.coefs_
[i
];
109 dest
.coefs_
.push (temp
);
115 Polynomial::scalarmultiply (Real fact
)
117 for (int i
= 0; i
<= degree (); i
++)
122 Polynomial::subtract (const Polynomial
& p1
, const Polynomial
& p2
)
125 int tempord
= p2
.degree () >? p1
.degree ();
127 for (int i
= 0; i
<= tempord
; i
++)
129 Real temp
= 0.0; // can't store result directly.. a=a-b
130 if (i
<= p1
.degree ())
131 temp
+= p1
.coefs_
[i
];
132 if (i
<= p2
.degree ())
133 temp
-= p2
.coefs_
[i
];
134 dest
.coefs_
.push (temp
);
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
)
154 for (int k
= u
.degree () - v
.degree () - 1; k
>= 0; k
-= 2)
155 coefs_
[k
] = -coefs_
[k
];
157 for (int k
= u
.degree () - v
.degree (); k
>= 0; k
--)
158 for (int j
= v
.degree () + k
- 1; j
>= k
; j
--)
159 coefs_
[j
] = -coefs_
[j
] - coefs_
[v
.degree () + k
] * v
.coefs_
[j
- k
];
161 for (int k
= u
.degree () - v
.degree (); k
>= 0; k
--)
162 for (int j
= v
.degree () + k
- 1; j
>= k
; j
--)
163 coefs_
[j
] -= coefs_
[v
.degree () + k
] * v
.coefs_
[j
- k
];
166 int k
= v
.degree () - 1;
167 while (k
>= 0 && coefs_
[k
] == 0.0)
170 coefs_
.set_size (1+ ((k
< 0) ? 0 : k
));
175 Polynomial::check_sol (Real x
) const
178 Polynomial
p (*this);
182 if ( abs (f
) > abs (d
) * FUDGE
)
185 warning ("x=%f is not a root of polynomial\n"
186 "f (x)=%f, f' (x)=%f \n", x, f, d); */
190 Polynomial::check_sols (Array
<Real
> roots
) const
192 for (int i
=0; i
< roots
.size (); i
++)
193 check_sol (roots
[i
]);
196 Polynomial::Polynomial (Real a
, Real b
)
204 inline Real
cubic_root (Real x
)
207 return pow (x
, 1.0/3.0) ;
209 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 */
240 Real cb_p
= p
* p
* p
;
241 Real D
= q
* q
+ cb_p
;
244 if (iszero (q
)) { /* one triple solution */
248 } else { /* one single and one double solution */
249 Real u
= cubic_root (-q
);
254 } else if (D
< 0) { /* Casus irreducibilis: three real solutions */
255 Real phi
= 1.0 / 3 * acos (-q
/ sqrt (-cb_p
));
256 Real t
= 2 * sqrt (-p
);
258 sol
.push (t
* cos (phi
));
259 sol
.push (-t
* cos (phi
+ M_PI
/ 3));
260 sol
.push ( -t
* cos (phi
- M_PI
/ 3));
261 } else { /* one real solution */
262 Real sqrt_D
= sqrt (D
);
263 Real u
= cubic_root (sqrt_D
- q
);
264 Real v
= -cubic_root (sqrt_D
+ q
);
270 Real sub
= 1.0 / 3 * A
;
272 for (int i
= sol
.size (); i
--;)
277 assert (fabs (eval (sol
[i
]) ) < 1e-8);
285 Polynomial::lc () const
287 return coefs_
.top ();
293 return coefs_
.top ();
297 Polynomial::degree ()const
299 return coefs_
.size () -1;
302 all roots of quadratic eqn.
305 Polynomial::solve_quadric ()const
308 /* normal form: x^2 + px + q = 0 */
309 Real p
= coefs_
[1] / (2 * coefs_
[2]);
310 Real q
= coefs_
[0] / coefs_
[2];
323 /* solve linear equation */
325 Polynomial::solve_linear ()const
329 s
.push ( -coefs_
[0] / coefs_
[1]);
335 Polynomial::solve () const
337 Polynomial
* me
= (Polynomial
*) this;
343 return solve_linear ();
345 return solve_quadric ();
347 return solve_cubic ();
354 Polynomial:: operator *= (Polynomial
const &p2
)
356 *this = multiply (*this,p2
);
360 Polynomial::operator += (Polynomial
const &p
)
362 *this = add ( *this, p
);
366 Polynomial::operator -= (Polynomial
const &p
)
368 *this = subtract (*this, p
);