Re: revised quadratic.fs

On Aug 6, 12:03 am, Krishna Myneni <krishna.myn...@xxxxxxxxxxx> wrote:
It's surprising that something as simple as code for solving the
quadratic equation would have to be revisited, but the straightforward
approach is prone to problems with accuracy for certain parameter
sets. The book, Numerical Recipes, by Press, et al.,  has a discussion
about this.

I have revised my implementation, with the aim of getting it ready for
submission to the FSL. The latest version may be found at

ftp://ccreweb.org/software/fsl/extras/quadratic.fs

For comparison, the previous version may be found in

ftp://ccreweb.org/software/fsl/obsolete/

Any problem reports, suggestions for improvements, or suggestions for
other test cases and/ or test results would be welcome.

Krishna


Thanks to all for the feedback on this code. I have made another round
of significant revisions to quadratic.fs, to factor the code in a
sensible way. The changes largely make use of the factoring proposed
by humptydumpty, with the exception of his scaling transformation,
which I feel is unnecessary. Version 1.2 of quadratic.fs also makes
use of a couple of complex number words from the FSL complex library
(the library itself need not be loaded) to further clarify the
calculation. There is an overall performance penalty for the
factoring, but it is small, at least on an optimized Forth system such
as gforth-fast.

V1.2 of quadratic.fs is located at

ftp://ccreweb.org/software/fsl/extras/

It is also reproduced below.

Krishna

--

\ quadratic.fs
\
\ Solve for the two complex solutions of the quadratic equation:
\
\ a*x^2 + b*x + c = 0
\
\ for real coefficients a, b, and c. The two complex roots are
\ z1 = x1r + i*x1i and z2 = x2r + i*x2i
\
\ Copyright 2002--2009, K. Myneni, krishna.myneni@xxxxxxxxxxx
\ Permission is granted to modify and use this code
\ for any application, provided this notice is preserved.
\
\ Notes:
\
\ 1. Ensure that argument "a" is not zero, or an infinity will
result;
\ the correct solution of the simple linear equation will not be
\ given.
\
\ 2. The complex number roots, returned by solve_quadratic, are
\ ordered as follows:
\
\ a) For real roots, r1 and r2, with r1 < r2, the return
\ order is ( F: -- z2 z1 ), where z1 = r1 +i0, z2 = r2 + i0
\
\ b) For complex conjugate roots, z1 = e - i*f, z2 = e + i*f,
\ (f > 0), the ordering on the stack is ( F: -- z2 z1 ).
\
\ Revisions:
\
\ 2002-11-02 km; first version.
\ 2003-10-25 Christopher Brannon; fixed problem with calculation
\ of complex roots.
\ 2007-11-04 km; revised comments; added test code; save and
restore base.
\ 2009-08-05 km; revised to preserve accuracy when the product a*c
is
\ much less than b^2; see [1]. Added new test case.
\ 2009-08-11 km; added test case for purely imaginary roots.
\ 2009-08-18 km; factored the code with modification of example
posted
\ to comp.lang.forth by "humptydumpty", and used
words from
\ the jvn-dnw complex library; revised comments;
V1.2
\ References:
\
\ 1. W.H. Press, et. al., Numerical Recipes in C, 2nd ed., pp.
183--184,
\ eqns. 5.6.4 and 5.6.5.

CR .( QUADRATIC V1.2 18 August 2009 KM )
BASE @ DECIMAL

\ Useful words from the complex library
[UNDEFINED] zdup [IF] : zdup fover fover ; [THEN]
[UNDEFINED] zconj [IF] : conjg fnegate ; [THEN]

fvariable qa
fvariable 2qa
fvariable qb
fvariable qc

: q_discriminant ( F: -- d )
qb f@ fdup f* 4e qa f@ f* qc f@ f* f- ;

: q_real_roots ( F: d -- r1 r2 )
fsqrt qb f@ fdup f0< IF f- ELSE f+ fnegate THEN 2e f/
qc f@ fover f/ fswap qa f@ f/ ;

: q_complex_roots ( F: d -- z1 z2 )
fabs fsqrt 2qa f@ f/ \ imaginary part
qb f@ fnegate 2qa f@ f/ \ real part
fswap zdup conjg \ complex conjugate
;

: solve_quadratic ( F: a b c -- z1 z2 )
qc f! qb f! fdup qa f! 2e f* 2qa f!
q_discriminant
fdup f0< IF \ complex conjugate roots
q_complex_roots
ELSE \ two real roots
q_real_roots \ F: -- r1 r2
0e fswap 0e \ promote to z1 z2
THEN
;

BASE !

TEST-CODE? [IF] \ test code
==============================================
[undefined] T{ [IF] include ttester.fs [THEN]
BASE @ DECIMAL

1e-15 rel-near F!
1e-256 abs-near F!
set-near

\ Examples from:
\
\ http://www.purplemath.com/modules/quadform.htm

-2e 3e F/ fconstant -2/3
-3e 2e F/ fconstant -3/2
5e fsqrt fconstant sqrt{5}
2e fsqrt 3e F/ fconstant sqrt{2}/3
3e fsqrt 2e F/ fconstant sqrt{3}/2
10e fsqrt 2e F/ fconstant sqrt{10}/2

: zz}t rrrr}t ;

CR
TESTING SOLVE_QUADRATIC
t{ 1e 3e -4e solve_quadratic -> 1e 0e
-4e 0e zz}t
t{ 2e -4e -3e solve_quadratic -> 1e sqrt{10}/2 F- 0e 1e sqrt{10}/2
F+ 0e zz}t
t{ 1e -2e -4e solve_quadratic -> 1e sqrt{5} F- 0e 1e sqrt{5} F
+ 0e zz}t
t{ 9e 12e 4e solve_quadratic -> -2/3 0e
-2/3 0e zz}t
t{ 3e 4e 2e solve_quadratic -> -2/3 sqrt{2}/3 -2/3 sqrt{2}/3
fnegate zz}t
t{ 1e 3e 3e solve_quadratic -> -3/2 sqrt{3}/2 -3/2 sqrt{3}/2
fnegate zz}t
t{ 1e 0e 1e solve_quadratic -> 0e 1e
0e -1e zz}t

\ Test case which loses accuracy with ordinary quadratic formula:
\
\ x^2 + x + c = 0
\
\ when c << 1, the approximate solution is x1 = -c, x2 = -1 + c
\
t{ 1e 1e 1e-17 solve_quadratic -> -1e-17 0e -1e 1e-17 f
+ 0e zz}t


BASE !
[THEN]
.

Relevant Pages

  • Re: A dead subject
    ... > more like Spock, control your emotions. ... roots with the equation in any number of forms. ... all use the standard form exclusively and we don't always ...
    (sci.math)
  • Re: Ummm... Whats a "-0" in Spinlandia?
    ... it can only have 2 roots. ... You deny that quadratic have 2 roots? ... and yes this is a quadratic equation. ... "it aint what you know that gets you in trouble ...
    (talk.origins)
  • Re: revised quadratic.fs
    ... With respect to the ordering of the roots, ... For the real roots, I'm inclined to just add the extra code ... for the following reason: A quadratic equation can ... check, if he/she is extra paranoid (and when it comes to numerics, one ...
    (comp.lang.forth)
  • Help. Where is my error?
    ... program below on quadratic equation and will be most grateful if someone could ... quadratic equation has two distinct real roots: ... quadratic equation has two distinct imaginary roots. ... printf ("the values of the imaginary roots ...
    (comp.lang.c)
Loading