!********o*********o*********o*********o*********o*********o*********o**
! This program, given a 5th degree polynomial, attempts to find a real
! root.
!
! Method
!   The method used is Newton's method.  Given a point x at which f(x)
!   is non-zero, we speculate that f is almost a straight line at x,
!   and we find where that straight line (which is the tangent to f at
!   x) will be zero.  Mathematically, we have
!      new x = x - f(x)/f'(x)
!   (where f' is the derivative of f).  We do this repeatedly until we
!   find an x for which f(x) is close enough to zero.
!
! For more Information, see
!   www.ugrad.math.ubc.ca/coursedoc/math100/notes/approx/newton.html
!   www.cs.tcd.ie/John.Byrne/Roots/newton.html
!---------------------
! Inputs
!   c5,c4,c3,c2,c1,c0: The coefficients of the polynomial
!                      i.e. f(x) = c5*x**5 + c4*x**4 + ... + c1*x + c0
!
! Outputs
!   x:  A value for which f(x) = 0
!
! Constants
!   goodEnough: A small number to control how close we need to get to a
!               solution.
!
! Other Variables
!   d4,d3,d2,d1,d0: The coefficients of the derivative of the polynomial
!   fx:             The value of f, at x
!   dx:             The value of the derivative of f, at x
!---------------------
! Author:   Bob Harris
! Class:    CMPSC 201F
! Section:  1234
! Date Due: 8 Oct 2003
!********o*********o*********o*********o*********o*********o*********o**

PROGRAM Newton

IMPLICIT NONE

REAL, PARAMETER :: goodEnough = 0.000001
REAL :: c5,c4,c3,c2,c1,c0
REAL ::    d4,d3,d2,d1,d0
REAL :: x,fx,dx

! fetch the polynomial from the user

PRINT *, "Please give me the 6 coefficients, from x**5 down to x**0"
READ  *, c5,c4,c3,c2,c1,c0

! compute the coefficients of the derivative

d4 = 5*c5
d3 = 4*c4
d2 = 3*c3
d1 = 2*c2
d0 =   c1

! show the user the polynomial and derivative

PRINT *, "f(x) = ",c5,"x**5+",c4,"x**4+",c3,"x**3+",c2,"x**2+",c1,"x+",c0
PRINT *, "d(x) = ",d4,"x**4+",d3,"x**3+",d2,"x**2+",d1,"x+",d0

! pick and initial value

x = 0

! find the root, using newton's method

DO
    ! compute the value of f(x) and f'(x)

    fx = ((((c5*x+c4)*x+c3)*x+c2)*x+c1)*x+c0
    dx =  (((d4*x+d3)*x+d2)*x+d1)*x+d0

    ! see if this is close enough to being a root
    ! (shortcoming 1 is that this statement might never cause an exit)

    if (abs(fx) < goodEnough) EXIT

    ! report progress to our user

    PRINT *, "f(",x,") = ",fx

    ! make a new approximation, according to
    !   new x = x - f(x(n))/f'(x(n))
    ! (shortcoming 2 is this statement might cause a divide by zero)

    x = x - fx/dx

END DO

! tell the user our answer

PRINT *, "f(",x,") = ",fx

END PROGRAM Newton

