{-# LANGUAGE NoImplicitPrelude #-}
{-# LANGUAGE RebindableSyntax #-}

module Integration.Simpson (
  simpson,
  simpson_1 )
where

import Misc ( partition )

import NumericPrelude hiding ( abs )
import qualified Algebra.RealField as RealField ( C )
import qualified Algebra.ToInteger as ToInteger ( C )
import qualified Algebra.ToRational as ToRational ( C )

-- | Use the Simpson's rule to numerically integrate @f@ over the
--   interval [@a@, @b@].
--
--   Examples:
--
--   >>> let f x = 1
--   >>> simpson_1 f (-1) 1
--   2.0
--
--   >>> let f x = x
--   >>> simpson_1 f (-1) 1
--   0.0
--
--   >>> import Algebra.Absolute (abs)
--   >>> let f x = x^2
--   >>> let area = simpson_1 f (-1) 1
--   >>> abs (area - (2/3)) < 1/10^12
--   True
--
--   >>> let f x = x^3
--   >>> simpson_1 f (-1) 1
--   0.0
--
--   >>> let f x = x^3
--   >>> simpson_1 f 0 1
--   0.25
--
simpson_1 :: (RealField.C a, ToRational.C a, RealField.C b)
            => (a -> b) -- ^ The function @f@
            -> a       -- ^ The \"left\" endpoint, @a@
            -> a       -- ^ The \"right\" endpoint, @b@
            -> b
simpson_1 f a b =
  coefficient * ((f a) + 4*(f midpoint) + (f b))
  where
    coefficient = fromRational' $ (toRational (b - a)) / 6
    midpoint = (a + b) / 2


-- | Use the composite Simpson's rule to numerically integrate @f@
--   over @n@ subintervals of [@a@, @b@].
--
--   Examples:
--
--   >>> import Algebra.Absolute (abs)
--   >>> let f x = x^4
--   >>> let area = simpson 10 f (-1) 1
--   >>> abs (area - (2/5)) < 0.0001
--   True
--
--   Note that the convergence here is much faster than the Trapezoid
--   rule!
--
--   >>> import Algebra.Absolute (abs)
--   >>> let area = simpson 10 sin 0 pi
--   >>> abs (area - 2) < 0.00001
--   True
--
simpson :: (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
simpson n f a b =
  sum $ map simpson_pairs pieces
  where
    pieces = partition n a b
    -- Convert the simpson_1 function into one that takes pairs
    -- (a,b) instead of individual arguments 'a' and 'b'.
    simpson_pairs = uncurry (simpson_1 f)