{-|
    Module      :  Data.Number.ER.Real.Approx.Sequence
    Description :  exact reals via convergent sequences
    Copyright   :  (c) Michal Konecny
    License     :  BSD3

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

    Types and methods related to explicit 
    convergent sequences of real number approximations.
-}
module Data.Number.ER.Real.Approx.Sequence 
(
    ConvergRealSeq(..),
    makeFastConvergRealSeq,
    convertFuncRA2Seq,
    convertBinFuncRA2Seq,
    convergRealSeqElem,
    showConvergRealSeq,
    showConvergRealSeqAuto
)
where

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

import Data.Maybe
import Data.Ratio

{-|
  A converging sequence of real number approximations.
  
  * Every finite subsequence has a non-empty intersection.
  
  * The limit should be a singleton.
-}
data ConvergRealSeq ra =
    ConvergRealSeq (EffortIndex -> ra)

convergRealSeqElem :: (ConvergRealSeq ra) -> EffortIndex -> ra
convergRealSeqElem (ConvergRealSeq sq) ix = sq ix
        
{-| 
    Using this operator, a unary funtion working over
    approximations can be converted to one that works
    over exact numbers represented through a sequence
    of approximations.
-}
convertFuncRA2Seq ::
    (EffortIndex -> ra -> ra) ->
    (ConvergRealSeq ra) ->
    (ConvergRealSeq ra)
convertFuncRA2Seq f (ConvergRealSeq argSeq) = 
    ConvergRealSeq resultSeq
    where
    resultSeq ix =
        f ix (argSeq ix)
        
{-|
    The same as above, where f is binary
-}            
convertBinFuncRA2Seq :: 
    (EffortIndex -> ra -> ra -> ra) -> 
    (ConvergRealSeq ra) -> 
    (ConvergRealSeq ra) -> 
    (ConvergRealSeq ra)
    
convertBinFuncRA2Seq f (ConvergRealSeq arg1) (ConvergRealSeq arg2) = 
    ConvergRealSeq resultSeq
    where
    resultSeq ix =
        f ix (arg1 ix) (arg2 ix)

{-|
    Turn an arbitrary convergent sequence into one with
    a guaranteed convergence rate - the precision (as defined
    by 'RA.ERApprox.RA.getPrecision') of x_ix is at least ix.
-}
makeFastConvergRealSeq :: 
    (RA.ERApprox ra) => 
    (ConvergRealSeq ra) -> 
    (ConvergRealSeq ra)
makeFastConvergRealSeq (ConvergRealSeq argSeq) = 
    ConvergRealSeq fastSeq
    where
    fastSeq ix =
        head $ catMaybes $ map (precisionOK . argSeq) indexSeries
        where
        indexSeries =
    --        take 5 $ -- upper bound on iteration - for testing
            binGeomSeries (max 1 ix)
        precisionOK ra
            | RA.getPrecision ra >= (effIx2prec ix) = Just ra
            | otherwise = Nothing

{-| 
    binGeomSeries n is the geometric series
    [ n, 2n, 4n, 8n, ...]
-}    
binGeomSeries
    :: (Num a)
    => a
    -> [a]
binGeomSeries n =
    n : (binGeomSeries (2 * n))

instance (RA.ERApprox ra) => Show (ConvergRealSeq ra) 
    where
    show = showConvergRealSeq 6 True False 10 -- cheating here, should throw an error


{-|
    Show function for ConvergRealSeq's with full arguments.
-}    
showConvergRealSeq
    :: (RA.ERApprox ra)
    => Int
    -> Bool
    -> Bool
    -> Precision
    -> (ConvergRealSeq ra)
    -> String

showConvergRealSeq numDigits showGran showComponents prec r =
    RA.showApprox numDigits showGran showComponents $
         convergRealSeqElem (makeFastConvergRealSeq r) (prec2effIx prec)


{-|
    Show function for ConvergRealSeq's with all parameters fixed
    except for number of digits
-}
showConvergRealSeqAuto 
    :: (RA.ERApprox ra)
    => Int
    -> (ConvergRealSeq ra)
    -> String
showConvergRealSeqAuto numDigits argSeq =
    showConvergRealSeq numDigits True False prec argSeq
    where
    prec = effIx2prec $ ceiling $ (fromInteger $ toInteger numDigits) * 3.3219280948873626



instance
    (RA.ERApprox ra)
    => Eq (ConvergRealSeq ra)
    where
    r1 == r2 = 
        iterateRA_A raEq 2 [r1, r2]
        where
        raEq _ ([a1,a2]) = RA.equalReals a1 a2
                
instance
    (RA.ERApprox ra)
    => Ord (ConvergRealSeq ra)
    where
    compare r1 r2 = 
        iterateRA_A eraComp 2 [r1, r2]
        where
        eraComp _ ([a1,a2]) = RA.compareReals a1 a2
            
pointwiseConvergRealSeq1 f (ConvergRealSeq sq) =
    ConvergRealSeq (f . sq)
pointwiseConvergRealSeq2 f (ConvergRealSeq sq1) (ConvergRealSeq sq2) =
    ConvergRealSeq (\ix -> f (sq1 ix) (sq2 ix))
            
instance 
    (RA.ERApprox ra)
    => Num (ConvergRealSeq ra)
    where
    fromInteger n = ConvergRealSeq sq
        where
        sq ix =
            RA.setMinGranularityOuter (effIx2gran ix) $ fromInteger n
    abs = pointwiseConvergRealSeq1 $ abs
    signum = pointwiseConvergRealSeq1 $ signum
    negate = pointwiseConvergRealSeq1 $ negate
    (+) = pointwiseConvergRealSeq2 $ (+)
    (*) = pointwiseConvergRealSeq2 $ (*)
    
instance
    (RA.ERApprox ra)
    => Fractional (ConvergRealSeq ra)
    where
    fromRational q = ConvergRealSeq sq
        where
        sq ix =
            (RA.setMinGranularityOuter (effIx2gran ix) num) / denom
        num = fromInteger $ numerator q
        denom = fromInteger $ denominator q
    recip = pointwiseConvergRealSeq1 $ recip

{-|
    Take a converging sequence of partial functions F_i that operate on 
    real approximations and turn it into a function F that operates on converging sequences. 
    F looks for some members of the real approximation sequences 
    and an i so that F_i is defined for the chosen approximations
    and returns its result.  
-}
iterateRA_A
    :: (EffortIndex -> [ra] -> Maybe a) 
        -- ^ a sequence of partial functions based on approximations
    -> EffortIndex -- ^ a starting index to use when searching sequences
    -> ([ConvergRealSeq ra] -> a) 
        -- ^ a total function based on sequences

iterateRA_A fn_RA startIx args =
    head $ catMaybes $ map ((uncurry fn_RA) . args_Prec) indexSeries
    where
    indexSeries =
--        take 5 $ -- upper bound on iteration - for testing
        binGeomSeries $ max 1 startIx
        -- [(max 1 startIx)..]
    args_Prec currentIndex =
        (currentIndex, map (\ arg -> convergRealSeqElem arg currentIndex) args)