```{-|
Module      :  Data.Number.ER.Real.Arithmetic.Integration
Description :  simple integration methods

Maintainer  :  mik@konecny.aow.cz
Stability   :  experimental
Portability :  portable

Simple integration methods for Haskell functions operating
on real number approximations.
-}
module Data.Number.ER.Real.Arithmetic.Integration
(
integrateCont,
--    integrateDiff,
integrateCont_R,
)
where

import qualified Data.Number.ER.Real.Approx as RA
import Data.Number.ER.BasicTypes
import Data.Number.ER.Real.Approx.Sequence
import Data.Number.ER.Real.Arithmetic.Elementary

testIntegr1 ::
(RA.ERIntApprox ira) =>
(ConvergRealSeq ira)
testIntegr1 = integrateCont erExp_IR 0 1

integrateCont ::
(RA.ERIntApprox ira) =>
(EffortIndex -> ira -> ira) ->
(ConvergRealSeq ira) -> (ConvergRealSeq ira) -> (ConvergRealSeq ira)

integrateCont f = convertBinFuncRA2Seq \$ integrateContAdapt_R f

integrateDiff ::
(RA.ERIntApprox ira) =>
(EffortIndex -> ira -> ira, EffortIndex -> ira -> ira) ->
(ConvergRealSeq ra) -> (ConvergRealSeq ra) -> (ConvergRealSeq ra)

integrateDiff f = convertBinFuncRA2Seq \$ integrateDiffAdapt_RA f

{-|
naive integration, using a partition of 2 * prec equally sized intervals
-}
integrateCont_R ::
(RA.ERIntApprox ira) =>
(EffortIndex -> ira -> ira) ->
EffortIndex -> (ira) -> (ira) -> (ira)
integrateCont_R f ix a b =
sum \$ map rectArea rectangles
where
rectArea (width, height) = width * height
rectangles =
zip (repeat width) \$ map (f ix) covering
width = (b - a) / (fromInteger rectCount)
rectCount = 2 * (fromInteger \$ toInteger gran)
gran = effIx2gran ix
covering = getCoveringIntervals division
getCoveringIntervals ( pt1 : pt2 : rest ) =
((pt1) RA.\/ (pt2)) : (getCoveringIntervals \$ pt2 : rest)
getCoveringIntervals _ = []
division = map getEndPoint \$ [0..rectCount]
getEndPoint n =
a + ((fromInteger n) * width)

{-|
integration using divide and conquer adaptive partitioning
-}
(RA.ERIntApprox ira) =>
(EffortIndex -> ira -> ira) ->
EffortIndex -> (ira) -> (ira) -> (ira)
integrateContAdapt_R f ix a b =
sum rectangleAreas
where
rectangleAreas =
getRs ix a b
getRs subix l r
| RA.getPrecision area >= prec = [area]
| otherwise =
(getRs nsubix l m) ++ (getRs nsubix m r)
where
prec = foldl1 min [effIx2prec ix, RA.getPrecision l, RA.getPrecision r]
area = (r - l) * (f subix (l RA.\/ r))
nsubix = subix + 1
m = (l + r)/2

{-|
integration using divide and conquer adaptive partitioning
making use of the derivative of the integrated function
-}
(RA.ERIntApprox ira) =>
(EffortIndex -> ira -> ira, EffortIndex -> ira -> ira) ->
EffortIndex -> (ra) -> (ra) -> (ra)
integrateDiffAdapt_RA f prec a b =
error "TODO"

{-
sum rectangleAreas
where
rectangleAreas =
getRs prec a b
getRs p l r
| getPrecision area >= prec = [area]
| otherwise =
(getRs np l m) ++ (getRs np m r)
where
np = p + 1
m = (l + r)/2
-- merge the information given by the rectangle method
-- with the information given by the derivative method
areaRect = w * fVal -- same as in integrateContAdapt_R
(fVal, fDeriv) = applyRdiffR f p (l \/ r)
w = r - l
| isExact fDeriv = w * (fl + fr) / 2 -- derivative is constant and perfectly known
| otherwise = areaLow \/ areaHigh
fl = fst \$ applyRdiffR f (2 * p) l
fr = fst \$ applyRdiffR f (2 * p) r
-- interestingly, we have to request fl, fr with higher precision than
-- we requested fDeriv so that the derivative would be of any use
-- with these values - replace (2 * p) by p and it will not converge!
-- area computed by a scary formula:
areaLow = t + w * (fl * dHigh - fr * dLow) / dDiff
areaHigh = - t - w * (fl * dLow - fr * dHigh) / dDiff -- swap dHigh and dLow
t = (w^2 * dLow * dHigh + (fr - fl)^2)/(2 * dDiff)
dDiff = dHigh - dLow
(dLow, dHigh) = getBounds fDeriv
-}

```