program prak

real :: c, d, h, sigma, eps=0.0001, lmdb, sm
integer :: N, n0, m=2, i, k, j
real, allocatable, dimension (:) :: X, Y, polinom, gradient, V, tmp
logical :: flag

    write (*, '(A)', ADVANCE='NO') 'c:  '
    read *, c
    write (*, '(A)', ADVANCE='NO') 'd:  '
    read *, d
    write (*, '(A)', ADVANCE='NO') 'N:  '
    read *, N
    write (*, '(A)', ADVANCE='NO') 'no:  '
    read *, n0
    print *,

    h=(d-c)/N
    allocate (X(N+1))

    do i=1, N+1
        X(i)=c+(i-1)*h
    end do

    allocate (Y(N+1))
    write (*, '(A)') 'function:'
    write (*, '(A)') 'test 1) y = 6*x^3-0.533*x^2+0.88*x+2.4'
    write (*, '(A)') 'prod 2) y = (1+x)*exp(-2*x)'
    read *, i
    print *,

    select case (i)
        case(1)
            call func_test (X, Y, N+1)
        case(2)
            call func_1 (X, Y, N+1)
    end select

    allocate (polinom(n0+m))
    allocate (gradient(n0+m))
    allocate (V(n0+m))
    allocate (tmp(n0+m))



    do k=n0, n0+m-1
        i=0
        flag=.true.
        sigma=1000
    do j=1, k+1
        polinom(j)=0
        gradient(j)=0
        V(j)=0
        tmp(j)=0
    end do

    write (*, '(A15, i2)') 'n: ', k

        do while (sigma>eps)
            tmp=polinom
            i=i+1

            call get_V()

            lmdb = lambda()
            polinom=polinom+lmdb*V
            if (flag.eqv..true.) then
                flag=.false.
            else
                sigma=get_sigma()
            end if

            sm=0
            do l=1, N+1
                sm = sm + (Y(l)-P(X(l), polinom))**2
            end do
            sm = sqrt(sm/(N+1))

            write(*,'(I6, A12, F14.8, A12, F14.8)') i, ': |grad(F)|=', sqrt(dot_product(gradient,gradient)), 'dev: ', sm

        end do

        print *,

        write (*, '(A2, A15, A15, A15, A15)') 'it', 'X(i)   ', 'Y(i)    ','P(i)    ','d    '

        do i=1, N+1
            write (*, '(I2, F14.8, F14.8, F14.8, F14.8)') i-1, X(i), Y(i), P(X(i), polinom), abs(Y(i)-P(X(i), polinom))
        end do

        print *,
        write (*, '(A)') 'a0..an'
        do i=0, k
            write (*, '(A1, I2, A1, F12.6)') 'a',i, ' =', polinom(i+1)
        end do
        print *,

   end do

   deallocate (X)
   deallocate (Y)
   deallocate (polinom)
   deallocate (V)
   deallocate (gradient)
   deallocate (tmp)

   read *,

contains

    real function P (point, array)

        real :: point
        integer :: i
        real, dimension (:) :: array

        P=0

        do i=0, k
            P= array(k+1-i)+P*point
        end do

    end function P

    subroutine get_gradient ()
        integer :: i, j

        do j=1, k+1
            gradient(j)=0
            do i=1, N+1
                gradient(j)=gradient(j)+2*(X(i)**(j-1))*(P(X(i), polinom)-Y(i))
            end do
        end do

    end subroutine get_gradient

    real function G (array)
        real, dimension (:) :: array
        integer :: i

        G=0

        do i=1, N+1
            G=G+(Y(i)-P(X(i), array))**2
        end do

    end function G

    subroutine get_V ()
        real :: tmp

        if (flag.eqv..TRUE.) then
                call get_gradient ()
                V=(-1)*gradient
            else
                tmp=dot_product(gradient, gradient)
                call get_gradient ()
                V=(-1)*gradient+(dot_product(gradient, gradient)/tmp)*V
        end if

    end subroutine get_V

    real function lambda ()
        real :: m1, m2, m3

        m1=G(polinom-V)
        m2=G(polinom)
        m3=G(polinom+V)

        lambda=(m1-m3)/(2*(m1-2*m2+m3))

    end function lambda

    real function get_sigma ()
        integer :: i
        real :: q

        tmp=tmp-polinom
        get_sigma=0

        do i=1, k+1
            if (polinom(i)/=0) then
                q=abs(tmp(i)/polinom(i))
                if (q>get_sigma) then
                    get_sigma=q
                end if
            end if
        end do

    end function get_sigma

end program prak

subroutine func_1 (X, Y, N)

    integer :: N, i
    real, dimension (N) :: X, Y

    do i=1, N
        Y(i)=(1+X(i))*exp(-2*X(i))
    end do

end subroutine func_1

subroutine func_test (X, Y, N)

    integer :: N, i
    real, dimension (N) :: X, Y

    do i=1, N
        Y(i)=6*X(i)**3-0.533*X(i)**2+0.88*X(i)+2.4
    end do

end subroutine func_test