{-# LANGUAGE PartialTypeSignatures #-}
{-# OPTIONS_GHC -Wno-missing-signatures #-}
{-# OPTIONS_GHC -Wno-partial-type-signatures #-}
{-|
    Module      :  AERN2.Real.Examples.Introduction
    Description :  aern2-real introductory examples
    Copyright   :  (c) Michal Konecny
    License     :  BSD3

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

    Introductory examples for packages aern2-mp and aern2-real.

    Please see aern2-real/README.md for explanations.

    You can run the following examples in ghci.
    If you installed AERN2 using the official instructions,
    you can start ghci using the following command in the base
    folder:

    @
    stack repl aern2-real/examples/AERN2/Real/Examples/Introduction.hs
    @
-}
module AERN2.Real.Examples.Introduction where

import MixedTypesNumPrelude

import qualified Numeric.CollectErrors as CN

import AERN2.MP
import AERN2.MP.WithCurrentPrec
import AERN2.Real

-- import Debug.Trace

------------------------------
-- real numbers
------------------------------

-- Start with a simple real number:

sine1 :: SinCosType Integer
sine1 = forall t. CanSinCos t => t -> SinCosType t
sin Integer
1

sine1_run1 :: ExtractedApproximation CReal Precision
sine1_run1 = CReal
sine1 forall e q.
CanExtractApproximation e q =>
e -> q -> ExtractedApproximation e q
? (Integer -> Precision
prec Integer
120)
-- result: [0.84147098480789650665250232... ± ~4.6644e-35 ~2^(-114)]
sine1_run2 :: ExtractedApproximation CReal Accuracy
sine1_run2 = CReal
sine1 forall e q.
CanExtractApproximation e q =>
e -> q -> ExtractedApproximation e q
? (forall t. ConvertibleExactly t Accuracy => t -> Accuracy
bits Integer
120)
-- result: [0.84147098480789650665250232... ± ~2.2431e-55 ~2^(-181)]

-- Next, do a bit more work:

sumSines1 :: Integer -> CReal
sumSines1 :: Integer -> CReal
sumSines1 Integer
n = forall t.
(CanAddSameType t, ConvertibleExactly Integer t) =>
[t] -> t
sum [forall t. CanSinCos t => t -> SinCosType t
sin Integer
i | Integer
i <- [Integer
1..Integer
n]]

-- Request the above expression with n = 100 using roughly 100 significant binary digits:
sumSines1_run1 :: CN MPBall
sumSines1_run1 :: CN MPBall
sumSines1_run1 = (Integer -> CReal
sumSines1 Integer
100) forall e q.
CanExtractApproximation e q =>
e -> q -> ExtractedApproximation e q
? (Integer -> Precision
prec Integer
120)
{- ghci log:

*AERN2.Real.Introduction> sumSines1_run1
[-0.12717101366042011543675217... ± ~2.8393e-33 ~2^(-108)]
(0.03 secs, 26,203,776 bytes)
-}

-- Same as above but request guaranteed 100 bits of accuracy:
sumSines1_run2 :: ExtractedApproximation CReal Accuracy
sumSines1_run2 = (Integer -> CReal
sumSines1 Integer
100) forall e q.
CanExtractApproximation e q =>
e -> q -> ExtractedApproximation e q
? (forall t. ConvertibleExactly t Accuracy => t -> Accuracy
bits Integer
100)
{- ghci log:

*AERN2.Real.Introduction> sumSines1_run2
[-0.12717101366042011543675217... ± ~2.8393e-33 ~2^(-108)]
(0.19 secs, 319,789,600 bytes)

This is considetably slower because there is some backtracking when target accuracy is not reached.  
-}

------------------------------
-- real number comparisons
------------------------------

{-
  First consider comparisons of real number approximations.
  These may be decided or undecided, using a 'Kleenean'.
-}

pi100 :: CN MPBall
pi100 :: CN MPBall
pi100 = CReal
piforall e q.
CanExtractApproximation e q =>
e -> q -> ExtractedApproximation e q
?(forall t. ConvertibleExactly t Accuracy => t -> Accuracy
bits Integer
100)

compare_run1 :: CN Kleenean
compare_run1 :: CN Kleenean
compare_run1 = CN MPBall
pi100 forall a b.
HasOrderAsymmetric a b =>
a -> b -> OrderCompareType a b
> Integer
0
-- returns: CertainTrue

compare_run2 :: CN Kleenean
compare_run2 :: CN Kleenean
compare_run2 = CN MPBall
pi100 forall a b. HasEqAsymmetric a b => a -> b -> EqCompareType a b
== CN MPBall
pi100
-- returns: TrueOrFalse

compare_run3 :: CKleenean
compare_run3 :: CSequence Kleenean
compare_run3 = CReal
pi forall a b.
HasOrderAsymmetric a b =>
a -> b -> OrderCompareType a b
> Integer
0
-- in ghci prints: {?(prec 36): CertainTrue}
-- (evaluated using default precision 36)

compare_run4 :: EqCompareType CReal (AddType CReal (PowType Integer Integer))
compare_run4 = CReal
pi forall a b. HasEqAsymmetric a b => a -> b -> EqCompareType a b
== CReal
pi forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Integer
2forall t1 t2. CanPow t1 t2 => t1 -> t2 -> PowType t1 t2
^(-Integer
100)
-- in ghci prints: {?(prec 36): TrueOrFalse}

compare_run5 :: ExtractedApproximation
  (EqCompareType CReal (AddType CReal (PowType Integer Integer)))
  Precision
compare_run5 = (CReal
pi forall a b. HasEqAsymmetric a b => a -> b -> EqCompareType a b
== CReal
pi forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ Integer
2forall t1 t2. CanPow t1 t2 => t1 -> t2 -> PowType t1 t2
^(-Integer
100)) forall e q.
CanExtractApproximation e q =>
e -> q -> ExtractedApproximation e q
? (Integer -> Precision
prec Integer
1000)
-- returns: CertainFalse

compare_run6 :: EqCompareType CReal Integer
compare_run6 = (forall t. CanBeCReal t => t -> CReal
creal Integer
0) forall a b. HasEqAsymmetric a b => a -> b -> EqCompareType a b
== Integer
0
-- in ghci prints: {?(prec 36): CertainTrue}
-- this is decided in finite time because 0 is represented exactly

compare_run7 :: ExtractedApproximation (CSequence Kleenean) Precision
compare_run7 = CReal
pi forall a b. HasEqAsymmetric a b => a -> b -> EqCompareType a b
== CReal
pi forall e q.
CanExtractApproximation e q =>
e -> q -> ExtractedApproximation e q
? (Integer -> Precision
prec Integer
10000)
-- returns: TrueOrFalse
-- (cannot confirm pi=pi in finite time)

------------------------------
-- checking partial functions
------------------------------

partialfn_bad1 :: SqrtType Integer
partialfn_bad1 = forall t. CanSqrt t => t -> SqrtType t
sqrt (-Integer
1)
{- ghci log:

*AERN2.Real.Introduction> partialfn_bad1 ? (bits 100)
{{ERROR: out of domain: negative sqrt argument}}

-}

a_third :: CReal
a_third = forall t. CanBeCReal t => t -> CReal
creal (Integer
1forall t1 t2. CanDiv t1 t2 => t1 -> t2 -> DivType t1 t2
/Integer
3)

partialfn_bad2 :: DivType Integer CReal
partialfn_bad2 = Integer
1forall t1 t2. CanDiv t1 t2 => t1 -> t2 -> DivType t1 t2
/(CReal
a_thirdforall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
-CReal
a_third)
{- ghci log:

*AERN2.Real.Introduction> partialfn_bad2 ? (prec 100)
{{POTENTIAL ERROR: division by 0}}

*AERN2.Real.Introduction> partialfn_bad2 ? (bits 100)
{{POTENTIAL ERROR: numeric error: failed to find an approximation with sufficient accuracy}}

-}

partialfn_bad3 :: DivType Integer CReal
partialfn_bad3 = Integer
1forall t1 t2. CanDiv t1 t2 => t1 -> t2 -> DivType t1 t2
/(CReal
piforall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
-CReal
pi)
{- ghci log:

*AERN2.Real.Introduction> partialfn_bad3 ? (prec 100)
{{POTENTIAL ERROR: division by 0}}

*AERN2.Real.Introduction> partialfn_bad3 ? (bits 100)
-- TAKES A VERY LONG TIME

-}

{-
 When computing on approximations which do not have enough information
 to check whether an error occurs, we get a *potential* error:
-}

partialfn_ok4 :: SqrtType CReal
partialfn_ok4 = forall t. CanSqrt t => t -> SqrtType t
sqrt (CReal
piforall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
-CReal
pi)
{- ghci log:

*AERN2.Real.Introduction> partialfn_ok4 ? (prec 100)
[0.00000000000000000061331736... ± ~6.1332e-19 ~2^(-60)]{{POTENTIAL ERROR: out of domain: negative sqrt argument}}
-}

partialfn_ok5 :: SqrtType CReal
partialfn_ok5 = forall cnt. CanClearPotentialErrors cnt => cnt -> cnt
clearPotentialErrors (forall t. CanSqrt t => t -> SqrtType t
sqrt (CReal
piforall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
-CReal
pi))
{- ghci log:

*AERN2.Real.Introduction> partialfn_ok5 ? (prec 100)
[0.00000000000000000061331736... ± ~6.1332e-19 ~2^(-60)]

-}

partialfn_bad6 :: SqrtType (SubType CReal Integer)
partialfn_bad6 = forall cnt. CanClearPotentialErrors cnt => cnt -> cnt
clearPotentialErrors (forall t. CanSqrt t => t -> SqrtType t
sqrt (CReal
piforall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
-CReal
piforall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
-Integer
1))
{- ghci log:

*AERN2.Real.Introduction> partialfn_bad6 ? (prec 100) 
{{ERROR: out of domain: negative sqrt argument}}

-}

partialfn_bad7 :: SqrtType (SubType CReal (PowType Integer Integer))
partialfn_bad7 = forall cnt. CanClearPotentialErrors cnt => cnt -> cnt
clearPotentialErrors (forall t. CanSqrt t => t -> SqrtType t
sqrt (CReal
piforall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
-CReal
piforall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
-Integer
2forall t1 t2. CanPow t1 t2 => t1 -> t2 -> PowType t1 t2
^(-Integer
1000)))
{- ghci log:

*AERN2.Real.Introduction> partialfn_bad7 ? (prec 100)
[0.00000000000000000061331736... ± ~6.1332e-19 ~2^(-60)]

*AERN2.Real.Introduction> partialfn_bad7 ? (prec 1000)
{{ERROR: out of domain: negative sqrt argument}}

-}

detectCN :: CN.CanTestErrorsPresent a => a -> Maybe a
detectCN :: forall a. CanTestErrorsPresent a => a -> Maybe a
detectCN a
r = if forall t. CanNeg t => t -> NegType t
not (forall es. CanTestErrorsPresent es => es -> Bool
CN.hasError a
r) then forall a. a -> Maybe a
Just a
r else forall a. Maybe a
Nothing
{- ghci log:

*AERN2.Real.Introduction> detectCN (sqrt (-1) ? (prec 100))
Nothing

*AERN2.Real.Introduction> detectCN (sqrt 0 ? (prec 100))
Just [0 ± 0]

-}

---------------------------------
-- Computing limits
--------------------------------

fact :: Integer -> CReal
fact :: Integer -> CReal
fact Integer
n = forall t. CanBeCReal t => t -> CReal
creal forall a b. (a -> b) -> a -> b
$ forall t.
(CanMulSameType t, ConvertibleExactly Integer t) =>
[t] -> t
product [Integer
1..Integer
n]

e_sum :: Integer -> CReal
e_sum :: Integer -> CReal
e_sum Integer
n = forall t.
(CanAddSameType t, ConvertibleExactly Integer t) =>
[t] -> t
sum forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (forall t. CanRecip t => t -> DivType Integer t
recip forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> CReal
fact) [Integer
0..Integer
n]

my_e :: CReal
my_e :: CReal
my_e = forall ix s. HasLimits ix s => (ix -> s) -> LimitType ix s
limit forall a b. (a -> b) -> a -> b
$ \(Integer
n :: Integer) -> Integer -> CReal
e_sum (Integer
nforall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+Integer
2)

{- ghci log:

*AERN2.Real.Introduction> my_e ? (prec 1000)
[2.71828182845904523536028747... ± ~0.0000 ~2^(-1217)]

-}

-- a faster version:

e_sum2 :: Integer -> CReal
e_sum2 :: Integer -> CReal
e_sum2 Integer
n = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl forall {t1} {t2}.
(CanAddAsymmetric Integer (DivType t1 t2), CanDiv t1 t2) =>
t1 -> t2 -> AddType Integer (DivType t1 t2)
aux (forall t. CanBeCReal t => t -> CReal
creal Integer
1) forall a b. (a -> b) -> a -> b
$ forall a. [a] -> [a]
reverse [Integer
1..Integer
n]
  where aux :: t1 -> t2 -> AddType Integer (DivType t1 t2)
aux t1
x t2
m = Integer
1 forall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+ t1
x forall t1 t2. CanDiv t1 t2 => t1 -> t2 -> DivType t1 t2
/ t2
m

my_e2 :: CReal
my_e2 :: CReal
my_e2 = forall ix s. HasLimits ix s => (ix -> s) -> LimitType ix s
limit forall a b. (a -> b) -> a -> b
$ \(Integer
n :: Integer) -> Integer -> CReal
e_sum2 (Integer
nforall t1 t2. CanAddAsymmetric t1 t2 => t1 -> t2 -> AddType t1 t2
+Integer
2)

---------------------------------
-- "parallel" branching for real numbers
--------------------------------

absQ :: Rational -> Rational
absQ :: Rational -> Rational
absQ Rational
x = if Rational
x forall a b.
HasOrderAsymmetric a b =>
a -> b -> OrderCompareType a b
< Integer
0 then -Rational
x else Rational
x

absR1 :: CReal -> CReal
absR1 :: CReal -> CReal
absR1 CReal
x = if CReal
x forall a b.
HasOrderAsymmetric a b =>
a -> b -> OrderCompareType a b
< Integer
0 then -CReal
x else CReal
x


pif_run1 :: CReal
pif_run1 = CReal -> CReal
absR1 (CReal
piforall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
-CReal
pi)
{- ghci log:

*AERN2.Real.Introduction> pif_run1
{?(prec 36): [0 ± ~2.9104e-11 ~2^(-35)]}

-}

-- pif_run2 = foldl1 (.) (replicate 100 (absR1 . (100*))) (pi-pi)

absR2_approx :: t
-> Rational
-> IfThenElseType (SelectType (OrderCompareType t Rational)) t
absR2_approx t
x (Rational
q :: Rational) = if forall k. CanSelect k => k -> k -> SelectType k
select (t
x forall a b.
HasOrderAsymmetric a b =>
a -> b -> OrderCompareType a b
> -Rational
q) (t
x forall a b.
HasOrderAsymmetric a b =>
a -> b -> OrderCompareType a b
< Rational
q) then t
x else -t
x

absR2 :: CReal -> CReal
absR2 :: CReal -> CReal
absR2 CReal
x = forall ix s. HasLimits ix s => (ix -> s) -> LimitType ix s
limit forall a b. (a -> b) -> a -> b
$ forall {t}.
(NegType t ~ t,
 HasIfThenElse (SelectType (OrderCompareType t Rational)) t,
 CanSelect (OrderCompareType t Rational),
 HasOrderAsymmetric t Rational, CanNeg t) =>
t
-> Rational
-> IfThenElseType (SelectType (OrderCompareType t Rational)) t
absR2_approx CReal
x

select_run1 :: CReal
select_run1 = CReal -> CReal
absR2 (CReal
piforall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
-CReal
pi)
{- ghci log:

*AERN2.Real.Introduction> select_run1
{?(prec 36): [0 ± ~4.3656e-11 ~2^(-34)]}

-}


logistic :: (RealNumber t) => Rational -> Integer -> t -> t
logistic :: forall t. RealNumber t => Rational -> Integer -> t -> t
logistic Rational
c Integer
n t
x0 =
  (forall (t :: * -> *) a. Foldable t => (a -> a -> a) -> t a -> a
foldl1 forall b c a. (b -> c) -> (a -> b) -> a -> c
(.) (forall n a. CanBeInteger n => n -> a -> [a]
replicate Integer
n t -> MulType t t
lg)) t
x0
  where
  lg :: t -> MulType t t
lg t
x = Rational
c forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
* t
x forall t1 t2. CanMulAsymmetric t1 t2 => t1 -> t2 -> MulType t1 t2
* (Integer
1forall t1 t2. CanSub t1 t2 => t1 -> t2 -> SubType t1 t2
-t
x)

logistic1 :: (RealNumber t) => Integer -> t
logistic1 :: forall t. RealNumber t => Integer -> t
logistic1 Integer
n = forall t. RealNumber t => Rational -> Integer -> t -> t
logistic Rational
3.82 Integer
n (forall t1 t2. ConvertibleExactly t1 t2 => t1 -> t2
convertExactly Rational
0.5)

logistic1_CReal_run :: Integer -> CReal
logistic1_CReal_run :: Integer -> CReal
logistic1_CReal_run Integer
n = 
  forall t. RealNumber t => Integer -> t
logistic1 Integer
n

logistic1_WithCurrentPrec_run :: Integer -> CReal
logistic1_WithCurrentPrec_run :: Integer -> CReal
logistic1_WithCurrentPrec_run Integer
n = 
  (forall (p :: Nat). KnownNat p => WithCurrentPrec p (CN MPBall))
-> CReal
crealFromWithCurrentPrec forall a b. (a -> b) -> a -> b
$ forall t. RealNumber t => Integer -> t
logistic1 Integer
n

logistic1_WithCurrentPrec_p_run :: Integer -> Precision -> CN MPBall
logistic1_WithCurrentPrec_p_run :: Integer -> Precision -> CN MPBall
logistic1_WithCurrentPrec_p_run Integer
n Precision
p = 
  forall t.
Precision
-> (forall (p :: Nat). KnownNat p => WithCurrentPrec p t) -> t
runWithPrec Precision
p forall a b. (a -> b) -> a -> b
$ forall t. RealNumber t => Integer -> t
logistic1 Integer
n

{-  Example uses:

*AERN2.Real.Examples.Introduction> logistic1_CReal_run 100 ? (bits 100)
[0.95087585116480286419338875... ± ~2.9792e-32 ~2^(-104)]
  
*AERN2.Real.Examples.Introduction> logistic1_CReal_run 10000 ? (bits 100)
[0.20775682944252359241450861... ± ~0.0000 ~2^(-2566)]
(2.06 secs, 2,970,188,704 bytes)

*AERN2.Real.Examples.Introduction> logistic1_WithCurrentPrec_run 10000 ? (bits 100)
[0.20775682944252359241450861... ± ~0.0000 ~2^(-2566)]
(2.02 secs, 2,897,034,816 bytes)

*AERN2.Real.Examples.Introduction> logistic1_WithCurrentPrec_p_run 10000 (prec 20000)
[0.20775682944252359241450861... ± ~1.0317e-200 ~2^(-664)]
(1.05 secs, 1,421,858,848 bytes)

-}

{-
  Recommended further reading:  ClosestPairDist.hs
-}