Re: Ferrari's method to solve a quartic

From	mhx@iae.nl (Marcel Hendrix)
Subject	Re: Ferrari's method to solve a quartic
Newsgroups	comp.lang.forth
Message-ID	<76929310018435@frunobulax.edu> (permalink)
Date	2012-02-25 10:11 +0200
References	<38b01bd7-59ac-49aa-9be8-14995a5dc273@a15g2000yqf.googlegroups.com>
Organization	Wanadoo

Show all headers | View raw

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 *)

Thread

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

csiph-web