Limit.hs

{-# 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 = limitMPBall (\ac -> 0.5^(fromAccuracy ac))

instance HasLimits Integer (CN MPBall -> CN MPBall) where
  type LimitType Integer (CN MPBall -> CN MPBall) = (CN MPBall -> CN MPBall)
  limit = limitMPBall (\ac -> (fromAccuracy 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 s = limit (s . c)
    where
    c :: Integer -> Int
    c = int

  
limitMPBall :: (Accuracy -> ix) -> (ix -> (CN MPBall -> CN MPBall)) -> (CN MPBall -> CN MPBall)
limitMPBall ac2ix fs 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) $
  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]
  xPNext = setPrecision (increaseP $ getPrecision x) x
  increaseP p =
    prec $ (integer p) + 10
    -- prec $ ((101 * (integer p)) `div` 100) + 1

  tryAccuracies [] =
    CN.noValueNumErrorPotential $ CN.NumError "limit (CN MPBall -> CN MPBall) failed"
  tryAccuracies (ac : rest) =
    let result = withAccuracy ac in
    case CN.hasError result of
      False -> result
      _ -> tryAccuracies rest

  withAccuracy ac =
    maybeTrace ("limit (CN MPBall -> CN MPBall): withAccuracy: ac = " ++ show ac) $
    (fs ix xPNext) +- e
    where
    ix = ac2ix ac
    e = 0.5^(fromAccuracy 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