module functions
    implicit none
    contains
    
    function f(x) result(y)
        real(8), intent(in) :: x
        real(8) :: y
        y = ((((x * x) - 3.5d0) * x + 2.5d0) * x - 7d0) * x - 6.4d0
    end function
    
    
    function F1(x) result(y)
        real(8), intent(in) :: x
        real(8) :: y

!        y = (1d0/6d0)*(x**6) - (3.5d0/4d0)*(x**4) + (2.5d0/3d0)*(x**3) - (7d0/2d0)*(x**2) - 6.4d0*x
        y = (((((1d0/6d0) * x * x -  3.5d0/4d0) * x + 2.5d0/3d0) * x - 7d0/2d0) * x - 6.4d0) * x
    end function
    
    
    function trap(grid, N, h) result(int1)
        integer, intent(in) :: N
        integer :: i
        real(8), dimension(:) :: grid(N+1)
        real(8) :: h, int1
        int1 = 0d0
        do i = 2, N
            int1 = int1 + grid(i)
        end do
        int1 = (2d0*int1 + grid(1) + grid(N+1))*h/2d0
    end function
    
    
    subroutine intrect(a, b, eps, err, h, N)
        integer :: N, i
        real(8) :: a, b, eps, err, h
        real(8) :: int1, prevint, k, r
        prevint = 0d0
        int1 = 100d0
        N = 1
        do while (dabs(int1 - prevint)/3d0 > (eps))
            N = N * 2
            h = (b - a) / N
            prevint = int1
            int1 = 0
            do i = 0, N - 1
                int1 = int1 + f(a + (i + 0.5d0) * h)
            end do
            int1 = int1 * h
        end do
        err = dabs((F1(b) - F1(a)) - int1)
        write(*,*) N, eps, err
    end subroutine

    subroutine inttrap(a, b, eps, err, h, N)
        integer :: N, i
        real(8) :: a, b, eps, err, h
        real(8) :: int1, prevint
        prevint = 0d0
        int1 = 100d0
        N = 1
        do while (dabs(int1 - prevint) / 3 > (eps))
            N = N * 2
            h = (b - a) / N
            prevint = int1
            int1 = 0.5d0 * (f(a) + f(b))
            do i = 1, N - 1
                int1 = int1 + f(a + i * h)
            end do
            int1 = int1 * h
        end do
        err = dabs((F1(b) - F1(a)) - int1)
        write(*,*) N, eps, err
    end subroutine

    subroutine intsimps(a, b, eps, err, h, N)
        integer :: N, i
        real(8) :: a, b, eps, err, h
        real(8) :: int1, inta, intb, prevint
        prevint = 0d0
        int1 = 100d0
        N = 1
        do while (dabs(int1 - prevint) / 15 > eps)
            N = N * 2
            h = (b - a) / N

            inta = 0.5d0 * (f(a) + f(b))
            do i = 1, N - 1
                inta = inta + f(a + i * h)
            end do
            intb = 0
            do i = 0, N - 1
                intb = intb + f(a + (i + 0.5d0) * h)
            end do

            prevint = int1
            int1 = ( inta + 2 * intb) * h / 3.0d0

        end do

        err = dabs((F1(b) - F1(a)) - int1)
        write(*,*) N, eps, err
    end subroutine



end module