Hi all,
For several years now, the quadratic formula has come up in discussions of
stack shuffling.
All of the implementations that have come up, however, implemented the naive
formula:
x =(*-b +/- sqrt(b^2-4ac))/(2a)
The problem here is that if a and c are small, you're going to run into
accuracy problems. A better implementation is as follows, where x_1 and x_2
are the roots:
q = -1/2(b + sgn(b) * sqrt(b^2 - 4ac))
x_1 = q/a
x_2 = c/q
This is from Numerical Recipes in C++, page 227.
For fun, can anyone come up with a cleaner implementation of this?
: discriminant ( c a b -- d )
sq -rot * 4 * - sqrt ; inline
: q ( c a b -- q )
[ discriminant ] [ sgn ] [ ] tri
[ * ] dip + -0.5 * ; inline
: roots ( c a q -- x1 x2 )
tuck [ / ] [ swap / ] 2bi* ; inline
: quadratic ( c b a -- x1 x2 )
swap [ drop ] [ q ] 3bi roots ;
For comparison, here is my original version:
: monic ( c b a -- c' b' ) tuck [ / ] 2bi@ ;
: discriminant ( c b -- b d ) tuck sq 4 / swap - sqrt ;
: critical ( b d -- -b/2 d ) [ -2 / ] dip ;
: +- ( x y -- x+y x-y ) [ + ] [ - ] 2bi ;
: quadratic ( c b a -- alpha beta ) monic discriminant critical +- ;
-------------------------------------------------------------------------
This SF.Net email is sponsored by the Moblin Your Move Developer's challenge
Build the coolest Linux based applications with Moblin SDK & win great prizes
Grand prize is a trip for two to an Open Source event anywhere in the world
http://moblin-contest.org/redirect.php?banner_id=100&url=/
_______________________________________________
Factor-talk mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/factor-talk
