Path: csiph.com!x330-a1.tempe.blueboxinc.net!usenet.pasdenom.info!weretis.net!feeder4.news.weretis.net!ecngs!feeder2.ecngs.de!newsfeed.freenet.ag!news2.euro.net!postnews2.euro.net!news.wanadoo.nl!not-for-mail From: mhx@iae.nl (Marcel Hendrix) Subject: Re: Ferrari's method to solve a quartic Newsgroups: comp.lang.forth Message-ID: <76929310018435@frunobulax.edu> Date: Sat, 25 Feb 2012 10:11:59 +0200 References: <38b01bd7-59ac-49aa-9be8-14995a5dc273@a15g2000yqf.googlegroups.com> X-Newsreader: iForth 2.0 console (October 21, 2006) Lines: 287 Organization: Wanadoo NNTP-Posting-Date: 25 Feb 2012 09:11:31 GMT NNTP-Posting-Host: s529d937f.adsl.wanadoo.nl X-Trace: 1330161091 dr6.euro.net 233 82.157.147.127:54322 X-Complaints-To: abuse@wanadoo.nl Xref: x330-a1.tempe.blueboxinc.net comp.lang.forth:9707 Krishna Myneni writes Re: Ferrari's method to solve a quartic > On Feb 23, 8:39 pm, Paul Rubin 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 *)