!       ニュートン法(減速法)
!       f(x)=x**5.0e0-7.0e0*X**4.0e0+16.0e0*x**3.0e0-8.0e0*x**2.0e0-16.0e0*x+16.0e0
!       f(x)=0の解
!       伊理テキストP64の表9-1(右)の作成
!       作成者    福士 邦彦
!       2011/06/07 新規作成
!       2011/06/09 最終更新
!       2012/07/12 増田 和孝 最終更新

!---------------メインプログラム-----------------------

program gensoku_newton1

    implicit none
    integer :: loop     ! ループ回数のカウンター
    real :: myu      ! 減速法のmyuの値
    real :: kekka    ! 減速法の終了判定に使う値
    real :: xs       ! xの初期値
    real :: xe       ! xの計算値
    real :: lim      ! 計算終了判定の基準値
    real :: hantei   ! 計算終了判定の判定値
    real :: ff       ! 関数f(x)
    real :: dfdx     ! f'(x)の値


    ! xの初期値を代入
    write(*,*) '最初のxの値を入力してください >'
    read(*,*) xs


    !---------解を求める計算------------------------

        
    !write(*,*) dfdx(xs)
    ! 計算終了判定の基準値を設定
    lim=epsilon(0.0e0)

    ! ループカウンターの初期値を設定
    loop=0

    !xeの初期値を設定
    xe=0.0e0

    ! ループに入るためにhanteiに適当な値を代入
    hantei=abs((xs-xe)/xs)

    ! abs((xs-xe)/xs)<limとなるまで計算を続けるループ
    do while (hantei>lim)

        ! ループカウンターの更新
        loop=loop+1

        !myuの初期値を設定
        myu=1.0e0

        !ループに入るためのkekkaの値
        kekka=abs(ff(xs-myu*ff(xs)/dfdx(xs)))/abs(ff(xs))
        !write(*,*) myu, kekka
        
        ! ------myuの値を決定するループ-----------------
        do while (kekka>=(1.0e0-myu/4.0e0))

           ! myuの値を更新
           myu=myu/2.0e0

           !kekkaの値を更新
           kekka=abs(ff(xs-myu*ff(xs)/dfdx(xs)))/abs(ff(xs))
        !   write(*,*) myu, kekka
        end do
        !ループ終了---------------------------------
    
        ! xeを計算
        xe=xs-myu*ff(xs)/dfdx(xs)

        ! 判定値を更新
        hantei=abs((xs-xe)/xs)

        ! xsを更新
        xs=xe

        ! ループ回数とxsとff(xs)の値を表示
        write(*,*) loop,xs,ff(xs),myu

    end do

end program gensoku_newton1



!-------------function　f(x)の値の計算--関数名ff---------------

function ff(x)

    implicit none
    real,intent(in) :: x        ! xの値
    real :: ff       ! f(x)の値    

    ! f(x)の計算
    ff=x**5.0e0-7.0e0*X**4.0e0+16.0e0*x**3.0e0-8.0e0*x**2.0e0-16.0e0*x+16.0e0

    return

end function ff





!-------------function　f'(x)の値の計算--関数名dfdx---------------

function dfdx(x)

    implicit none
    real,intent(in) :: x        ! xの値
    real :: dfdx     ! f'(x)の値

    ! f'(x)の計算
    dfdx=5.0e0*x**4.0e0-28.0e0*x**3.0e0+48.0e0*x**2.0e0-16.0e0*x-16.0e0

    return

end function dfdx
