{-# LANGUAGE RebindableSyntax #-}

module Integration.Trapezoid
where

import Misc (partition)

import NumericPrelude hiding (abs)
import Algebra.Absolute (abs)
import qualified Algebra.Field as Field
import qualified Algebra.RealField as RealField
import qualified Algebra.RealRing as RealRing
import qualified Algebra.ToInteger as ToInteger
import qualified Algebra.ToRational as ToRational

-- | 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 :: (Field.C a, ToRational.C a, Field.C 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 rule 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 0 pi
--   >>> abs (area - 2) < 0.0001
--   True
--
trapezoid :: (RealField.C a,
              ToRational.C a,
              RealField.C b,
              ToInteger.C c,
              Enum 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)