!
! NR_polinomial.f90
!
! compiled: gfortran NR_polinomial.f90 -o b -fdefault-double-8 -fdefault-real-8
!           ifort NR_polinomial.f90 -o b -no-complex-limited-range -real-size 64
!
! Autor: Pedro Garcia Freitas <sawp@sawp.com.br>
! License: Creative Commons
!          <http://creativecommons.org/licenses/by-nc-nd/2.5/br/>
!
program FP
	implicit none;
	complex, external :: polinomial, diffpolinomial;
	complex, external :: newtonRaphson_c;
	complex result;

	real :: x=1.2d0, errto=0.0001d0;
	integer :: imax=1000;

	result = newtonRaphson_c(polinomial, diffpolinomial, x, errto, imax);

	print *, "The root of f(x) = -40. - 2.*x + 27.*x**2 - 10.*x**3 + & 
                        & x**4 : ", result; 
end program FP

! -
! newtonRaphson_c
! return the [complex] root of a function (using the Newton-Raphson Method)
! G: function pointer
! x: some initial point
! imax: max of iterations allowed
! errto: toleratec error 
!
! Author: Pedro Garcia sawp@sawp.com.br
! see: http:!www.sawp.com.br
!
! License: Creative Commons
!          <http://creativecommons.org/licenses/by-nc-nd/2.5/br/>
! ago 2009
! 
complex function newtonRaphson_c(G, diffG, x, errto, imax)
	implicit none;
	complex, external :: G, diffG;
	real, intent(in) :: x, errto;
	integer :: imax;

	! iteration counter ! 
	integer :: iterCount;

	! relative error ! 
	real(kind=8) :: errno;

	! xi+1 ! 
	complex :: x1,x0;

	x0 = (0.,0.);
	x1 = cmplx(x, 0.d0);
	errno = 100.d0;
	iterCount = 0;


	do
		x0 = x1;
		x1 = x0 - G(x0)/diffG(x0);
		iterCount = iterCount + 1;

		if(x1 /= (0.d0,0.d0) )then			
			errno = abs( (cabs(x1) - cabs(x0))/cabs(x1) ) * 100.d0
		end if
			
		if(errno < errto .or. iterCount > imax)then
			exit;
		end if	
	end do

	newtonRaphson_c = x1;
	return;
end function newtonRaphson_c

! - 
! Test Function
complex function polinomial(x)
	implicit none;
	complex, intent(in) :: x;

	polinomial = (-40. - 2.*x + 27.*x**2 - 10.*x**3 + x**4);
	return;
end function polinomial

complex function diffpolinomial(x)
	implicit none;
	complex, intent(in) :: x;

	diffpolinomial = (-2. + 54.*x - 30.*x**2 + 4.*x**3);
end function diffpolinomial