{-# LANGUAGE CPP #-}
-- #define DEBUG

{-# OPTIONS_GHC -Wno-orphans #-}
{-|
    Module      :  AERN2.MP.Ball.Limit
    Description :  limits of MPBall sequences
    Copyright   :  (c) Michal Konecny
    License     :  BSD3

    Maintainer  :  mikkonecny@gmail.com
    Stability   :  experimental
    Portability :  portable

    Limits of MPBall sequences.
-}
module AERN2.MP.Ball.Limit where

#ifdef DEBUG
import Debug.Trace (trace)
#define maybeTrace trace
#define maybeTraceIO putStrLn
#else
#define maybeTrace (\ (_ :: String) t -> t)
#define maybeTraceIO (\ (_ :: String) -> return ())
#endif

import MixedTypesNumPrelude

import qualified Numeric.CollectErrors as CN

-- import Math.NumberTheory.Logarithms (integerLog2)

import AERN2.Limit
import AERN2.MP.Ball.Type
import AERN2.MP.Ball.Field ()

{-
  The following strategies are inspired by
  Mueller: The iRRAM: Exact Arithmetic in C++, Section 10.1
  https://link.springer.com/chapter/10.1007/3-540-45335-0_14
-}

instance HasLimits Rational (CN MPBall -> CN MPBall) where
  type LimitType Rational (CN MPBall -> CN MPBall) = (CN MPBall -> CN MPBall)
  limit :: (Rational -> CN MPBall -> CN MPBall)
-> LimitType Rational (CN MPBall -> CN MPBall)
limit = (Accuracy -> Rational)
-> (Rational -> CN MPBall -> CN MPBall) -> CN MPBall -> CN MPBall
forall ix.
(Accuracy -> ix)
-> (ix -> CN MPBall -> CN MPBall) -> CN MPBall -> CN MPBall
limitMPBall (\Accuracy
ac -> Rational
0.5Rational -> Integer -> PowType Rational Integer
forall t1 t2. CanPow t1 t2 => t1 -> t2 -> PowType t1 t2
^(Accuracy -> Integer
fromAccuracy Accuracy
ac))

instance HasLimits Integer (CN MPBall -> CN MPBall) where
  type LimitType Integer (CN MPBall -> CN MPBall) = (CN MPBall -> CN MPBall)
  limit :: (Integer -> CN MPBall -> CN MPBall)
-> LimitType Integer (CN MPBall -> CN MPBall)
limit = (Accuracy -> Integer)
-> (Integer -> CN MPBall -> CN MPBall) -> CN MPBall -> CN MPBall
forall ix.
(Accuracy -> ix)
-> (ix -> CN MPBall -> CN MPBall) -> CN MPBall -> CN MPBall
limitMPBall (\Accuracy
ac -> (Accuracy -> Integer
fromAccuracy Accuracy
ac))
  -- limit s = limit (s . getN)
  --   where
  --   -- q > 0 => 0 < 0.5^(getN q) <= q
  --   getN :: Rational -> Integer
  --   getN r = integer $ integerLog2 (max 1 $ (ceiling $ 1/r) - 1) + 1

instance HasLimits Int (CN MPBall -> CN MPBall) where
  type LimitType Int (CN MPBall -> CN MPBall) = (CN MPBall -> CN MPBall)
  limit :: (Int -> CN MPBall -> CN MPBall)
-> LimitType Int (CN MPBall -> CN MPBall)
limit Int -> CN MPBall -> CN MPBall
s = (Integer -> CN MPBall -> CN MPBall)
-> LimitType Integer (CN MPBall -> CN MPBall)
forall ix s. HasLimits ix s => (ix -> s) -> LimitType ix s
limit (Int -> CN MPBall -> CN MPBall
s (Int -> CN MPBall -> CN MPBall)
-> (Integer -> Int) -> Integer -> CN MPBall -> CN MPBall
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Int
c)
    where
    c :: Integer -> Int
    c :: Integer -> Int
c = Integer -> Int
forall t. CanBeInt t => t -> Int
int

  
limitMPBall :: (Accuracy -> ix) -> (ix -> (CN MPBall -> CN MPBall)) -> (CN MPBall -> CN MPBall)
limitMPBall :: (Accuracy -> ix)
-> (ix -> CN MPBall -> CN MPBall) -> CN MPBall -> CN MPBall
limitMPBall Accuracy -> ix
ac2ix ix -> CN MPBall -> CN MPBall
fs CN MPBall
x =
  maybeTrace ("limit (CN MPBall -> CN MPBall): x = " ++ show x) $
  maybeTrace ("limit (CN MPBall -> CN MPBall): xPNext = " ++ show xPNext) $
  maybeTrace ("limit (CN MPBall -> CN MPBall): accuracies = " ++ show accuracies) $
  [Accuracy] -> CN MPBall
tryAccuracies [Accuracy]
accuracies
  where
  acX :: Accuracy
acX = CN MPBall -> Accuracy
forall a. HasAccuracy a => a -> Accuracy
getFiniteAccuracy CN MPBall
x
  accuracies :: [Accuracy]
accuracies = DivIType Integer Integer -> [Accuracy]
forall t1.
(HasOrderAsymmetric (DivIType t1 Integer) Integer,
 ConvertibleExactly (DivIType t1 Integer) Accuracy,
 CanDivIMod t1 Integer,
 CanMulAsymmetric Integer (DivIType t1 Integer),
 MulType Integer (DivIType t1 Integer) ~ t1,
 OrderCompareType (DivIType t1 Integer) Integer ~ Bool) =>
DivIType t1 Integer -> [Accuracy]
aux (Accuracy -> Integer
fromAccuracy Accuracy
acX)
    where
    aux :: DivIType t1 Integer -> [Accuracy]
aux DivIType t1 Integer
a
      | DivIType t1 Integer
a DivIType t1 Integer
-> Integer -> OrderCompareType (DivIType t1 Integer) Integer
forall a b.
HasOrderAsymmetric a b =>
a -> b -> OrderCompareType a b
>= Integer
4 = DivIType t1 Integer -> Accuracy
forall t. ConvertibleExactly t Accuracy => t -> Accuracy
bits DivIType t1 Integer
a Accuracy -> [Accuracy] -> [Accuracy]
forall a. a -> [a] -> [a]
: DivIType t1 Integer -> [Accuracy]
aux ((Integer
100 Integer
-> DivIType t1 Integer -> MulType Integer (DivIType t1 Integer)
forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
* DivIType t1 Integer
a) t1 -> Integer -> DivIType t1 Integer
forall t1 t2. CanDivIMod t1 t2 => t1 -> t2 -> DivIType t1 t2
`divI` Integer
105)
      | Bool
otherwise = [DivIType t1 Integer -> Accuracy
forall t. ConvertibleExactly t Accuracy => t -> Accuracy
bits DivIType t1 Integer
a]
  xPNext :: CN MPBall
xPNext = Precision -> CN MPBall -> CN MPBall
forall t. CanSetPrecision t => Precision -> t -> t
setPrecision (Precision -> Precision
forall t. ConvertibleExactly t Integer => t -> Precision
increaseP (Precision -> Precision) -> Precision -> Precision
forall a b. (a -> b) -> a -> b
$ CN MPBall -> Precision
forall t. HasPrecision t => t -> Precision
getPrecision CN MPBall
x) CN MPBall
x
  increaseP :: t -> Precision
increaseP t
p =
    Integer -> Precision
prec (Integer -> Precision) -> Integer -> Precision
forall a b. (a -> b) -> a -> b
$ (t -> Integer
forall t. CanBeInteger t => t -> Integer
integer t
p) Integer -> Integer -> AddType Integer Integer
forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Integer
10
    -- prec $ ((101 * (integer p)) `div` 100) + 1

  tryAccuracies :: [Accuracy] -> CN MPBall
tryAccuracies [] =
    NumError -> CN MPBall
forall v. NumError -> CN v
CN.noValueNumErrorPotential (NumError -> CN MPBall) -> NumError -> CN MPBall
forall a b. (a -> b) -> a -> b
$ String -> NumError
CN.NumError String
"limit (CN MPBall -> CN MPBall) failed"
  tryAccuracies (Accuracy
ac : [Accuracy]
rest) =
    let result :: CN MPBall
result = Accuracy -> CN MPBall
withAccuracy Accuracy
ac in
    case CN MPBall -> Bool
forall es. CanTestErrorsPresent es => es -> Bool
CN.hasError CN MPBall
result of
      Bool
False -> CN MPBall
result
      Bool
_ -> [Accuracy] -> CN MPBall
tryAccuracies [Accuracy]
rest

  withAccuracy :: Accuracy -> CN MPBall
withAccuracy Accuracy
ac =
    maybeTrace ("limit (CN MPBall -> CN MPBall): withAccuracy: ac = " ++ show ac) $
    (ix -> CN MPBall -> CN MPBall
fs ix
ix CN MPBall
xPNext) CN MPBall -> Rational -> PlusMinusType (CN MPBall) Rational
forall t1 t2. CanPlusMinus t1 t2 => t1 -> t2 -> PlusMinusType t1 t2
+- Rational
PowType Rational Integer
e
    where
    ix :: ix
ix = Accuracy -> ix
ac2ix Accuracy
ac
    e :: PowType Rational Integer
e = Rational
0.5Rational -> Integer -> PowType Rational Integer
forall t1 t2. CanPow t1 t2 => t1 -> t2 -> PowType t1 t2
^(Accuracy -> Integer
fromAccuracy Accuracy
ac)

-- data WithLipschitz f = WithLipschitz f f

-- instance HasLimits Rational (WithLipschitz (MPBall -> CN MPBall)) where
--   type LimitType Rational (WithLipschitz (MPBall -> CN MPBall)) = (MPBall -> CN MPBall)
--   limit ffs xPre =
--     maybeTrace ("limit (MPBall -> CN MPBall): x = " ++ show x) $
--     maybeTrace ("limit (MPBall -> CN MPBall): xC = " ++ show xC) $
--     maybeTrace ("limit (MPBall -> CN MPBall): xE = " ++ show xE) $
--     maybeTrace ("limit (MPBall -> CN MPBall): accuracies = " ++ show accuracies) $
--     tryAccuracies accuracies
--     where
--     acX = getFiniteAccuracy x
--     accuracies = aux (fromAccuracy acX)
--       where
--       aux a
--         | a >= 4 = bits a : aux ((100 * a) `divI` 105)
--         | otherwise = [bits a]
--     x = increasePrec xPre
--       where
--       increasePrec z = setPrecision (inc $ getPrecision xPre) z
--       inc p =
--           prec $ (integer p) + 10
--           -- prec $ ((101 * (integer p)) `div` 100) + 1
--     xC = centreAsBall x
--     xE = radius x

--     tryAccuracies [] =
--       noValueNumErrorPotentialECN (Nothing :: Maybe MPBall) $
--         NumError "limit (MPBall -> CN MPBall) failed"
--     tryAccuracies (ac : rest) =
--       let result = withAccuracy ac in
--       case getErrorsCN result of
--         [] -> result
--         _ -> tryAccuracies rest

--     withAccuracy ac =
--       maybeTrace ("limit (MPBall -> CN MPBall): withAccuracy: ac = " ++ show ac) $
--       lift1CE (updateRadius (+ (errorBound e))) fx
--       where
--       e = 0.5^!(fromAccuracy ac)
--       WithLipschitz f f' = ffs e
--       fxC_CN = f xC
--       f'x_CN = f' x
--       fx :: CN MPBall
--       fx =
--           case (ensureNoCN fxC_CN, ensureNoCN f'x_CN) of
--             ((Just fxC, []), (Just f'x, [])) ->
--               cn $ updateRadius (+ (xE * (errorBound f'x))) fxC
--             _ ->
--               f x -- fallback

-- instance HasLimits Rational (WithLipschitz (CReal -> CReal)) where
--   type LimitType Rational (WithLipschitz (CReal -> CReal)) = (CReal -> CReal)
--   limit ffs x = limit fs x
--     where
--     fs e = f
--       where
--       WithLipschitz f _ = ffs e