!
! muller.f90
!
! compiled: gfortran muller.f90 -o m -fdefault-double-8 -fdefault-real-8
!           ifort muller.f90 -o m -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;
        complex, external :: muller;
        complex result;

        real :: x=1.2d0;
        real :: errto=0.0001d0;
        integer :: imax=100;

        result = muller(polinomial, x, errto, imax);

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

! -
! muller
! return the [complex] root of a function (using the Muller Method)
!
! F: 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 muller(F, x, errto, imax)
        implicit none;
        complex, external :: F;
        real, intent(in) :: errto;
        real, intent(in) :: x;
        integer, intent(in) :: imax;

        ! iteration counter ! 
        integer :: iterCount;

        ! relative error ! 
        real :: dx;

        !
        complex :: a,b,c;
        complex :: x0,x1,x2,x3;
        complex h0, h1, d0, d1, delta, denominator;

        ! initialization block
        call random_seed();
        call random_number(dx);
        x2 = cmplx(x,0.);        
        x1 = cmplx(x + dx*x,0.);
        x0 = cmplx(x - dx*x, 0.);
        dx = 100.d0;
        iterCount = 0;

        do
                iterCount = iterCount + 1;
                h0 = x1 - x0;
                h1 = x2 - x1;
                d0 = (f(x1) - f(x2))/h0;
                d1 = (f(x2) - f(x1))/h1;

                a = (d1 - d0)/(h1+h0);
                b = a*h1 + d1;
                c = f(x2);
                
                delta = sqrt(b*b - 4.*a*c);

                if (abs(b+delta) > abs(b-delta)) then
                        denominator = b + delta;
                else
                        denominator = b - delta;
                end if

                dx = -2.*c/denominator;
                x3 = x2 + dx;

                if( (abs(dx) < errto*real(x) ) .or. iterCount >imax )then
                        exit;
                else
                        x0 = x1;
                        x1 = x2;
                        x2 = x3;
                end if
        end do
        
        muller = x3;

        return;
end function muller

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

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