MODULE functions

  real(kind=8), private, parameter :: poly_coeffs(11) =            &
       [-11.0, 15.6, -4.6, 22.5, 8.1, 5.1, -0.3, -8.0, 0.0, -9.9, 2.2]
  
CONTAINS

  FUNCTION my_exp(x) result(y)
    real(kind=8), intent(in) :: x
    real(kind=8) :: y
    y = exp(x)
  END FUNCTION my_exp

  FUNCTION my_sin(x) result(y)
    real(kind=8), intent(in) :: x
    real(kind=8) :: y
    y = sin(x)
  END FUNCTION my_sin

  FUNCTION my_poly(x) result(y)
    real(kind=8), intent(in) :: x
    real(kind=8) :: y
    integer(kind=4) :: i
    y = sum(poly_coeffs(:) * [1.0_8, (x ** [(i, i = 1, 10)])])
  END FUNCTION my_poly
  
  FUNCTION my_exp_int(ibeg, iend) result(y)
    real(kind=8), intent(in) :: ibeg, iend
    real(kind=8) :: y
    y = exp(iend) - exp(ibeg)
  END FUNCTION my_exp_int
  
  FUNCTION my_sin_int(ibeg, iend) result(y)
    real(kind=8), intent(in) :: ibeg, iend
    real(kind=8) :: y
    y = -cos(iend) + cos(ibeg)
  END FUNCTION my_sin_int

  FUNCTION my_poly_int_indefinite(x) result(y)
    real(kind=8), intent(in) :: x
    real(kind=8) :: y
    integer(kind=4) :: i, j
    y = sum(poly_coeffs(:) * (1 / real([(j, j = 1, 11)])) *        &
         (x ** [(i, i = 1, 11)]))
  END FUNCTION my_poly_int_indefinite
  
  FUNCTION my_poly_int(ibeg, iend) result(y)
    real(kind=8), intent(in) :: ibeg, iend
    real(kind=8) :: y
    y = my_poly_int_indefinite(iend) - my_poly_int_indefinite(ibeg)
  END FUNCTION my_poly_int
  
END MODULE functions