module Integration.Trapezoid
where


-- | Partition the interval [@a@, @b@] into @n@ subintervals, which we
--   then return as a list of pairs.
partition :: (RealFrac a, Integral b)
          => b -- ^ The number of subintervals to use, @n@
          -> a -- ^ The \"left\" endpoint of the interval, @a@
          -> a -- ^ The \"right\" endpoint of the interval, @b@
          -> [(a,a)]
 -- Somebody asked for zero subintervals? Ok.
partition 0 _ _ = []
partition n a b
  | n < 0 = error "partition: asked for a negative number of subintervals"
  | otherwise =
    [ (xi, xj) | k <- [0..n-1],
                 let k' = fromIntegral k,
                 let xi = a + k'*h,
                 let xj = a + (k'+1)*h ]
    where
      h = fromRational $ (toRational (b-a))/(toRational n)


-- | Use the trapezoid rule to numerically integrate @f@ over the
--   interval [@a@, @b@].
--
--   Examples:
--
--   >>> let f x = x
--   >>> trapezoid_1 f (-1) 1
--   0.0
--
--   >>> let f x = x^3
--   >>> trapezoid_1 f (-1) 1
--   0.0
--
--   >>> let f x = 1
--   >>> trapezoid_1 f (-1) 1
--   2.0
--
--   >>> let f x = x^2
--   >>> trapezoid_1 f (-1) 1
--   2.0
--
trapezoid_1 :: (RealFrac a, Fractional b, Num b)
            => (a -> b) -- ^ The function @f@
            -> a       -- ^ The \"left\" endpoint, @a@
            -> a       -- ^ The \"right\" endpoint, @b@
            -> b
trapezoid_1 f a b =
  (((f a) + (f b)) / 2) * (fromRational $ toRational (b - a))


-- | Use the composite trapezoid tule to numerically integrate @f@
--   over @n@ subintervals of [@a@, @b@].
--
--   Examples:
--
--   >>> let f x = x^2
--   >>> let area = trapezoid 1000 f (-1) 1
--   abs (area - (2/3)) < 0.00001
--   True
--
--   >>> let area = trapezoid 1000 sin (-1) 1
--   >>> abs (area - 2) < 0.00001
--   True
--
trapezoid :: (RealFrac a, Fractional b, Num b, Integral c)
          => c -- ^ The number of subintervals to use, @n@
          -> (a -> b) -- ^ The function @f@
          -> a       -- ^ The \"left\" endpoint, @a@
          -> a       -- ^ The \"right\" endpoint, @b@
          -> b
trapezoid n f a b =
  sum $ map trapezoid_pairs pieces
  where
    pieces = partition n a b
    -- Convert the trapezoid_1 function into one that takes pairs
    -- (a,b) instead of individual arguments 'a' and 'b'.
    trapezoid_pairs = uncurry (trapezoid_1 f)