Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.forth > #9708

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 <13718810018435@frunobulax.edu> (permalink)
Date 2012-02-25 15:32 +0200
References <38b01bd7-59ac-49aa-9be8-14995a5dc273@a15g2000yqf.googlegroups.com>
Organization Wanadoo

Show all headers | View raw

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

Back to comp.lang.forth | Previous | Next — Previous in thread | Next in thread | Find similar | Unroll thread

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