PROGRAM Newton
    USE double
    IMPLICIT NONE
    REAL(KIND=DP) :: target, current
    REAL(KIND=DP), DIMENSION(5) :: previous
    INTEGER :: i, k
mainloop: DO
        PRINT *, 'Type in a real number'
        READ (*, *, IOSTAT=k) target
        IF (k < 0) THEN
            STOP
        ELSE IF (k > 0) THEN
            PRINT *, 'Some sort of horrible I/O error'
            STOP
        END IF 
        IF (target .EQ. 0.0_DP) THEN
            PRINT *, 'The cube root of zero is, er, zero'
            CYCLE
        END IF
!
! This is cheating, but the hardest part of Newton-Raphson solution of
! Nth roots is getting the starting value, and doing so properly would
! merely be confusing.  So use a horrible hack to get an approximation.
!
        current = 1.1*CMPLX(target,KIND=KIND(0.0))**0.3
        DO i = 1,5
            previous(i) = 0.0_DP
        END DO
loop:   DO
            current = current - &
                value(current,target)/derivative(current)
            PRINT *, current
            DO i = 1,5
                if (current .EQ. previous(i)) EXIT loop
            END DO
            DO i = 1,4
                previous(i+1) = previous(i)
            END DO
            previous(1) = current
        END DO loop
        PRINT *, current**3, current**3.0_DP
    END DO mainloop

CONTAINS

FUNCTION value (arg, targ)
    USE double
    IMPLICIT NONE
    REAL(KIND=DP) :: value, arg, targ
    value = arg**3-targ
END FUNCTION value

FUNCTION derivative (arg)
    USE double
    IMPLICIT NONE
    REAL(KIND=DP) :: derivative, arg
    derivative = 3*arg**2
END FUNCTION derivative

END PROGRAM Newton