Ferrari's method to solve a quartic

From	mhx@iae.nl (Marcel Hendrix)
Subject	Ferrari's method to solve a quartic
Newsgroups	comp.lang.forth
Message-ID	<14830111018435@frunobulax.edu> (permalink)
Date	2012-02-24 02:52 +0200
Organization	Wanadoo

Show all headers | View raw

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

Back to comp.lang.forth | Previous | Next — 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