C    Copyright 1984 by ABComputing				   May 15, 1984
C--------------------------------------------------------------------------
C
C    ÚÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ¿
C    ³									 ³
C    ³	 EDITOR'S NOTE:                                                  ³
C    ³									 ³
C    ³	 THIS IS THE MAIN PROGRAM  FOR SOLVING CUBIC AND FOURTH	 ORDER	 ³
C    ³	 POLYNOMIALS, WITH REAL	COEFFICIENTS.				 ³
C    ³									 ³
C    ³	 THIS CODE WAS WRITTEN BY  THE AUTHOR ON A VAX,	 AND PRESENTED	 ³
C    ³	 TO ME IN HARDCOPY  FORM.  I HAVE REARRANGED  CERTAIN SECTIONS	 ³
C    ³	 OF THE	 CODE, TYPED  THIS SOURCE  CODE, AND  MUST ACCEPT FULL	 ³
C    ³	 BLAME FOR ANY ERRORS.						 ³
C    ³									 ³
C    ÀÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÙ
C
C--------------------------------------------------------------------------



      LOGICAL	 L
      COMPLEX	 P3, P4, ROOTS,	Z
      DIMENSION	 ROOTS(4)

C     --NESTED FORM OF CUBIC AND QUARTIC
      P3(Z) = C0 + Z*(C1+ Z*(C2+Z))
      P4(Z) = C0 + Z*(C1+ Z*(C2+Z*(C3+Z)))


C     --START OF MAIN PROGRAM.

10    PRINT*, 'ENTER .FALSE. FOR CUBIC; .TRUE. FOR QUARTIC'
      ACCEPT*,L
      IF (L) GOTO 30

C ---------------

      PRINT*, 'ENTER COEFFICIENTS OF CUBIC: C0,C1,C2'
      ACCEPT*,C0,C1,C2

      CALL CUBIC(C0,C1,C2,ROOTS)

C     --PRINT ROOTS OF CUBIC AND VALUES	OF P3 AT ROOTS.

      DO 20 I =	1,3
	   Z = ROOTS(I)
	   V = CABS(P3(Z))
20	   PRINT*,I,Z,V
      GOTO 10

C ---------------

30    PRINT*, 'ENTER COEFFICIENTS OF QUARTIC: C0,C1,C2,C3'
      ACCEPT*,C0,C1,C2,C3

      CALL QURTC(C0,C1,C2,C3,ROOTS)

C     --PRINT ROOTS OF QUARTIC AND VALUES OF P4	AT ROOTS.

      DO 40 I =	1,4
	   Z = ROOTS(I)
	   V = CABS(P4(Z))
40	   PRINT*,I,Z,V
      GOTO 10

      END
C--------------------------------------------------------------------------


C--------------------------------------------------------------------------
      FUNCTION CBRT(X)

C     --THIS FUNCTION RETURNS THE CUBE ROOT OF ITS ARGUMENT, EVEN IF X <0.

      IF ( X .EQ. 0.) THEN
	   CBRT	= 0.
      ELSE
	   CBRT	= SIGN(ABS(X)**.3333333,X)
      ENDIF

      RETURN
      END
C--------------------------------------------------------------------------


C--------------------------------------------------------------------------
      SUBROUTINE CUBIC(C0,C1,C2,ROOTS)

C     --RETURNS	THE ROOTS OF THE CUBIC:	C0 + C1*X + C2*X**2 + X**3 = 0.


      COMPLEX	 ROOTS
      DIMENSION	 ROOTS(3)

C     --CUR RETURNS A REAL ROOT	OF THE CUBIC
      R	= CUR(C0,C1,C2)
      ROOTS(1) = R

C     --SOLVE QUADRATICS REMAINING AFTER FACTORING OUT ROOT
      CALL QDRTC(C1+R*(C2+R),C2+R,ROOTS(2),ROOTS(3))

      RETURN
      END
C--------------------------------------------------------------------------


C--------------------------------------------------------------------------
      FUNCTION CUR(C0,C1,C2)

C     --THIS FUNCTION RETURNS A	REAL ROOT OF P3.

      COMPLEX U

      G	= ( 2.*C2*C2*C2	- 9.C2*C1 + 27.*C0)/54.
      H	= ( 3.*C1 -C2*C2)/9.
      ARG = G*G	+ H*H*H

C     --ARG IS ESSENTIALLY THE DISCRIMINANT(D).	IF D < 0 THERE ARE
C     THREE REAL ROOTS AND THE CUBE ROOT OF A COMPLEX NUMBER IS	NEEDED.

      IF (ARG .LT. 0. )	 THEN
	   U = (CSQRT(CMPLX(ARG,0.)) - G )** .3333333
	   CUR = U + CONJG(U) -	C2/3.
	   RETURN
      ENDIF

C     --A POSITIVE DISCRIMINANT	MEANS ONLY ONE REAL ROOT. A ZERO DISCRIMINANT
C     MEANS MULTIPLE REAL ROOTS.

      IF (ARG .GT. 0. )	 THEN
	   V = CBRT( -G	+ SIGN(SQRT(ARG), -G ))
	   CUR = V - H/V - C2/3
      ELSE
	   CUR = 2.*CBRT(-G) - C2/3.
      ENDIF

      RETURN
      END
C--------------------------------------------------------------------------


C--------------------------------------------------------------------------
      SUBROUTINE QDRTC(C0,C1,R1,R2)

C     --RETURNS	THE ROOTS OF THE QUADRATIC: C0 + C1*X +	X**2 = 0.

      COMPLEX  R1, R2

      B	= -C1/2.
      D	= B*B -	C0

C     --IF D = 0, THEN WE HAVE A DOUBLE	ROOT. IF D < 0,	THEN HAVE COMPLEX
C     ROOTS.

      IF( D .EQ. 0.) THEN
	   R1 =	B
	   R2 =	B
	   RETURN
      ENDIF

      IF( D .LT. 0.) THEN
	   R1 =	CMPLX(B,SQRT(-D))
	   R2 =	CONJG(R1)
	   RETURN
      ENDIF

C     --CALCULATE REAL ROOTS WITH MINIMUM ROUNDING

      IF( B .GE. 0. ) THEN
	   R1 =	B + SQRT(D)
      ELSE
	   R1 =	B - SQRT(D)
      ENDIF

      R2 = C0/R1

      RETURN
      END
C--------------------------------------------------------------------------


C--------------------------------------------------------------------------
      SUBROUTINE QURTC(C0,C1,C2,C3,ROOTS)

C     --FIND THE ROOTS OF  C0 +	C1*X + C2*X**2 + C3*X**3 + X**4	= 0.

      COMPLEX  ROOTS
      DIMENSION	ROOTS(4)

C     --RETURNS	A REAL ROOT OF THE RESOLVENT CUBIC
      R	= CUR((	4*C2*C0-C3*C3*C0-C1*C1)/8,C3*C1/4-C0,-C2/2)
      B	= R*R -	C0

C     --THE IF CONSTRUCT INSURES THE "+" ARGUMENT FOR SQRT

      IF ( B .GT. 1E-6)	THEN
	   B = SQRT(B)
	   A = ( R*C3-C1)/2./B
      ELSE
	   A = SQRT( 2.*R + C3*C3/4. - C2)
	   B = (R*C3 - C1)/2./A
      ENDIF

C     --SOLVE THE TWO QUADRATIC	FACTORS

      CALL QDRTC( R-B,C3/2-A,ROOTS(1),ROOTS(2))
      CALL QDRTC( R+B,C3/2+A,ROOTS(3),ROOTS(4))

      RETURN
      END




		     ÚÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄ¿
		     ³ File Name:  ÛÛ	list1.for   ÛÛ	³
		     ÀÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÄÙ