{-# LANGUAGE RebindableSyntax #-}

module Integration.Simpson
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 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
--
--   >>> 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:
--
--   >>> 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!
--
--   >>> 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)