-- | Compute definite integrals.
module Numerics.Integrate (
      integrateMidpoint
    , integrateSimpson13
) where

import Data.List

sum' :: Num a => [a] -> a
sum' :: forall a. Num a => [a] -> a
sum' = (a -> a -> a) -> a -> [a] -> a
forall b a. (b -> a -> b) -> b -> [a] -> b
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' a -> a -> a
forall a. Num a => a -> a -> a
(+) a
0
{-# INLINE sum' #-}

-- | Numerical integration with the midpoint method. This is one of the simplest
-- integration algorithms.
--
-- https://en.wikipedia.org/wiki/Midpoint_method
integrateMidpoint
    :: Fractional a
    => (a -> a) -- ^ \(f\)
    -> a        -- ^ \(a\)
    -> a        -- ^ \(b\)
    -> Int      -- ^ Number of interval subdivisions
    -> a        -- ^ \(\int_a^b\mathrm dt\, f(t)\)
integrateMidpoint :: forall a. Fractional a => (a -> a) -> a -> a -> Int -> a
integrateMidpoint a -> a
f a
a a
b Int
n =
    a
h a -> a -> a
forall a. Num a => a -> a -> a
* ( a -> a
f a
a a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
2
        a -> a -> a
forall a. Num a => a -> a -> a
+ [a] -> a
forall a. Num a => [a] -> a
sum' [a -> a
f (a
a a -> a -> a
forall a. Num a => a -> a -> a
+ Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
ka -> a -> a
forall a. Num a => a -> a -> a
*a
h) | Int
k <- [Int
1..Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]]
        a -> a -> a
forall a. Num a => a -> a -> a
+ a -> a
f a
b a -> a -> a
forall a. Fractional a => a -> a -> a
/ a
2
        )
  where
    h :: a
h = (a
ba -> a -> a
forall a. Num a => a -> a -> a
-a
a)a -> a -> a
forall a. Fractional a => a -> a -> a
/Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n

-- | Numerical integration with Simpson’s ⅓ rule, which approximates curves with
-- quadratic polynomials. It’s not that complicated but still does very well on
-- not-so-crazy functions.
--
-- https://en.wikipedia.org/wiki/Simpson%27s_rule
--
-- Useful as a parameter to 'Numerics.ConvergentRecursion.recurseUntilPrecision'.
integrateSimpson13
    :: Fractional a
    => (a -> a) -- ^ \(f\)
    -> a        -- ^ \(a\)
    -> a        -- ^ \(b\)
    -> Int      -- ^ Number of interval subdivisions
    -> a        -- ^ \(\int_a^b\mathrm dt\, f(t)\)
integrateSimpson13 :: forall a. Fractional a => (a -> a) -> a -> a -> Int -> a
integrateSimpson13 a -> a
f a
a a
b Int
n
    | Int -> Bool
forall a. Integral a => a -> Bool
odd Int
n = (a -> a) -> a -> a -> Int -> a
forall a. Fractional a => (a -> a) -> a -> a -> Int -> a
integrateSimpson13 a -> a
f a
a a
b (Int -> Int
forall a. Enum a => a -> a
succ Int
n)
integrateSimpson13 a -> a
f a
a a
b Int
n =
    a
ha -> a -> a
forall a. Fractional a => a -> a -> a
/a
3 a -> a -> a
forall a. Num a => a -> a -> a
* ( a -> a
f a
a
          a -> a -> a
forall a. Num a => a -> a -> a
+ a
2 a -> a -> a
forall a. Num a => a -> a -> a
* [a] -> a
forall a. Num a => [a] -> a
sum' [a -> a
f (Int -> a
x (Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
j  )) | Int
j <- [Int
1..Int
nInt -> Int -> Int
forall a. Integral a => a -> a -> a
`div`Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1]]
          a -> a -> a
forall a. Num a => a -> a -> a
+ a
4 a -> a -> a
forall a. Num a => a -> a -> a
* [a] -> a
forall a. Num a => [a] -> a
sum' [a -> a
f (Int -> a
x (Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
jInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)) | Int
j <- [Int
1..Int
nInt -> Int -> Int
forall a. Integral a => a -> a -> a
`div`Int
2]]
          a -> a -> a
forall a. Num a => a -> a -> a
+ a -> a
f a
b
          )
  where
    x :: Int -> a
x Int
i = a
a a -> a -> a
forall a. Num a => a -> a -> a
+ Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
i :: Int)a -> a -> a
forall a. Num a => a -> a -> a
*a
h
    h :: a
h = (a
ba -> a -> a
forall a. Num a => a -> a -> a
-a
a)a -> a -> a
forall a. Fractional a => a -> a -> a
/Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n