Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]
Groups > comp.lang.forth > #9688 > unrolled thread
| Started by | mhx@iae.nl (Marcel Hendrix) |
|---|---|
| First post | 2012-02-24 02:52 +0200 |
| Last post | 2012-02-23 20:25 -0800 |
| Articles | 7 — 3 participants |
Back to article view | Back to comp.lang.forth
Ferrari's method to solve a quartic mhx@iae.nl (Marcel Hendrix) - 2012-02-24 02:52 +0200
Re: Ferrari's method to solve a quartic Paul Rubin <no.email@nospam.invalid> - 2012-02-23 18:39 -0800
Re: Ferrari's method to solve a quartic Krishna Myneni <krishna.myneni@ccreweb.org> - 2012-02-23 18:58 -0800
Re: Ferrari's method to solve a quartic mhx@iae.nl (Marcel Hendrix) - 2012-02-25 10:11 +0200
Re: Ferrari's method to solve a quartic Krishna Myneni <krishna.myneni@ccreweb.org> - 2012-02-28 18:59 -0800
Re: Ferrari's method to solve a quartic mhx@iae.nl (Marcel Hendrix) - 2012-02-25 15:32 +0200
Re: Ferrari's method to solve a quartic Krishna Myneni <krishna.myneni@ccreweb.org> - 2012-02-23 20:25 -0800
| From | mhx@iae.nl (Marcel Hendrix) |
|---|---|
| Date | 2012-02-24 02:52 +0200 |
| Subject | Ferrari's method to solve a quartic |
| Message-ID | <14830111018435@frunobulax.edu> |
Here is Ferrari's method to solve a quartic equation.
I only tested it with Cardano's problem, where the absolute
error of the roots is about 10^-14 to 10^-15. This high
error is caused by the use of complex numbers (implemented
with standard double precision in iForth).
-marcel
--
(*
* LANGUAGE : ANS Forth with extensions
* PROJECT : Forth Environments
* DESCRIPTION : Ferrari's method to solve a quartic equation
* CATEGORY : Numeric Utility
* AUTHOR : Marcel Hendrix
* LAST CHANGE : February 23, 2012, Marcel Hendrix
*)
NEEDS -miscutil
NEEDS -fcbrt
NEEDS -cplx_fsl
REVISION -ferrari "--- Solve a Quartic Version 0.01 ---"
PRIVATES
DOC
(*
The quartic is the highest order polynomial equation that can be solved by
radicals in the general case (i.e., one where the coefficients can take any
value).
Lodovico Ferrari is attributed with the discovery of the solution to the
quartic in 1540, but since this solution, like all algebraic solutions of
the quartic, requires the solution of a cubic to be found, it couldn't be
published immediately. The solution of the quartic was published together
with that of the cubic by Ferrari's mentor Gerolamo Cardano in the book Ars
Magna (1545).
The proof that four is the highest degree of a general polynomial for which
such solutions can be found was first given in the Abel-Ruffini theorem in
1824, proving that all attempts at solving the higher order polynomials
would be futile. The notes left by Evariste Galois prior to dying in a duel
in 1832 later led to an elegant complete theory of the roots of polynomials,
of which this theorem was one result.
( http://en.wikipedia.org/wiki/Quartic_function )
*)
ENDDOC
1e FVALUE A PRIVATE
0e FVALUE B PRIVATE
6e FVALUE C PRIVATE
-60e FVALUE D PRIVATE
36e FVALUE E PRIVATE
0e FVALUE Alpha PRIVATE
0e FVALUE Beta PRIVATE
0e FVALUE Gamma PRIVATE
0e FVALUE P PRIVATE
0e FVALUE Q PRIVATE
0e FVALUE +R PRIVATE
0+0i ZVALUE +U PRIVATE
0+0i ZVALUE y PRIVATE
0+0i ZVALUE W PRIVATE
CREATE qsol PRIVATE 4 ZFLOATS ALLOT
: Alpha! ( -- )
B A F/ FSQR -3e F* 8e F/
C A F/ F+ TO Alpha ; PRIVATE
: Beta! ( -- )
B A F/ FCUBE 8e F/
B C F* A FSQR F2* F/ F-
D A F/ F+ TO Beta ; PRIVATE
: Gamma! ( -- )
B A F/ FQUAD -3e F* 256e F/
C B FSQR F* A FCUBE 16e F* F/ F+
B D F* A FSQR 4e F* F/ F-
E A F/ F+ TO Gamma ; PRIVATE
: Beta=0? ( -- bool )
Beta F0<> IF FALSE EXIT ENDIF
Alpha FSQR Gamma 4e F* F- FSQRT FLOCAL a2
B A -4e F* F/ FLOCAL a1
Alpha FNEGATE a2 F+ F2/ FSQRT a1 F+ 0e R,I->Z qsol 0 COMPLEX[] Z!
Alpha FNEGATE a2 F- F2/ FSQRT a1 F+ 0e R,I->Z qsol 1 COMPLEX[] Z!
Alpha FNEGATE a2 F+ F2/ FSQRT a1 F- 0e R,I->Z qsol 2 COMPLEX[] Z!
Alpha FNEGATE a2 F- F2/ FSQRT a1 F- 0e R,I->Z qsol 3 COMPLEX[] Z!
TRUE ; PRIVATE
: P! ( -- ) Alpha FSQR -12e F/ Gamma F- TO P ; PRIVATE
: Q! ( -- ) Alpha FCUBE -108e F/ Alpha Gamma F* 3e F/ F+ Beta FSQR 8e F/ F- TO Q ; PRIVATE
: +R! ( -- ) Q -0.5e F* ( a1) Q FSQR 4e F/ P FCUBE 27e F/ F+ FSQRT ( a2) F+ TO +R ; PRIVATE
: +U! ( -- ) +R 0e R,I->Z 3e 1/F 0e R,I->Z Z** TO +U ; PRIVATE
: y! ( -- )
+U 0+0i Z= IF Alpha -5e F* 6e F/ Q FCBRT F- 0e R,I->Z TO y EXIT ENDIF
Alpha -5e F* 6e F/ 0e R,I->Z
+U Z+
P 3e F/ 0e R,I->Z +U Z/ Z- TO y ; PRIVATE
: W! ( -- ) Alpha 0e R,I->Z y 2e Z*F Z+ ZSQRT TO W ; PRIVATE
: SETUP-QUARTIC ( F: a b c d e -- )
TO E TO D TO C TO B TO A
Alpha! Beta! Gamma!
Beta=0? ?EXIT
P! Q! +R! +U! y! W! ; PRIVATE
: COMPUTE-QUARTIC ( -- addr )
Beta F0= IF qsol EXIT ENDIF
Alpha 3e F* 0e R,I->Z y 2e Z*F Z+ ZLOCAL term1
Beta F2* 0e R,I->Z W Z/ ZLOCAL term2
B A F/ -4e F/ 0e R,I->Z ZLOCAL a1
W 2e Z/F ZLOCAL a2
term1 term2 Z+ ZNEGATE ZSQRT 2e Z/F ZLOCAL a3+
term1 term2 Z- ZNEGATE ZSQRT 2e Z/F ZLOCAL a3-
a1 a2 Z+ a3+ Z+ qsol 0 COMPLEX[] Z! \ the sign distribution is tricky
a1 a2 Z+ a3+ Z- qsol 1 COMPLEX[] Z!
a1 a2 Z- a3- Z- qsol 2 COMPLEX[] Z!
a1 a2 Z- a3- Z+ qsol 3 COMPLEX[] Z!
qsol ; PRIVATE
: QUARTIC[] ( ix -- addr ) ( F: a b c d -- ) SETUP-QUARTIC COMPUTE-QUARTIC ;
: .QUARTIC ( addr -- ) 4 0 ?DO Z@+ CR ." x" I 1+ 0 .R ." = " ZS. LOOP DROP ;
NESTING @ 1 =
[IF]
: xeval ZLOCAL x x 4 Z^n x ZSQR 6e Z*F Z+ x 60e Z*F Z- 36e 0e R,I->Z Z+ ;
: test-Cardano \ -- Cardano's problem x^4 + 6x^2 - 60x + 36 = 0
1e 0e 6e -60e 36e SETUP-QUARTIC
CR ." Beta == 0 -> " Beta F0= IF ." TRUE" ELSE ." FALSE" ENDIF
CR ." +R = " +R F.
CR ." +U = " +U ZS.
CR ." W = " W ZS.
CR ." y = " y ZS.
COMPUTE-QUARTIC
4 0 ?DO Z@+ CR ." x" I 1+ 0 .R ." = " ZDUP ZS. ." -> " xeval ZS. LOOP DROP ;
DOC
(*
FORTH> test-cardano
Beta == 0 -> FALSE
+R = 374.1276730966858224403
+U = 7.2056518963939844012e0 + i0.0000000000000000000e0
W = 3.7442732882456315480e0 + i0.0000000000000000000e0
y = 4.0097912285348771276e0 + i0.0000000000000000000e0
x1 = 3.0998744240188162990e0 + i0.0000000000000000000e0 -> 1.9012569296705805754e-14 + i0.0000000000000000000e0
x2 = 6.4439886422681535992e-1 + i0.0000000000000000000e0 -> 7.2615524704389144972e-15 + i0.0000000000000000000e0
x3 = -1.8721366441228157740e0 - i3.8101353367982659924e0 -> 3.6776137690708310402e-15 + i3.3903435614490717854e-14
x4 = -1.8721366441228157740e0 + i3.8101353367982659924e0 -> 3.6776137690708310402e-15 - i3.3903435614490717854e-14 ok
*)
ENDDOC
[THEN]
:ABOUT CR ." Try: 1e 0e 6e -60e 36e QUARTIC[] .QUARTIC -- Solve x^4 + 6x^2 - 60x + 36 = 0"
CR ." Should give:"
CR ." x1 = 3.0998744240188162990e+0 + i0.0000000000000000000e0"
CR ." x2 = 6.4439886422681535992e-1 + i0.0000000000000000000e0"
CR ." x3 = -1.8721366441228157740e+0 - i3.8101353367982659924e0"
CR ." x4 = -1.8721366441228157740e+0 + i3.8101353367982659924e0" ;
.ABOUT -ferrari CR
DEPRIVE
(* End of Source *)
[toc] | [next] | [standalone]
| From | Paul Rubin <no.email@nospam.invalid> |
|---|---|
| Date | 2012-02-23 18:39 -0800 |
| Message-ID | <7xy5rtx8ct.fsf@ruckus.brouhaha.com> |
| In reply to | #9688 |
mhx@iae.nl (Marcel Hendrix) writes: > Here is Ferrari's method to solve a quartic equation. Wow! I've never seen that actually coded before, outside of symbolic algebra packages that figure out whether a given polynomial (maybe degree > 4) has a solvable Galois group and crunch out the formulas. I thought everyone wanting actual numbers from cubic or higher polynomials just uses numerical rootfinders even when an algebraic solution exists, because the algebra is so complicated. Cool!
[toc] | [prev] | [next] | [standalone]
| From | Krishna Myneni <krishna.myneni@ccreweb.org> |
|---|---|
| Date | 2012-02-23 18:58 -0800 |
| Message-ID | <38b01bd7-59ac-49aa-9be8-14995a5dc273@a15g2000yqf.googlegroups.com> |
| In reply to | #9690 |
On Feb 23, 8:39 pm, Paul Rubin <no.em...@nospam.invalid> wrote: > m...@iae.nl (Marcel Hendrix) writes: > > Here is Ferrari's method to solve a quartic equation. > > Wow! I've never seen that actually coded before, outside of symbolic > algebra packages that figure out whether a given polynomial (maybe > degree > 4) has a solvable Galois group and crunch out the formulas. I > thought everyone wanting actual numbers from cubic or higher polynomials > just uses numerical rootfinders even when an algebraic solution exists, > because the algebra is so complicated. Cool! There is a cubic root finder in the FSL (algorithm #6, cubic.x), written by the late Julian Noble, a quadratic root finder (algorithm #63, quadratic.x), and a of course, a numerical root finder for polynomial roots (algorithm #61, lagroots.x), which may be found at the FSL site. Marcel's quartic root finder looks like it would be a great addition. I wonder if the algorithm extends to a quartic with complex coefficients? A revision of the original FSL cubic solver, which is more amenable to use in applications (the original FSL version only displays the roots, but does not return them), may be found at the link below. I've used it, for example, to find eigenvalues of real 3x3 matrices. ftp://ccreweb.org/software/fsl/cubic.fs Krishna
[toc] | [prev] | [next] | [standalone]
| From | mhx@iae.nl (Marcel Hendrix) |
|---|---|
| Date | 2012-02-25 10:11 +0200 |
| Message-ID | <76929310018435@frunobulax.edu> |
| In reply to | #9691 |
Krishna Myneni <krishna.myneni@ccreweb.org> writes Re: Ferrari's method to solve a quartic
> On Feb 23, 8:39 pm, Paul Rubin <no.em...@nospam.invalid> wrote:
>> m...@iae.nl (Marcel Hendrix) writes:
>> > Here is Ferrari's method to solve a quartic equation.
[..]
> I wonder if the algorithm extends to a quartic with
> complex coefficients?
[..]
Ferrari's algorithm (now) also works for complex coefficients.
I fixed some bugs related to Beta being (nearly) zero, and in
case of multiple poles. There are now three test cases: Cardano's,
(z-1)^4, and a case picked up from a discussion/question
on stackoverflow.
Maybe the result can be improved with root polishing?
The previous version (real coefficients) produces better
results for (z-1)^4, i.e. the complex part of the roots
is exactly zero, while it is ~ 1e-10 for zferrari.
Maybe I should code it using the MPC library.
-marcel
-- --------------
(*
* LANGUAGE : ANS Forth with extensions
* PROJECT : Forth Environments
* DESCRIPTION : Ferrari's method to solve a quartic equation
* CATEGORY : Numeric Utility
* AUTHOR : Marcel Hendrix
* LAST CHANGE : February 23, 2012, Marcel Hendrix
* LAST CHANGE : February 24, 2012, Marcel Hendrix; support complex coefficients, handles Beta~0
*)
NEEDS -miscutil
NEEDS -fcbrt
NEEDS -cplx_fsl
REVISION -zferrari "--- Solve a Quartic Version 0.02 ---"
PRIVATES
DOC
(*
The quartic is the highest order polynomial equation that can be solved by
radicals in the general case (i.e., one where the coefficients can take any
value).
Lodovico Ferrari is attributed with the discovery of the solution to the
quartic in 1540, but since this solution, like all algebraic solutions of
the quartic, requires the solution of a cubic to be found, it couldn't be
published immediately. The solution of the quartic was published together
with that of the cubic by Ferrari's mentor Gerolamo Cardano in the book Ars
Magna (1545).
The proof that four is the highest degree of a general polynomial for which
such solutions can be found was first given in the Abel-Ruffini theorem in
1824, proving that all attempts at solving the higher order polynomials
would be futile. The notes left by Evariste Galois prior to dying in a duel
in 1832 later led to an elegant complete theory of the roots of polynomials,
of which this theorem was one result.
( http://en.wikipedia.org/wiki/Quartic_function )
*)
ENDDOC
0+0i ZVALUE A PRIVATE
0+0i ZVALUE B PRIVATE
0+0i ZVALUE C PRIVATE
0+0i ZVALUE D PRIVATE
0+0i ZVALUE E PRIVATE
0+0i ZVALUE Alpha PRIVATE
0+0i ZVALUE Beta PRIVATE
0+0i ZVALUE Gamma PRIVATE
1e-14 FVALUE qeps
0+0i ZVALUE P PRIVATE
0+0i ZVALUE Q PRIVATE
0+0i ZVALUE +R PRIVATE
0+0i ZVALUE +U PRIVATE
0+0i ZVALUE y PRIVATE
0+0i ZVALUE W PRIVATE
CREATE qsol PRIVATE 4 ZFLOATS ALLOT
: Alpha! ( -- ) B A Z/ ZSQR -0.375e Z*F C A Z/ Z+ TO Alpha ; PRIVATE
: Beta! ( -- )
B A Z/ ZCUBE 0.125e Z*F
B C Z* A ZSQR Z*2 Z/ Z-
D A Z/ Z+ TO Beta ; PRIVATE
: Gamma! ( -- )
B A Z/ ZQUAD -3e Z*F 256e Z/F
B ZSQR C Z* A ZCUBE 16e Z*F Z/ Z+
B D Z* A ZSQR 4e Z*F Z/ Z-
E A Z/ Z+ TO Gamma ; PRIVATE
: Beta=0? ( -- bool )
Beta ZABS qeps F> IF FALSE EXIT ENDIF
Alpha ZSQR Gamma 4e Z*F Z- ZSQRT Alpha Z- Z/2 ZSQRT ZLOCAL a2+
Alpha ZSQR Gamma 4e Z*F Z- ZSQRT ZNEGATE Alpha Z- Z/2 ZSQRT ZLOCAL a2-
B A -4e Z*F Z/ ZLOCAL a1
a1 a2+ Z+ qsol 0 COMPLEX[] Z!
a1 a2+ Z- qsol 1 COMPLEX[] Z!
a1 a2- Z+ qsol 2 COMPLEX[] Z!
a1 a2- Z- qsol 3 COMPLEX[] Z!
CLEAR P CLEAR Q CLEAR +R CLEAR +U CLEAR y CLEAR W
TRUE ; PRIVATE
: set-P ( -- ) Alpha ZSQR -12e Z/F Gamma Z- TO P ; PRIVATE
: set-Q ( -- ) Alpha ZCUBE -108e Z/F Alpha Gamma Z* 3e Z/F Z+ Beta ZSQR 8e Z/F Z- TO Q ; PRIVATE
: +R! ( -- ) Q -0.5e Z*F ( a1) Q ZSQR 4e Z/F P ZCUBE 27e Z/F Z+ ZSQRT ( a2) Z+ TO +R ; PRIVATE
: +U! ( -- ) +R ZCBRT TO +U ; PRIVATE
: y! ( -- )
+U 0+0i Z= IF Alpha -5e Z*F 6e Z/F Q ZCBRT Z- TO y EXIT ENDIF
Alpha -5e Z*F 6e Z/F
+U Z+
P 3e Z/F +U Z/ Z- TO y ; PRIVATE
: W! ( -- ) Alpha y Z*2 Z+ ZSQRT TO W ; PRIVATE
: SETUP-QUARTIC ( F: za zb zc zd ze -- )
TO E TO D TO C TO B TO A
Alpha! Beta! Gamma!
Beta=0? ?EXIT
set-P set-Q +R! +U! y! W! ; PRIVATE
: COMPUTE-QUARTIC ( -- addr )
Beta ZABS qeps F< IF qsol EXIT ENDIF
Alpha 3e Z*F y Z*2 Z+ ZLOCAL term1
Beta 2e Z*F W Z/ ZLOCAL term2
B A Z/ -4e Z/F ZLOCAL a1
W Z/2 ZLOCAL a2
term1 term2 Z+ ZNEGATE ZSQRT Z/2 ZLOCAL a3+
term1 term2 Z- ZNEGATE ZSQRT Z/2 ZLOCAL a3-
a1 a2 Z+ a3+ Z+ qsol 0 COMPLEX[] Z! \ the sign distribution is tricky
a1 a2 Z+ a3+ Z- qsol 1 COMPLEX[] Z!
a1 a2 Z- a3- Z- qsol 2 COMPLEX[] Z!
a1 a2 Z- a3- Z+ qsol 3 COMPLEX[] Z!
qsol ; PRIVATE
: QUARTIC[] ( ix -- addr ) ( F: za zb zc zd -- ) SETUP-QUARTIC COMPUTE-QUARTIC ;
: .QUARTIC ( addr -- ) 4 0 ?DO Z@+ CR ." x" I 1+ 0 .R ." = " #18 +ZE.R LOOP DROP ;
NESTING @ 1 =
[IF]
: xeval ( F: z1 -- z2 )
ZLOCAL x
x ZQUAD A Z*
x ZCUBE B Z* Z+
x ZSQR C Z* Z+
x D Z* Z+
E Z+ ;
: .QUARTIC+ ( addr -- ) CR 4 0 ?DO Z@+ CR ." x" I 1+ 0 .R ." = " ZDUP #18 +ZE.R xeval ." --> " #18 +ZE.R LOOP DROP ;
: .INPUTS ( -- )
CR ." A = " A ZS.
CR ." B = " B ZS.
CR ." C = " C ZS.
CR ." D = " D ZS.
CR ." E = " E ZS.
CR ." Alpha = " Alpha ZS.
CR ." Beta = " Beta ZS.
CR ." Gamma = " Gamma ZS.
CR ." Beta == 0 -> " Beta ZABS qeps F< IF ." TRUE" ELSE ." FALSE" ENDIF
CR ." P = " P ZS.
CR ." Q = " Q ZS.
CR ." +R = " +R ZS.
CR ." +U = " +U ZS.
CR ." y = " y ZS.
CR ." W = " W ZS. ;
: test-Cardano \ -- Cardano's problem x^4 + 6x^2 - 60x + 36 = 0
1+0i 0+0i 6e 0e -60e 0e 36e 0e SETUP-QUARTIC
.INPUTS
COMPUTE-QUARTIC
.QUARTIC+ ;
: test-beta \ Beta nearly zero: 0.9604000000000001 x^4 - 5.997600000000001 x^3 + 13.951750054511718 x^2 - 14.326264455924333 x + 5.474214401412618
0.9604000000000001e 0e -5.997600000000001e 0e 13.951750054511718e 0e
-14.326264455924333e 0e 5.474214401412618e 0e SETUP-QUARTIC
.INPUTS
COMPUTE-QUARTIC
.QUARTIC+ ;
: test-(z-1)^4 \ z^4 - 4 z^3 + 6 z^2 - 4 z + 1 = 0
1+0i -4e 0e 6e 0e -4e 0e 1+0i SETUP-QUARTIC
.INPUTS
COMPUTE-QUARTIC
.QUARTIC+ ;
DOC
(*
FORTH> 1+0i 0+0i 6e 0e -60e 0e 36e 0e QUARTIC[] .QUARTIC
x1 = ( 3.099874424018816299e+0000 0.000000000000000000e+0000 )
x2 = ( 6.443988642268153600e-0001 0.000000000000000000e+0000 )
x3 = ( -1.872136644122815774e+0000 -3.810135336798265992e+0000 )
x4 = ( -1.872136644122815774e+0000 3.810135336798265992e+0000 ) ok
FORTH> test-cardano
A = 1.0000000000000000000e0 + i0.0000000000000000000e0
B = 0.0000000000000000000e0 + i0.0000000000000000000e0
C = 6.0000000000000000000e0 + i0.0000000000000000000e0
D = -6.0000000000000000000e1 + i0.0000000000000000000e0
E = 3.6000000000000000000e1 + i0.0000000000000000000e0
Alpha = 6.0000000000000000000e0 + i0.0000000000000000000e0
Beta = -6.0000000000000000000e1 + i0.0000000000000000000e0
Gamma = 3.6000000000000000000e1 + i0.0000000000000000000e0
Beta == 0 -> FALSE
P = -3.9000000000000000000e1 - i0.0000000000000000000e0
Q = -3.8000000000000000000e2 + i0.0000000000000000000e0
+R = 3.7412767309668583948e2 - i0.0000000000000000000e0
+U = 7.2056518963939844012e0 + i0.0000000000000000000e0
y = 4.0097912285348771276e0 + i0.0000000000000000000e0
W = 3.7442732882456315480e0 + i0.0000000000000000000e0
x1 = ( 3.099874424018816299e+0000 0.000000000000000000e+0000 ) --> ( 1.901256929670580576e-0014 0.000000000000000000e+0000 )
x2 = ( 6.443988642268153600e-0001 0.000000000000000000e+0000 ) --> ( 7.261552470438914497e-0015 0.000000000000000000e+0000 )
x3 = ( -1.872136644122815774e+0000 -3.810135336798265992e+0000 ) --> ( 3.677613769070831040e-0015 3.390343561449071785e-0014 )
x4 = ( -1.872136644122815774e+0000 3.810135336798265992e+0000 ) --> ( 3.677613769070831040e-0015 -3.390343561449071785e-0014 ) ok
FORTH> test-beta
A = 9.6040000000000014258e-1 + i0.0000000000000000000e0
B = -5.9976000000000011524e0 + i0.0000000000000000000e0
C = 1.3951750054511718347e1 + i0.0000000000000000000e0
D = -1.4326264455924333063e1 + i0.0000000000000000000e0
E = 5.4742144014126177252e0 + i0.0000000000000000000e0
Alpha = -9.7511396801629485200e-2 + i0.0000000000000000000e0
Beta = 5.9440299904345295092e-15 + i0.0000000000000000000e0
Gamma = -3.4174439518693429388e-3 + i0.0000000000000000000e0
Beta == 0 -> TRUE
P = 0.0000000000000000000e0 + i0.0000000000000000000e0
Q = 0.0000000000000000000e0 + i0.0000000000000000000e0
+R = 0.0000000000000000000e0 + i0.0000000000000000000e0
+U = 0.0000000000000000000e0 + i0.0000000000000000000e0
y = 0.0000000000000000000e0 + i0.0000000000000000000e0
W = 0.0000000000000000000e0 + i0.0000000000000000000e0
x1 = ( 1.914604907169397796e+0000 -0.000000000000000000e+0000 ) --> ( 2.010544508657119422e-0015 0.000000000000000000e+0000 )
x2 = ( 1.207844072422439074e+0000 0.000000000000000000e+0000 ) --> ( -2.020952849512980265e-0015 0.000000000000000000e+0000 )
x3 = ( 1.561224489795918435e+0000 1.654276959321655804e-0001 ) --> ( 0.000000000000000000e+0000 9.434727304968859585e-0016 )
x4 = ( 1.561224489795918435e+0000 -1.654276959321655804e-0001 ) --> ( 0.000000000000000000e+0000 -9.434727304968859585e-0016 ) ok
FORTH> test-(z-1)^4
A = 1.0000000000000000000e0 + i0.0000000000000000000e0
B = -4.0000000000000000000e0 + i0.0000000000000000000e0
C = 6.0000000000000000000e0 + i0.0000000000000000000e0
D = -4.0000000000000000000e0 + i0.0000000000000000000e0
E = 1.0000000000000000000e0 + i0.0000000000000000000e0
Alpha = 4.3368086899420177360e-19 + i0.0000000000000000000e0
Beta = 4.3368086899420177360e-19 + i0.0000000000000000000e0
Gamma = 0.0000000000000000000e0 + i0.0000000000000000000e0
Beta == 0 -> TRUE
P = 0.0000000000000000000e0 + i0.0000000000000000000e0
Q = 0.0000000000000000000e0 + i0.0000000000000000000e0
+R = 0.0000000000000000000e0 + i0.0000000000000000000e0
+U = 0.0000000000000000000e0 + i0.0000000000000000000e0
y = 0.0000000000000000000e0 + i0.0000000000000000000e0
W = 0.0000000000000000000e0 + i0.0000000000000000000e0
x1 = ( 1.000000000000000000e+0000 0.000000000000000000e+0000 ) --> ( 0.000000000000000000e+0000 0.000000000000000000e+0000 )
x2 = ( 1.000000000000000000e+0000 -0.000000000000000000e+0000 ) --> ( 0.000000000000000000e+0000 0.000000000000000000e+0000 )
x3 = ( 1.000000000000000000e+0000 6.585445079827192917e-0010 ) --> ( 0.000000000000000000e+0000 4.038967834731580444e-0028 )
x4 = ( 1.000000000000000000e+0000 -6.585445079827192917e-0010 ) --> ( 0.000000000000000000e+0000 -4.038967834731580444e-0028 ) ok
*)
ENDDOC
[THEN]
:ABOUT CR ." Try: 1+0i 0+0i 6e 0e -60e 0e 36e 0e QUARTIC[] .QUARTIC -- Solve x^4 + 6x^2 - 60x + 36 = 0"
CR ." Should give:"
CR ." x1 = ( 3.0998744240188163e+0000 0.0000000000000000e+0000 )"
CR ." x2 = ( 6.4439886422681536e-0001 0.0000000000000000e+0000 )"
CR ." x3 = ( -1.8721366441228158e+0000 -3.8101353367982660e+0000 )"
CR ." x4 = ( -1.8721366441228158e+0000 3.8101353367982660e+0000 )" ;
.ABOUT -zferrari CR
DEPRIVE
(* End of Source *)
[toc] | [prev] | [next] | [standalone]
| From | Krishna Myneni <krishna.myneni@ccreweb.org> |
|---|---|
| Date | 2012-02-28 18:59 -0800 |
| Message-ID | <947e2cb1-aeaa-484e-aed2-6bc821708687@p12g2000yqe.googlegroups.com> |
| In reply to | #9707 |
On Feb 25, 2:11 am, m...@iae.nl (Marcel Hendrix) wrote: > Krishna Myneni <krishna.myn...@ccreweb.org> writes Re: Ferrari's method to solve a quartic > > > On Feb 23, 8:39 pm, Paul Rubin <no.em...@nospam.invalid> wrote: > >> m...@iae.nl (Marcel Hendrix) writes: > >> > Here is Ferrari's method to solve a quartic equation. > [..] > > I wonder if the algorithm extends to a quartic with > > complex coefficients? > > [..] > > Ferrari's algorithm (now) also works for complex coefficients. > I fixed some bugs related to Beta being (nearly) zero, and in > case of multiple poles. There are now three test cases: Cardano's, > (z-1)^4, and a case picked up from a discussion/question > on stackoverflow. > > Maybe the result can be improved with root polishing? > > The previous version (real coefficients) produces better > results for (z-1)^4, i.e. the complex part of the roots > is exactly zero, while it is ~ 1e-10 for zferrari. > Maybe I should code it using the MPC library. > > -marcel ... Great! I'll update my standard Forth version (quartic.fs) soon. I'll also see if I can find some additional test cases, when time permits. Thanks for posting this and the MPC library version. Krishna
[toc] | [prev] | [next] | [standalone]
| From | mhx@iae.nl (Marcel Hendrix) |
|---|---|
| Date | 2012-02-25 15:32 +0200 |
| Message-ID | <13718810018435@frunobulax.edu> |
| In reply to | #9691 |
Marcel Hendrix wrote:
> Krishna Myneni <krishna.myneni@ccreweb.org> writes Re: Ferrari's method to solve a quartic
>> On Feb 23, 8:39 pm, Paul Rubin <no.em...@nospam.invalid> wrote:
>>> m...@iae.nl (Marcel Hendrix) writes:
>>> > Here is Ferrari's method to solve a quartic equation.
[..]
>> I wonder if the algorithm extends to a quartic with
>> complex coefficients?
[..]
> The previous version (real coefficients) produces better
> results for (z-1)^4, i.e. the complex part of the roots
> is exactly zero, while it is ~ 1e-10 for zferrari.
> Maybe I should code it using the MPC library.
Here it is.
-marcel
-- ----------------------------------
(*
* LANGUAGE : ANS Forth with extensions
* PROJECT : Forth Environments
* DESCRIPTION : Ferrari's method to solve a quartic equation using MPC
* CATEGORY : Numeric Utility
* AUTHOR : Marcel Hendrix
* LAST CHANGE : February 23, 2012, Marcel Hendrix
* LAST CHANGE : February 24, 2012, Marcel Hendrix; support complex coefficients, handles Beta~0
* LAST CHANGE : Saturday, February 25, 2012, 12:36, Marcel Hendrix; for MPC
*)
NEEDS -miscutil
NEEDS -mpc
REVISION -zferrari "--- Solve a Quartic Version 0.04 ---"
PRIVATES
DOC
(*
The quartic is the highest order polynomial equation that can be solved by
radicals in the general case (i.e., one where the coefficients can take any
value).
Lodovico Ferrari is attributed with the discovery of the solution to the
quartic in 1540, but since this solution, like all algebraic solutions of
the quartic, requires the solution of a cubic to be found, it couldn't be
published immediately. The solution of the quartic was published together
with that of the cubic by Ferrari's mentor Gerolamo Cardano in the book Ars
Magna (1545).
The proof that four is the highest degree of a general polynomial for which
such solutions can be found was first given in the Abel-Ruffini theorem in
1824, proving that all attempts at solving the higher order polynomials
would be futile. The notes left by Evariste Galois prior to dying in a duel
in 1832 later led to an elegant complete theory of the roots of polynomials,
of which this theorem was one result.
( http://en.wikipedia.org/wiki/Quartic_function )
*)
ENDDOC
0+0i Z#VALUE A PRIVATE
0+0i Z#VALUE B PRIVATE
0+0i Z#VALUE C PRIVATE
0+0i Z#VALUE D PRIVATE
0+0i Z#VALUE E PRIVATE
0+0i Z#VALUE Alpha PRIVATE
0+0i Z#VALUE Beta PRIVATE
0+0i Z#VALUE Gamma PRIVATE
S" 1e-40" F#IN F#VALUE qeps
S" -0.375" F#IN F#CONSTANT F#-3/8 PRIVATE
0+0i Z#VALUE P PRIVATE
0+0i Z#VALUE Q PRIVATE
0+0i Z#VALUE +R PRIVATE
0+0i Z#VALUE +U PRIVATE
0+0i Z#VALUE y PRIVATE
0+0i Z#VALUE W PRIVATE
0+0i Z#VALUE term1 PRIVATE
0+0i Z#VALUE term2 PRIVATE
0+0i Z#VALUE a1 PRIVATE
0+0i Z#VALUE a2 PRIVATE
0+0i Z#VALUE a2+ PRIVATE
0+0i Z#VALUE a2- PRIVATE
0+0i Z#VALUE a3+ PRIVATE
0+0i Z#VALUE a3- PRIVATE
4 Z#BLOCK qsol PRIVATE
: Alpha! ( -- ) B A Z#/ Z#SQR F#-3/8 Z#*F C A Z#/ Z#+ TO Alpha ; PRIVATE
: Beta! ( -- )
B A Z#/ Z#CUBE 8 Z#U/
B C Z#* A Z#SQR Z#2* Z#/ Z#-
D A Z#/ Z#+ TO Beta ; PRIVATE
: Gamma! ( -- )
B A Z#/ Z#QUAD -3 Z#N* #256 Z#U/
B Z#SQR C Z#* A Z#CUBE 16 Z#N* Z#/ Z#+
B D Z#* A Z#SQR 4 Z#N* Z#/ Z#-
E A Z#/ Z#+ TO Gamma ; PRIVATE
: Beta=0? ( -- bool )
Beta Z#ABS qeps F#> IF FALSE EXIT ENDIF
Alpha Z#SQR Gamma 4 Z#N* Z#- Z#SQRT Alpha Z#- Z#2/ Z#SQRT TO a2+
Alpha Z#SQR Gamma 4 Z#N* Z#- Z#SQRT Z#NEGATE Alpha Z#- Z#2/ Z#SQRT TO a2-
B A -4 Z#N* Z#/ TO a1
a1 a2+ Z#+ qsol 0 Z#[] Z#!
a1 a2+ Z#- qsol 1 Z#[] Z#!
a1 a2- Z#+ qsol 2 Z#[] Z#!
a1 a2- Z#- qsol 3 Z#[] Z#!
CLEAR P CLEAR Q CLEAR +R CLEAR +U CLEAR y CLEAR W
TRUE ; PRIVATE
: set-P ( -- ) Alpha Z#SQR #12 Z#U/ Z#NEGATE Gamma Z#- TO P ; PRIVATE
: set-Q ( -- ) Alpha Z#CUBE #108 Z#U/ Z#NEGATE Alpha Gamma Z#* 3 Z#U/ Z#+ Beta Z#SQR 8 Z#U/ Z#- TO Q ; PRIVATE
: +R! ( -- ) Q Z#2/ Z#NEGATE ( a1) Q Z#SQR 4 Z#U/ P Z#CUBE 27 Z#U/ Z#+ Z#SQRT ( a2) Z#+ TO +R ; PRIVATE
: +U! ( -- ) +R Z#CBRT TO +U ; PRIVATE
: y! ( -- )
+U Z#ABS qeps F#< IF Alpha -5 Z#N* 6 Z#U/ Q Z#CBRT Z#- TO y EXIT ENDIF
Alpha -5 Z#N* 6 Z#U/
+U Z#+
P 3 Z#U/ +U Z#/ Z#- TO y ; PRIVATE
: W! ( -- ) Alpha y Z#2* Z#+ Z#SQRT TO W ; PRIVATE
: SETUP-QUARTIC ( F: za zb zc zd ze -- )
TO E TO D TO C TO B TO A
Alpha! Beta! Gamma!
Beta=0? ?EXIT
set-P set-Q +R! +U! y! W! ; PRIVATE
: COMPUTE-QUARTIC ( -- addr )
Beta Z#ABS qeps F#< IF qsol EXIT ENDIF
Alpha 3 Z#N* y Z#2* Z#+ TO term1
Beta Z#2* W Z#/ TO term2
B A Z#/ 4 Z#U/ Z#NEGATE TO a1
W Z#2/ TO a2
term1 term2 Z#+ Z#NEGATE Z#SQRT Z#2/ TO a3+
term1 term2 Z#- Z#NEGATE Z#SQRT Z#2/ TO a3-
a1 a2 Z#+ a3+ Z#+ qsol 0 Z#[] Z#! \ the sign distribution is tricky
a1 a2 Z#+ a3+ Z#- qsol 1 Z#[] Z#!
a1 a2 Z#- a3- Z#- qsol 2 Z#[] Z#!
a1 a2 Z#- a3- Z#+ qsol 3 Z#[] Z#!
qsol ; PRIVATE
: QUARTIC ( -- addr ) ( F: za zb zc zd -- ) SETUP-QUARTIC COMPUTE-QUARTIC ;
: .QUARTIC ( addr -- ) LOCAL q[] 4 0 ?DO q[] I Z#[] Z#@ CR ." x" I 1+ 0 .R ." = " Z#. LOOP ;
NESTING @ 1 =
[IF]
0+0i Z#VALUE x
: xeval ( F: z1 -- z2 )
TO x
x Z#QUAD A Z#*
x Z#CUBE B Z#* Z#+
x Z#SQR C Z#* Z#+
x D Z#* Z#+
E Z#+ ;
: .QUARTIC+ ( addr -- ) LOCAL q[] CR 4 0 ?DO q[] I Z#[] Z#@ CR ." x" I 1+ 0 .R ." = " Z#DUP Z#. xeval CR ." --> " Z#. LOOP ;
: .INPUTS ( -- )
CR ." A = " A Z#.
CR ." B = " B Z#.
CR ." C = " C Z#.
CR ." D = " D Z#.
CR ." E = " E Z#.
CR ." Alpha = " Alpha Z#.
CR ." Beta = " Beta Z#.
CR ." Gamma = " Gamma Z#.
CR ." Beta == 0 -> " Beta Z#ABS qeps F#< IF ." TRUE" ELSE ." FALSE" ENDIF
CR ." P = " P Z#.
CR ." Q = " Q Z#.
CR ." +R = " +R Z#.
CR ." +U = " +U Z#.
CR ." y = " y Z#.
CR ." W = " W Z#. ;
S" ( -4 0)" Z#IN Z#CONSTANT Z#-4.0
S" ( 6 0)" Z#IN Z#CONSTANT Z#6.0
S" (-60 0)" Z#IN Z#CONSTANT Z#-60.0
S" ( 36 0)" Z#IN Z#CONSTANT Z#36.0
S" ( 0.9604000000000001 0)" Z#IN Z#CONSTANT b1
S" ( -5.997600000000001 0)" Z#IN Z#CONSTANT b2
S" ( 13.951750054511718 0)" Z#IN Z#CONSTANT b3
S" (-14.326264455924333 0)" Z#IN Z#CONSTANT b4
S" ( 5.474214401412618 0)" Z#IN Z#CONSTANT b5
: test-Cardano \ -- Cardano's problem x^4 + 6x^2 - 60x + 36 = 0
1+0i 0+0i Z#6.0 Z#-60.0 Z#36.0 SETUP-QUARTIC
.INPUTS
COMPUTE-QUARTIC
.QUARTIC+ ;
: test-beta \ Beta nearly zero: 0.9604000000000001 x^4 - 5.997600000000001 x^3 + 13.951750054511718 x^2 - 14.326264455924333 x + 5.474214401412618
b1 b2 b3 b4 b5 SETUP-QUARTIC
.INPUTS
COMPUTE-QUARTIC
.QUARTIC+ ;
: test-(z-1)^4 \ z^4 - 4 z^3 + 6 z^2 - 4 z + 1 = 0
1+0i Z#-4.0 Z#6.0 Z#-4.0 1+0i SETUP-QUARTIC
.INPUTS
COMPUTE-QUARTIC
.QUARTIC+ ;
DOC
(*
FORTH> 1+0i 0+0i Z#6.0 Z#-60.0 Z#36.0 QUARTIC .QUARTIC
x1 = 3.0998744240188161015876924404485028457423e+0000 0.0000000000000000000000000000000000000000e-0001 i
x2 = 6.4439886422681550176229551792942104128268e-0001 0.0000000000000000000000000000000000000000e-0001 i
x3 = -1.8721366441228158016749939791889619435125e+0000 3.8101353367982661465104486477726884910560e+0000 i
x4 = -1.8721366441228158016749939791889619435125e+0000 -3.8101353367982661465104486477726884910560e+0000 i ok
FORTH> test-cardano
A = 1.0000000000000000000000000000000000000000e+0000 0.0000000000000000000000000000000000000000e-0001 i
B = 0.0000000000000000000000000000000000000000e-0001 0.0000000000000000000000000000000000000000e-0001 i
C = 6.0000000000000000000000000000000000000000e+0000 0.0000000000000000000000000000000000000000e-0001 i
D = -6.0000000000000000000000000000000000000000e+0001 0.0000000000000000000000000000000000000000e-0001 i
E = 3.6000000000000000000000000000000000000000e+0001 0.0000000000000000000000000000000000000000e-0001 i
Alpha = 6.0000000000000000000000000000000000000000e+0000 0.0000000000000000000000000000000000000000e-0001 i
Beta = -6.0000000000000000000000000000000000000000e+0001 0.0000000000000000000000000000000000000000e-0001 i
Gamma = 3.6000000000000000000000000000000000000000e+0001 0.0000000000000000000000000000000000000000e-0001 i
Beta == 0 -> FALSE
P = -3.9000000000000000000000000000000000000000e+0001 -0.0000000000000000000000000000000000000000e-0001 i
Q = -3.8000000000000000000000000000000000000000e+0002 0.0000000000000000000000000000000000000000e-0001 i
+R = 3.7412767309668582242564878268427340344589e+0002 -0.0000000000000000000000000000000000000000e-0001 i
+U = 7.2056518963939846930775491804063003724888e+0000 -0.0000000000000000000000000000000000000000e-0001 i
y = 4.0097912285348773231421386699582591268932e+0000 0.0000000000000000000000000000000000000000e-0001 i
W = 3.7442732882456316033499879583779238870250e+0000 0.0000000000000000000000000000000000000000e-0001 i
x1 = 3.0998744240188161015876924404485028457423e+0000 0.0000000000000000000000000000000000000000e-0001 i
--> 2.2108591501041778240989060768769022902056e-0075 0.0000000000000000000000000000000000000000e-0001 i
x2 = 6.4439886422681550176229551792942104128268e-0001 0.0000000000000000000000000000000000000000e-0001 i
--> 0.0000000000000000000000000000000000000000e-0001 0.0000000000000000000000000000000000000000e-0001 i
x3 = -1.8721366441228158016749939791889619435125e+0000 3.8101353367982661465104486477726884910560e+0000 i
--> 0.0000000000000000000000000000000000000000e-0001 6.6325774503125334722967182306307068706170e-0075 i
x4 = -1.8721366441228158016749939791889619435125e+0000 -3.8101353367982661465104486477726884910560e+0000 i
--> 0.0000000000000000000000000000000000000000e-0001 -6.6325774503125334722967182306307068706170e-0075 i ok
FORTH> test-beta
A = 9.6040000000000010000000000000000000000000e-0001 0.0000000000000000000000000000000000000000e-0001 i
B = -5.9976000000000009999999999999999999999999e+0000 0.0000000000000000000000000000000000000000e-0001 i
C = 1.3951750054511717999999999999999999999999e+0001 0.0000000000000000000000000000000000000000e-0001 i
D = -1.4326264455924333000000000000000000000000e+0001 0.0000000000000000000000000000000000000000e-0001 i
E = 5.4742144014126180000000000000000000000000e+0000 0.0000000000000000000000000000000000000000e-0001 i
Alpha = -9.7511396801629749551421474882146550948804e-0002 0.0000000000000000000000000000000000000000e-0001 i
Beta = 5.3631349419021963322480553222902379073137e-0015 0.0000000000000000000000000000000000000000e-0001 i
Gamma = -3.4174439518694905579812191628717883693476e-0003 0.0000000000000000000000000000000000000000e-0001 i
Beta == 0 -> FALSE
P = 2.6250712430190831787821175823593875669108e-0003 0.0000000000000000000000000000000000000000e-0001 i
Q = 1.1966495214908011499374682316970769466869e-0004 -0.0000000000000000000000000000000000000000e-0001 i
+R = 5.3587933955913037336767600334679800425037e-0006 0.0000000000000000000000000000000000000000e-0001 i
+U = 1.7499366938287110744897352993709535159728e-0002 0.0000000000000000000000000000000000000000e-0001 i
y = 4.8755698400814874775710738061551924589409e-0002 0.0000000000000000000000000000000000000000e-0001 i
W = 3.5227223822351018951482546263270123028078e-0014 0.0000000000000000000000000000000000000000e-0001 i
x1 = 1.5612244897959360787072371026316680470492e+0000 -1.6542769593216835969789894020584464029664e-0001 i
--> -4.2123274051525879873007970023884313331788e-0054 3.4544674220377778501545407451201598284464e-0077 i
x2 = 1.5612244897959360787072371026316680470492e+0000 1.6542769593216835969789894020584464029664e-0001 i
--> -4.2123274051525879873007970023884313331788e-0054 -3.4544674220377778501545407451201598284464e-0077 i
x3 = 1.2078440724224197532447709413299479764843e+0000 0.0000000000000000000000000000000000000000e-0001 i
--> -4.2123274051525879873010733597821943554068e-0054 0.0000000000000000000000000000000000000000e-0001 i
x4 = 1.9146049071693819497220585618954851525216e+0000 -0.0000000000000000000000000000000000000000e-0001 i
--> -4.2123274051525879873013497171759573776348e-0054 0.0000000000000000000000000000000000000000e-0001 i ok
FORTH> test-(z-1)^4
A = 1.0000000000000000000000000000000000000000e+0000 0.0000000000000000000000000000000000000000e-0001 i
B = -4.0000000000000000000000000000000000000000e+0000 0.0000000000000000000000000000000000000000e-0001 i
C = 6.0000000000000000000000000000000000000000e+0000 0.0000000000000000000000000000000000000000e-0001 i
D = -4.0000000000000000000000000000000000000000e+0000 0.0000000000000000000000000000000000000000e-0001 i
E = 1.0000000000000000000000000000000000000000e+0000 0.0000000000000000000000000000000000000000e-0001 i
Alpha = 0.0000000000000000000000000000000000000000e-0001 0.0000000000000000000000000000000000000000e-0001 i
Beta = 0.0000000000000000000000000000000000000000e-0001 0.0000000000000000000000000000000000000000e-0001 i
Gamma = 0.0000000000000000000000000000000000000000e-0001 0.0000000000000000000000000000000000000000e-0001 i
Beta == 0 -> TRUE
P = 0.0000000000000000000000000000000000000000e-0001 0.0000000000000000000000000000000000000000e-0001 i
Q = 0.0000000000000000000000000000000000000000e-0001 0.0000000000000000000000000000000000000000e-0001 i
+R = 0.0000000000000000000000000000000000000000e-0001 0.0000000000000000000000000000000000000000e-0001 i
+U = 0.0000000000000000000000000000000000000000e-0001 0.0000000000000000000000000000000000000000e-0001 i
y = 0.0000000000000000000000000000000000000000e-0001 0.0000000000000000000000000000000000000000e-0001 i
W = 0.0000000000000000000000000000000000000000e-0001 0.0000000000000000000000000000000000000000e-0001 i
x1 = 1.0000000000000000000000000000000000000000e+0000 0.0000000000000000000000000000000000000000e-0001 i
--> 0.0000000000000000000000000000000000000000e-0001 0.0000000000000000000000000000000000000000e-0001 i
x2 = 1.0000000000000000000000000000000000000000e+0000 -0.0000000000000000000000000000000000000000e-0001 i
--> 0.0000000000000000000000000000000000000000e-0001 0.0000000000000000000000000000000000000000e-0001 i
x3 = 1.0000000000000000000000000000000000000000e+0000 -0.0000000000000000000000000000000000000000e-0001 i
--> 0.0000000000000000000000000000000000000000e-0001 0.0000000000000000000000000000000000000000e-0001 i
x4 = 1.0000000000000000000000000000000000000000e+0000 0.0000000000000000000000000000000000000000e-0001 i
--> 0.0000000000000000000000000000000000000000e-0001 0.0000000000000000000000000000000000000000e-0001 i ok
*)
ENDDOC
[THEN]
:ABOUT CR ." Try: 1+0i 0+0i Z#6.0 Z#-60.0 Z#36.0 QUARTIC .QUARTIC -- Solve x^4 + 6x^2 - 60x + 36 = 0"
CR ." Should give:"
CR ." x1 = ( 3.0998744240188163e+0000 0.0000000000000000e+0000 )"
CR ." x2 = ( 6.4439886422681536e-0001 0.0000000000000000e+0000 )"
CR ." x3 = ( -1.8721366441228158e+0000 -3.8101353367982660e+0000 )"
CR ." x4 = ( -1.8721366441228158e+0000 3.8101353367982660e+0000 )" ;
.ABOUT -zferrari CR
DEPRIVE
(* End of Source *)
[toc] | [prev] | [next] | [standalone]
| From | Krishna Myneni <krishna.myneni@ccreweb.org> |
|---|---|
| Date | 2012-02-23 20:25 -0800 |
| Message-ID | <941ee05e-1c35-4f84-a564-eba9a6d6d6d0@f2g2000yqh.googlegroups.com> |
| In reply to | #9688 |
On Feb 23, 6:52 pm, m...@iae.nl (Marcel Hendrix) wrote: > Here is Ferrari's method to solve a quartic equation. > > I only tested it with Cardano's problem, where the absolute > error of the roots is about 10^-14 to 10^-15. This high > error is caused by the use of complex numbers (implemented > with standard double precision in iForth). > > -marcel ... Here's a version in standard Forth, using the FSL complex arithmetic library, (algorithm #60, complex.fs). ftp://ccreweb.org/software/fsl/extras/quartic.fs KM
[toc] | [prev] | [standalone]
Back to top | Article view | comp.lang.forth
csiph-web