```{-|

Module      :  Data.Number.ER.Real.Arithmetic.Taylor
Description :  interval Newton method
Copyright   :  (c) Amin Farjudian, Michal Konecny

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

Interval Newton's method for root finding.

To be used for obtaining functions out of their inverse(s) over various
intervals.
-}
module Data.Number.ER.Real.Arithmetic.Newton
(
erNewton_FullArgs,
erNewton_mdfd_FullArgs
)
where

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

erNewton_FullArgs
:: (RA.ERIntApprox ira)
=> (EffortIndex -> ira -> ira, EffortIndex -> ira -> ira) -- ^ a function and its derivative
-> ira -- ^ a starting point
-> ira -- ^ a lower bound of the absolute value of the derivative over the working interval
-> Int -- ^ number of iterations
-> EffortIndex  -- ^ the initial index to use for argument function and its derivative
-> ira -- ^ the result

erNewton_FullArgs (f ,df) startPt minDrv iterCnt i =
erNewton_FullArgs_aux startPt startOtherPt iterCnt
where
erNewton_FullArgs_aux newtonPt otherPt iterCnt
| (iterCnt <= 0 || RA.getPrecision result >= prec) =
result
| otherwise =
erNewton_FullArgs_aux newNewtonPt newOtherPt (iterCnt - 1)
where
result =
newtonPt RA.\/ otherPt
prec = effIx2prec i
newNewtonPt =
midPoint \$ RA.bounds \$
(newtonPt - ( (f i newtonPt) / (( df i newtonPt))))
--  /\ (ira2ra ((ra2ira minDrv) \/ 100000000)))))
newOtherPt = otherEndPoint newNewtonPt
startOtherPt = otherEndPoint startPt
otherEndPoint a =  a - ((f i a) / minDrv) --   /\ (0 \/ 10000000)

{-|
This auxiliary function returns the average of two ra's
-}
midPoint
:: (RA.ERIntApprox ira)
=> (ira ,ira)
-> ira
midPoint (x1, x2) = (x1 + x2) / 2

{-| Modified Newton Method
Notes:

1. It has a cubic convergence speed, as opposed to the original Newton's
square convergence speed.

2. It does not deal with multiple roots.

3. Per iteration, it makes two queries on the derivative, so it best
suits the cases where computation of the derivative is at most as
expensive as the function itself.
-}
erNewton_mdfd_FullArgs
:: (RA.ERIntApprox ira)
=> (EffortIndex -> ira -> ira, EffortIndex -> ira -> ira) -- ^ a function and its derivative
-> ira -- ^ a starting point
-> ira -- ^ The minimum of absolute value of derivative over the working interval
-> Int -- ^ number of iterations
-> EffortIndex  -- ^ It triggers the initial index to be called by the argument function and its derivative.
-> ira -- ^ the result

erNewton_mdfd_FullArgs (f ,df) startPt minDrv iterCnt i =
erNewton_FullArgs_aux startPt startOtherPt iterCnt
where
erNewton_FullArgs_aux newtonPt otherPt iterCnt
| iterCnt <= 0 = newtonPt RA.\/ otherPt
| otherwise = erNewton_FullArgs_aux newNewtonPt newOtherPt (iterCnt - 1)
where
valueAtNewOtherPt = f i newOtherPt
derivAtNewtonPt   = df i newOtherPt
unblurredDeriv = midPoint \$ RA.bounds \$ derivAtNewtonPt
intermediaryPt = midPoint \$ RA.bounds \$ newtonPt - valueAtNewOtherPt / (2 * derivAtNewtonPt)
derivAtIntermediaryPt = df i intermediaryPt
newNewtonPt =
midPoint \$ RA.bounds \$
(newtonPt - ( valueAtNewOtherPt / derivAtIntermediaryPt))
newOtherPt = otherEndPoint newNewtonPt
startOtherPt = otherEndPoint startPt
otherEndPoint a = a - ((f i a) / minDrv)

erNewton_mdfd
:: (RA.ERIntApprox ira)
=> (EffortIndex -> ira -> ira, EffortIndex -> ira -> ira) -- ^ a function and its derivative
-> ira -- ^ a starting point
-> ira -- ^ The minimum of absolute value of derivative over the working interval
-> EffortIndex  -- ^ It triggers the initial index to be called by the argument function and its derivative.
-> ira -- ^ the result

erNewton_mdfd (f ,df) startPt minDrv i =
erNewton_mdfd_FullArgs (f, df) startPt minDrv (fromInteger \$ toInteger \$ i) i

--apNewton_mdfd
--    :: (RA.ERIntApprox ira)
--    => (EffortIndex -> ra -> ra, EffortIndex -> ra -> ra) -- ^ a function and its derivative
--    -> ra -- ^ a starting point
--    -> ra -- ^ The minimum of absolute value of derivative over the working interval
--    -> EffortIndex  -- ^ It triggers the initial index to be called by the argument function and its derivative. Moreover, the number of iterations are predefined by this argument.
--    -> ra -- ^ the result
--
--apNewton_mdfd (f, df) startPt minDrv i =
--    erNewton_mdfd_FullArgs
--

--id_RA
--	:: (RA.ERIntApprox ira)
--	=> EffortIndex -> ira -> ira
--
--id_RA i x = x
--
--const_one_RA
--	:: (RA.ERIntApprox ira)
--	=> EffortIndex -> ira -> ira
--
--const_one_RA i x = (setMinGranularity (effIx2gran i) 1)
--
--
--test_erNewton_FullArgs_01_RA
--	:: (RA.ERIntApprox ira)
--	=> EffortIndex -> ira -> ira
--
--test_erNewton_FullArgs_01_RA i x = erNewton_FullArgs_01 ( id_RA, const_one_RA) x 10 i
--
--test_erNewton_FullArgs_01
--	:: (RA.ERIntApprox ira)
--	=> (ConvergRealSeq ira) -> (ConvergRealSeq ira)
--
--test_erNewton_FullArgs_01 = convertFuncRA2Seq test_erNewton_FullArgs_01_RA
--
--exp_Ra_minus_r_RA
--	:: (RA.ERIntApprox ira)
--	=> EffortIndex -> ira -> ira -> ira
--
--exp_Ra_minus_r_RA i r x = (erExp_RA i x) - r
--
--exp_Ra_minus_r
--	:: (RA.ERIntApprox ira)
--	=> (ConvergRealSeq ira) -> (ConvergRealSeq ira) -> (ConvergRealSeq ira)
--
--exp_Ra_minus_r = convertBinFuncRA2Seq exp_Ra_minus_r_RA
--
--logNewton_RA_02
--    :: (RA.ERIntApprox ira)
--    => EffortIndex -> ira -> ira
--
--logNewton_RA_02 i x =
--    erNewton_FullArgs_02
--        ( \ i y -> (erExp_RA i y) - x, erExp_RA)
--        (setMinGranularity (effIx2gran i) (2))
--        (setMinGranularity (effIx2gran i) 1)
--        i
--
--logNewton_02
--    :: (RA.ERIntApprox ira)
--    => (ConvergRealSeq ira) -> (ConvergRealSeq ira)
--
--logNewton_02 = convertFuncRA2Seq logNewton_RA_02

--logNewton_mdfd_RA
--    :: (RA.ERIntApprox ira)
--    => EffortIndex -> ira -> ira

--logNewton_mdfd_RA i x =
--    erNewton_mdfd_FullArgs
--        ( \ i y -> (erExp_RA i y) - x, erExp_RA)
--        (setMinGranularity (effIx2gran i) (2))
--        (setMinGranularity (effIx2gran i) 1)
--        i

--logNewton_mdfd
--    :: (RA.ERIntApprox ira)
--    => (ConvergRealSeq ira) -> (ConvergRealSeq ira)
--
--logNewton_mdfd = convertFuncRA2Seq logNewton_mdfd_RA
```