{-# LANGUAGE BangPatterns,GADTs,TypeSynonymInstances,FlexibleInstances #-}
{-|
= Computable Real Arithmetic
This module provides the data type 'CR' that implements the real closed field of computable real numbers.

== Centred Dyadic Approximations
The computable reals are realised as lists of rapidly shrinking intervals. The intervals used here are centred dyadic intervals, implemented here as the data type 'Approx'.

For more information on the theoretical aspects see <http://cs.swan.ac.uk/~csjens/pdf/centred.pdf>.
-}
module Data.CDAR.Approx (Approx(..)
                        ,CR(..)
--                        ,errorBits
--                        ,errorBound
--                        ,defaultPrecision
                        ,Precision
                        ,showA
                        ,showInBaseA
                        ,mBound
                        ,mapMB
                        ,setMB
                        ,enforceMB
                        ,approxAutoMB
                        ,approxMB
                        ,approxMB2
                        ,endToApprox
                        ,lowerBound
                        ,upperBound
                        ,lowerA
                        ,upperA
                        ,centre
                        ,centreA
                        ,radius
                        ,diameter
                        ,exact
                        ,approximatedBy
                        ,better
                        ,fromDyadic
                        ,fromDyadicMB
                        ,toApprox
                        ,toApproxMB
                        ,recipA
                        ,divAInteger
                        ,modA
                        ,divModA
                        ,toDoubleA
                        ,precision
                        ,significance
                        ,boundErrorTerm
                        ,boundErrorTermMB
                        ,limitSize
                        ,checkPrecisionLeft
                        ,limitAndBound
                        ,limitAndBoundMB
                        ,unionA
                        ,intersectionA
                        ,consistentA
                        ,poly
                        ,pow
                        ,powers
                        ,sqrtHeronA
                        ,sqrtA
                        ,sqrtRecA
                        ,findStartingValues
                        ,sqrtD
                        ,shiftD
                        ,sqrA
                        ,log2Factorials
                        ,taylorA
                        ,expA
                        ,expBinarySplittingA
                        ,expTaylorA
                        ,expTaylorA'
                        ,logA
                        ,logInternal
                        ,logBinarySplittingA
                        ,logTaylorA
                        ,sinTaylorA
                        ,sinTaylorRed1A
                        ,sinTaylorRed2A
                        ,sinA
                        ,cosA
                        ,atanA
                        ,sinBinarySplittingA
                        ,cosBinarySplittingA
                        ,atanBinarySplittingA
                        ,atanTaylorA
                        ,piRaw
                        ,piA
                        ,piMachinA
                        ,piBorweinA
                        ,piAgmA
                        ,log2A
                        ,lnSuperSizeUnknownPi
                        ,logAgmA
                        ,agmLn
                        ,showCRN
                        ,showCR
                        ,ok
                        ,require
                        ,toDouble
                        ,fromDouble
                        ,fromDoubleAsExactValue
                        ,polynomial
                        ,taylorCR
                        ,atanCR
                        ,piCRMachin
                        ,piMachinCR
                        ,piBorweinCR
                        ,piBinSplitCR
                        ,ln2
                        ,sinCRTaylor
                        ,sinCR
                        ,cosCR
                        ,sqrtCR
                        ,expCR
                        ) where

import           Control.Applicative (ZipList (..))
import           Control.DeepSeq
import           Control.Exception
import           Data.Bits
import           Data.Char (isDigit)
import           Data.CDAR.Classes
import           Data.CDAR.Dyadic
import           Data.CDAR.Extended
import           Data.CDAR.IntegerLog
import           Data.Char (intToDigit)
import           Data.List (findIndex, intersperse, transpose, unfoldr, zipWith4)
import           Data.Ratio

-- import Debug.Trace

-- |A type synonym. Used to denote number of bits after binary point.
type Precision = Int

{-|
= Centred Dyadic Approximations
There are two constructors for approximations:

- 'Approx' is encodes some finite interval with dyadic endpoints. A real
  number is /approximated/ by the approximation is it belongs to the interval.
- 'Bottom' is the trivial approximation that approximates all real numbers.

The four fields of an @Approx m e s@ should be thought of as:

[@mb@] the midpoint bound, ie maximum bits available for the midpoint
[@m@] the midpoint
[@e@] the error term
[@s@] the exponent

Thus, a value @Approx p m e s@ is to be interpreted as the interval
[(m-e)*2^s, (m+e)*2^s] where |m| <= 2^p.

== Centred intervals
We have opted for a centred representation of the intervals. It is also
possible to represent the endpoints as 'Dyadic' numbers. The rationale for a
centred repersentation is that we often normalise an approximation @Approx p m e
s@ so that @e@ is limited in size. This allows many operations to only work on
one large number @m@.

== Potential for overflow
Since the last field (the exponent) is only an 'Int' it may overflow. This is
an optimisation that was adopted since it seems unlikely that overflow in a 64
bit Int exponent would occur. In a 32 bit system, this is potentially an
issue.

The 'Integer' data type is unbonded, but is, of course, bounded by the
available memory available in the computer. No attempt has been made to check
for exhausted memory.

== Approximations as a Domain

Ordered by reverse inclusion the dyadic intervals encoded by the 'Approx'
approximations (including 'Bottom') constitute the compact elements of a Scott
domain /D/. (It is a substructure of the (algebraic) interval domain.)
We will identify our approximations with the compact elements of /D/.

Increasing sequences in /D/ have suprema. A sequence /converges/ if the length
of the approximations tend to zero. The supremum of a converging sequence is a
singleton set containing a real number. Let ρ be the map taking a converging
sequence to the unique real number in the supremum. The computations on
(computable) real numbers is via this representation map ρ.

There is no check that the sequences we have are in fact increasing, but we
are assuming that all sequences are pairwise consistent. We can thus create an
increasing sequence by considering the sequence of finite suprema. For
correctness, we have to ensure that all operations done on consistent
sequences result in consistent sequences. If non-consistent sequences are
somehow input we can make no guarantees at all about the computed value.

Note, that we cannot ensure that converging sequences are mapped to converging
sequences because of properties of computable real arithmetic. In particular,
at any discuntinuity, it is impossible to compute a converging sequence.
-}
data Approx = Approx Int Integer Integer Int
            | Bottom
              deriving (ReadPrec [Approx]
ReadPrec Approx
Int -> ReadS Approx
ReadS [Approx]
(Int -> ReadS Approx)
-> ReadS [Approx]
-> ReadPrec Approx
-> ReadPrec [Approx]
-> Read Approx
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
readListPrec :: ReadPrec [Approx]
$creadListPrec :: ReadPrec [Approx]
readPrec :: ReadPrec Approx
$creadPrec :: ReadPrec Approx
readList :: ReadS [Approx]
$creadList :: ReadS [Approx]
readsPrec :: Int -> ReadS Approx
$creadsPrec :: Int -> ReadS Approx
Read,Int -> Approx -> ShowS
[Approx] -> ShowS
Approx -> String
(Int -> Approx -> ShowS)
-> (Approx -> String) -> ([Approx] -> ShowS) -> Show Approx
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Approx] -> ShowS
$cshowList :: [Approx] -> ShowS
show :: Approx -> String
$cshow :: Approx -> String
showsPrec :: Int -> Approx -> ShowS
$cshowsPrec :: Int -> Approx -> ShowS
Show)

instance NFData Approx where
    rnf :: Approx -> ()
rnf Approx
Bottom = ()
    rnf (Approx Int
b Integer
m Integer
e Int
s) = Int -> ()
forall a. NFData a => a -> ()
rnf Int
b () -> () -> ()
`seq` Integer -> ()
forall a. NFData a => a -> ()
rnf Integer
m () -> () -> ()
`seq` Integer -> ()
forall a. NFData a => a -> ()
rnf Integer
e () -> () -> ()
`seq` Int -> ()
forall a. NFData a => a -> ()
rnf Int
s

instance Scalable Approx where
  scale :: Approx -> Int -> Approx
scale Approx
Bottom Int
_ = Approx
Bottom
  scale (Approx Int
b Integer
m Integer
e Int
s) Int
n = Int -> Integer -> Integer -> Int -> Approx
Approx Int
b Integer
m Integer
e (Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
n)

instance Scalable CR where
  scale :: CR -> Int -> CR
scale (CR ZipList Approx
x) Int
n = ZipList Approx -> CR
CR (ZipList Approx -> CR) -> ZipList Approx -> CR
forall a b. (a -> b) -> a -> b
$ (Approx -> Int -> Approx) -> Int -> Approx -> Approx
forall a b c. (a -> b -> c) -> b -> a -> c
flip Approx -> Int -> Approx
forall a. Scalable a => a -> Int -> a
scale Int
n (Approx -> Approx) -> ZipList Approx -> ZipList Approx
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> ZipList Approx
x

{-|
=The Computable Real data type

Computable reals are realised as infinite sequences of centred dyadic
representations.

All approximations in such a sequence should be pairwise consistent, i.e.,
have a non-empty intersection. However, there is no check that this is
actually the case.

If the diameter of the approximations tend to zero we say that the sequences
converges to the unique real number in the intersection of all intervals.
Given the domain /D/ of approximations described in 'Approx', we have a
representation (a retraction) ρ from the converging sequences in /D/ to ℝ.
Some operations on computable reals are partial, notably equality and
rnf m `seq`  there is no guarantee that a
rnf m `seq` 
rnf m `seq` 
rnf m `seq` there is a bound on how much effort is
rnf m `seq` mation. For involved computations it is
rnf m `seq` proximations are trivial, i.e.,
rnf m `seq` ually converge, it will generate proper
rnf m `seq` initial trivial approximations.
rnf m `seq` 
The amount of added effort in each iteration is rather substantial so the
expected precision of approximations increase very quickly.

==The actual data type

In fact, the type 'CR' is a newtype of 'ZipList' 'Approx' in the
implementation of infinite sequences of approximations, as that allows for
applicative style. Hopefully, it is not needed to access the internal
representation of 'CR' directly.
-}
newtype CR = CR {CR -> ZipList Approx
unCR :: ZipList Approx}

-- |Number of bits that error term is allowed to take up. A larger size allows
-- for more precise but slightly more costly computations. The value here is
-- suggested by test runs.
errorBits :: Int
errorBits :: Int
errorBits = Int
10

errorBound :: Integer
errorBound :: Integer
errorBound = Integer
2Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^Int
errorBits

errorBitsMB :: Int
errorBitsMB :: Int
errorBitsMB = Int
1

errorBoundMB :: Integer
errorBoundMB :: Integer
errorBoundMB = Integer
2Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^Int
errorBitsMB


-- |The default cutoff for diverging computations. May well be chosen much
-- smaller. 31 corresponds to about 10 decimal places.
defaultPrecision :: Precision
defaultPrecision :: Int
defaultPrecision = Int
31

{-|

Gives a decimal representation of an approximation. It tries to give as many
decimal digits as possible given the precision of the approximation. The
representation may be wrong by 1 ulp (unit in last place). If the value is not
exact the representation will be followed by @~@.

The representation is not always intuitive:

>>> showA (Approx 1 1 0)
"1.~"

The meaning of the above is that it is 1, but then the added @~@ (which must
be after the decimal point) means that the last position may be off by 1,
i.e., it could be down to 0 or up to 2. And [0,2] is indeed the range encoded
by the above approximation.
-}
showA :: Approx -> String
showA :: Approx -> String
showA = Int -> Approx -> String
showInBaseA Int
10

-- |Similar to 'showA' but can generate representations in other bases (<= 16).
{- am is the absolute value of the significand
   b corresponds to the value 1 with respect to the shift s -- this is used to find the digits in the auxiliary functions
   i is the integral part of am
   f is the fractional part of am
   i' and f' are the integral and fractional parts relevant for near zero approximations
   e' is the error term shifted appropriately when s positive, also set to at least 1
     (otherwise odd bases will yield infinite expansions)
-}
showInBaseA :: Int -> Approx -> String
showInBaseA :: Int -> Approx -> String
showInBaseA Int
_ Approx
Bottom = String
"⊥"
showInBaseA Int
base (Approx Int
_ Integer
m Integer
e Int
s)
    | Integer
e Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 Bool -> Bool -> Bool
&& (Int -> Bool
forall a. Integral a => a -> Bool
even Int
base Bool -> Bool -> Bool
|| Int
s Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0)
                     = String
sign String -> ShowS
forall a. [a] -> [a] -> [a]
++ Int -> Integer -> Integer -> Integer -> String
showExactA Int
base Integer
b Integer
i Integer
f
    | Integer
am Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
e         = String
"±" String -> ShowS
forall a. [a] -> [a] -> [a]
++ Int -> Integer -> Integer -> Integer -> String
showNearZeroA Int
base Integer
b Integer
i' Integer
f'
    | Bool
otherwise      = String
sign String -> ShowS
forall a. [a] -> [a] -> [a]
++ Int -> Integer -> Integer -> Integer -> Integer -> String
showInexactA Int
base Integer
b Integer
i Integer
f Integer
e'
    where b :: Integer
b = Int -> Integer
forall a. Bits a => Int -> a
bit (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
0 (-Int
s))
          am :: Integer
am = Integer -> Integer
forall a. Num a => a -> a
abs Integer
m
          i :: Integer
i = Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
shift Integer
am Int
s
          e' :: Integer
e' = Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
max Integer
1 (Integer -> Integer) -> Integer -> Integer
forall a b. (a -> b) -> a -> b
$ Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
shift Integer
e (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
0 Int
s)
          f :: Integer
f = Integer
am Integer -> Integer -> Integer
forall a. Bits a => a -> a -> a
.&. (Integer
bInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
1)
          i' :: Integer
i' = Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
shift (Integer
amInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
e) Int
s
          f' :: Integer
f' = (Integer
amInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
e) Integer -> Integer -> Integer
forall a. Bits a => a -> a -> a
.&. (Integer
bInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
1)
          sign :: String
sign = if Integer
m Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 then String
"-" else String
""

showExactA :: Int -> Integer -> Integer -> Integer -> String
showExactA :: Int -> Integer -> Integer -> Integer -> String
showExactA Int
base Integer
b Integer
i Integer
f =
    let g :: Integer -> Maybe (Char, Integer)
g Integer
i' = let (Integer
q,Integer
r) = Integer -> Integer -> (Integer, Integer)
forall a. Integral a => a -> a -> (a, a)
quotRem Integer
i' (Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
base)
               in if Integer
i' Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 then Maybe (Char, Integer)
forall a. Maybe a
Nothing
                  else (Char, Integer) -> Maybe (Char, Integer)
forall a. a -> Maybe a
Just (Int -> Char
intToDigit (Integer -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
r), Integer
q)
        ip :: String
ip = ShowS
forall a. [a] -> [a]
reverse ((Integer -> Maybe (Char, Integer)) -> Integer -> String
forall b a. (b -> Maybe (a, b)) -> b -> [a]
unfoldr Integer -> Maybe (Char, Integer)
g Integer
i)
        h :: Integer -> Maybe (Char, Integer)
h Integer
f' = let (Integer
q,Integer
r) = Integer -> Integer -> (Integer, Integer)
forall a. Integral a => a -> a -> (a, a)
quotRem ((Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
base)Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
f') Integer
b
               in if Integer
f' Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 then Maybe (Char, Integer)
forall a. Maybe a
Nothing
                  else (Char, Integer) -> Maybe (Char, Integer)
forall a. a -> Maybe a
Just (Int -> Char
intToDigit (Integer -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
q), Integer
r)
        fp :: String
fp = (Integer -> Maybe (Char, Integer)) -> Integer -> String
forall b a. (b -> Maybe (a, b)) -> b -> [a]
unfoldr Integer -> Maybe (Char, Integer)
h Integer
f
    in (if String -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null String
ip then String
"0" else String
ip)
       String -> ShowS
forall a. [a] -> [a] -> [a]
++ (if String -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null String
fp then String
"" else String
".")
       String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
fp

showNearZeroA :: Int -> Integer -> Integer -> Integer -> String
showNearZeroA :: Int -> Integer -> Integer -> Integer -> String
showNearZeroA Int
base Integer
b Integer
i Integer
f =
    let s :: String
s = Int -> Integer -> Integer -> Integer -> String
showExactA Int
base Integer
b Integer
i Integer
f
        t :: String
t = (Char -> Bool) -> ShowS
forall a. (a -> Bool) -> [a] -> [a]
takeWhile ((Char -> String -> Bool) -> String -> Char -> Bool
forall a b c. (a -> b -> c) -> b -> a -> c
flip Char -> String -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
elem String
"0.~") String
s
        u :: String
u = (Char -> Bool) -> ShowS
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Char -> Char -> Bool
forall a. Eq a => a -> a -> Bool
/= Char
'.') String
s
    in if String -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null String
t
       then Int -> Char -> String
forall a. Int -> a -> [a]
replicate (String -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length String
u) Char
'~'
       else String
t String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"~"

showInexactA :: Int -> Integer -> Integer -> Integer -> Integer -> String
showInexactA :: Int -> Integer -> Integer -> Integer -> Integer -> String
showInexactA Int
base Integer
b Integer
i Integer
f Integer
e =
    let g :: (Integer, Integer, Integer)
-> Maybe (Char, (Integer, Integer, Integer))
g (Integer
0,Integer
b',Integer
f') = if Integer
b' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
f'Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
e
                      then (Char, (Integer, Integer, Integer))
-> Maybe (Char, (Integer, Integer, Integer))
forall a. a -> Maybe a
Just (Char
'1', (Integer
0, (Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
base)Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
b', Integer
f'))
                      else Maybe (Char, (Integer, Integer, Integer))
forall a. Maybe a
Nothing
        g (Integer
n,Integer
b',Integer
f') = let (Integer
q,Integer
r) = Integer -> Integer -> (Integer, Integer)
forall a. Integral a => a -> a -> (a, a)
quotRem Integer
n (Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
base)
                          z :: (Integer, Integer, Integer)
z = (Integer
q, (Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
base)Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
b', Integer
rInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
b'Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
f')
                      in if Integer
eInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
f' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
b'
                         then (Char, (Integer, Integer, Integer))
-> Maybe (Char, (Integer, Integer, Integer))
forall a. a -> Maybe a
Just (Int -> Char
intToDigit (Integer -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
r), (Integer, Integer, Integer)
z)
                         else if Integer
e Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
min Integer
f' Integer
b'
                              then (Char, (Integer, Integer, Integer))
-> Maybe (Char, (Integer, Integer, Integer))
forall a. a -> Maybe a
Just (Int -> Char
intToDigit ((Integer -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
r Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`rem` (Int -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
base)), (Integer, Integer, Integer)
z)
                              else (Char, (Integer, Integer, Integer))
-> Maybe (Char, (Integer, Integer, Integer))
forall a. a -> Maybe a
Just (Char
'~', (Integer, Integer, Integer)
z)
        intRev :: String
intRev = ((Integer, Integer, Integer)
 -> Maybe (Char, (Integer, Integer, Integer)))
-> (Integer, Integer, Integer) -> String
forall b a. (b -> Maybe (a, b)) -> b -> [a]
unfoldr (Integer, Integer, Integer)
-> Maybe (Char, (Integer, Integer, Integer))
g (Integer
i,Integer
b,Integer
f)
        noFrac :: Bool
noFrac = case String
intRev of
                   [] -> Bool
False
                   (Char
x:String
_) -> Char
x Char -> Char -> Bool
forall a. Eq a => a -> a -> Bool
== Char
'~'
        int :: String
int = if String -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null String
intRev then String
"0" else ShowS
forall a. [a] -> [a]
reverse String
intRev
        h :: (Integer, Integer) -> Maybe (Char, (Integer, Integer))
h (Integer
f',Integer
err) = let (Integer
q,Integer
r) = Integer -> Integer -> (Integer, Integer)
forall a. Integral a => a -> a -> (a, a)
quotRem ((Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
base)Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
f') Integer
b
                         err' :: Integer
err' = (Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
base)Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
err
                         z :: (Integer, Integer)
z = (Integer
r, Integer
err')
                     in if Integer
err' Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
r Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
b
                        then (Char, (Integer, Integer)) -> Maybe (Char, (Integer, Integer))
forall a. a -> Maybe a
Just (Int -> Char
intToDigit (Integer -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
q), (Integer, Integer)
z)
                        else if Integer
err' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
min Integer
r Integer
b
                             then (Char, (Integer, Integer)) -> Maybe (Char, (Integer, Integer))
forall a. a -> Maybe a
Just (Int -> Char
intToDigit ((Integer -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
q Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`rem` (Int -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
base)), (Integer, Integer)
z)
                             else Maybe (Char, (Integer, Integer))
forall a. Maybe a
Nothing
        frac :: String
frac = ((Integer, Integer) -> Maybe (Char, (Integer, Integer)))
-> (Integer, Integer) -> String
forall b a. (b -> Maybe (a, b)) -> b -> [a]
unfoldr (Integer, Integer) -> Maybe (Char, (Integer, Integer))
h (Integer
f,Integer
e)
    in String
int String -> ShowS
forall a. [a] -> [a] -> [a]
++ if Bool
noFrac
              then String
""
              else String
"." String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
frac String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"~"

{-|
    Give the "bound on midpoint bit-size" component of an 'Approx'.

    The midpoint coponent should always be bounded by this as follows:
    @ abs m <= 2^mb@.
-}
mBound :: Approx -> Int
mBound :: Approx -> Int
mBound (Approx Int
mb Integer
_ Integer
_ Int
_) = Int
mb
mBound Approx
Bottom = String -> Int
forall a. HasCallStack => String -> a
error String
"mBound Bottom"

approxAutoMB :: Integer -> Integer -> Int -> Approx
approxAutoMB :: Integer -> Integer -> Int -> Approx
approxAutoMB Integer
m Integer
e Int
s = Int -> Integer -> Integer -> Int -> Approx
Approx Int
mb Integer
m Integer
e Int
s
    where
    ame :: Integer
ame = (Integer -> Integer
forall a. Num a => a -> a
abs Integer
m) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
e
    mb :: Int
mb | Integer
ame Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
4 = Int
2
       | Bool
otherwise = Int
1 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Integer -> Int
integerLog2 (Integer
ame Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
1)


mapMB :: (Int -> Int) -> Approx -> Approx
mapMB :: (Int -> Int) -> Approx -> Approx
mapMB Int -> Int
f (Approx Int
mb Integer
m Integer
e Int
s) = Int -> Integer -> Integer -> Int -> Approx
approxMB (Int -> Int
f Int
mb) Integer
m Integer
e Int
s
mapMB Int -> Int
_f Approx
Bottom = Approx
Bottom

setMB :: Int -> Approx -> Approx
setMB :: Int -> Approx -> Approx
setMB Int
mb = (Int -> Int) -> Approx -> Approx
mapMB (Int -> Int -> Int
forall a b. a -> b -> a
const Int
mb)


approxMB :: Int -> Integer -> Integer -> Int -> Approx
approxMB :: Int -> Integer -> Integer -> Int -> Approx
approxMB Int
mb Integer
m Integer
e Int
s = 
    Approx -> Approx
enforceMB (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Int -> Integer -> Integer -> Int -> Approx
Approx Int
mb Integer
m Integer
e Int
s

approxMB2 :: Int -> Int -> Integer -> Integer -> Int -> Approx
approxMB2 :: Int -> Int -> Integer -> Integer -> Int -> Approx
approxMB2 Int
mb1 Int
mb2 Integer
m Integer
e Int
s = 
    Approx -> Approx
enforceMB (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Int -> Integer -> Integer -> Int -> Approx
Approx (Int
mb1 Int -> Int -> Int
forall a. Ord a => a -> a -> a
`max` Int
mb2) Integer
m Integer
e Int
s

enforceMB :: Approx -> Approx
enforceMB :: Approx -> Approx
enforceMB Approx
Bottom = Approx
Bottom
enforceMB a :: Approx
a@(Approx Int
mb Integer
m Integer
e Int
s)
    | Int
m_size Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
mb = Approx
a
    | Integer -> Integer
forall a. Num a => a -> a
abs Integer
m Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
1 = Approx
a
    | Bool
otherwise = Int -> Integer -> Integer -> Int -> Approx
Approx Int
mb Integer
m' Integer
e'' (Int
s Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
d)
    where
    m_size :: Int
m_size = Int
1Int -> Int -> Int
forall a. Num a => a -> a -> a
+Integer -> Int
integerLog2 (Integer -> Integer
forall a. Num a => a -> a
abs Integer
m Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
1) --- |m| <= 2^m_size
    d :: Int
d = Int
m_size Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
mb
    m' :: Integer
m' = Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftR Integer
m Int
d -- we have: m' * 2^d <= m
    e' :: Integer
e' = Integer
1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ (Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftR (Integer
eInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
1) Int
d) -- we have: 0 <= e <= e' * 2^d
    e'' :: Integer
e''| Integer
m Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
m' Int
d  = Integer
e' -- no loss of information
       | Bool
otherwise = Integer
1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
e'


-- |Construct a centred approximation from the end-points.
endToApprox :: Int -> Extended Dyadic -> Extended Dyadic -> Approx
endToApprox :: Int -> Extended Dyadic -> Extended Dyadic -> Approx
endToApprox Int
mb (Finite Dyadic
l) (Finite Dyadic
u)
  | Dyadic
u Dyadic -> Dyadic -> Bool
forall a. Ord a => a -> a -> Bool
< Dyadic
l = Approx
Bottom -- Might be better with a signalling error.
  | Bool
otherwise =
    let a :: Dyadic
a@(Integer
m:^Int
s) = Dyadic -> Int -> Dyadic
forall a. Scalable a => a -> Int -> a
scale (Dyadic
l Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
+ Dyadic
u) (-Int
1)
        (Integer
n:^Int
t)   = Dyadic
uDyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
-Dyadic
a
        r :: Int
r        = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
s Int
t
        m' :: Integer
m'       = Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
m (Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
r)
        n' :: Integer
n'       = Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
n (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
r)
    in (Int -> Integer -> Integer -> Int -> Approx
approxMB Int
mb Integer
m' Integer
n' Int
r)
endToApprox Int
_ Extended Dyadic
_ Extended Dyadic
_ = Approx
Bottom

-- Interval operations
-- |Gives the lower bound of an approximation as an 'Extended' 'Dyadic' number.
lowerBound :: Approx -> Extended Dyadic
lowerBound :: Approx -> Extended Dyadic
lowerBound (Approx Int
_ Integer
m Integer
e Int
s) = Dyadic -> Extended Dyadic
forall a. a -> Extended a
Finite ((Integer
mInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
e)Integer -> Int -> Dyadic
:^Int
s)
lowerBound Approx
Bottom = Extended Dyadic
forall a. Extended a
NegInf

-- |Gives the upper bound of an approximation as an 'Extended' 'Dyadic' number.
upperBound :: Approx -> Extended Dyadic
upperBound :: Approx -> Extended Dyadic
upperBound (Approx Int
_ Integer
m Integer
e Int
s) = Dyadic -> Extended Dyadic
forall a. a -> Extended a
Finite ((Integer
mInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
e)Integer -> Int -> Dyadic
:^Int
s)
upperBound Approx
Bottom = Extended Dyadic
forall a. Extended a
PosInf

-- |Gives the lower bound of an 'Approx' as an exact 'Approx'.
lowerA :: Approx -> Approx
lowerA :: Approx -> Approx
lowerA Approx
Bottom = Approx
Bottom
lowerA (Approx Int
mb Integer
m Integer
e Int
s) = Int -> Integer -> Integer -> Int -> Approx
Approx Int
mb (Integer
mInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
e) Integer
0 Int
s

-- |Gives the upper bound of an 'Approx' as an exact 'Approx'.
upperA :: Approx -> Approx
upperA :: Approx -> Approx
upperA Approx
Bottom = Approx
Bottom
upperA (Approx Int
mb Integer
m Integer
e Int
s) = Int -> Integer -> Integer -> Int -> Approx
Approx Int
mb (Integer
mInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
e) Integer
0 Int
s

-- |Gives the mid-point of an approximation as a 'Maybe' 'Dyadic' number.
centre :: Approx -> Maybe Dyadic
centre :: Approx -> Maybe Dyadic
centre (Approx Int
_ Integer
m Integer
_ Int
s) = Dyadic -> Maybe Dyadic
forall a. a -> Maybe a
Just (Integer
mInteger -> Int -> Dyadic
:^Int
s)
centre Approx
_ = Maybe Dyadic
forall a. Maybe a
Nothing

-- |Gives the centre of an 'Approx' as an exact 'Approx'.
centreA :: Approx -> Approx
centreA :: Approx -> Approx
centreA Approx
Bottom = Approx
Bottom
centreA (Approx Int
mb Integer
m Integer
_e Int
s) = Int -> Integer -> Integer -> Int -> Approx
Approx Int
mb Integer
m Integer
0 Int
s

-- |Gives the radius of an approximation as a 'Dyadic' number. Currently a
-- partial function. Should be made to return an 'Extended' 'Dyadic'.
radius :: Approx -> Extended Dyadic
radius :: Approx -> Extended Dyadic
radius (Approx Int
_ Integer
_ Integer
e Int
s) = Dyadic -> Extended Dyadic
forall a. a -> Extended a
Finite (Integer
eInteger -> Int -> Dyadic
:^Int
s)
radius Approx
_ = Extended Dyadic
forall a. Extended a
PosInf

-- |Gives the lower bound of an approximation as an 'Extended' 'Dyadic' number.
diameter :: Approx -> Extended Dyadic
diameter :: Approx -> Extended Dyadic
diameter (Approx Int
_ Integer
_ Integer
e Int
s) = Dyadic -> Extended Dyadic
forall a. a -> Extended a
Finite (Dyadic -> Extended Dyadic) -> Dyadic -> Extended Dyadic
forall a b. (a -> b) -> a -> b
$ Dyadic
2 Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
* (Integer
eInteger -> Int -> Dyadic
:^Int
s)
diameter Approx
_ = Extended Dyadic
forall a. Extended a
PosInf

-- |Returns 'True' if the approximation is exact, i.e., it's diameter is 0.
exact :: Approx -> Bool
exact :: Approx -> Bool
exact (Approx Int
_ Integer
_ Integer
0 Int
_) = Bool
True
exact Approx
_ = Bool
False

-- |Checks if a number is approximated by an approximation, i.e., if it
-- belongs to the interval encoded by the approximation.
approximatedBy :: Real a => a -> Approx -> Bool
a
_ approximatedBy :: a -> Approx -> Bool
`approximatedBy` Approx
Bottom = Bool
True
a
r `approximatedBy` Approx
d =
    let r' :: Rational
r' = a -> Rational
forall a. Real a => a -> Rational
toRational a
r
    in Extended Dyadic -> Rational
forall a. Real a => a -> Rational
toRational (Approx -> Extended Dyadic
lowerBound Approx
d) Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
<= Rational
r' Bool -> Bool -> Bool
&& Rational
r' Rational -> Rational -> Bool
forall a. Ord a => a -> a -> Bool
<= Extended Dyadic -> Rational
forall a. Real a => a -> Rational
toRational (Approx -> Extended Dyadic
upperBound Approx
d)

-- |A partial order on approximations. The first approximation is better than
-- the second if it is a sub-interval of the second.
better :: Approx -> Approx -> Bool
Approx
d better :: Approx -> Approx -> Bool
`better` Approx
e = Approx -> Extended Dyadic
lowerBound Approx
d Extended Dyadic -> Extended Dyadic -> Bool
forall a. Ord a => a -> a -> Bool
>= Approx -> Extended Dyadic
lowerBound Approx
e Bool -> Bool -> Bool
&&
               Approx -> Extended Dyadic
upperBound Approx
d Extended Dyadic -> Extended Dyadic -> Bool
forall a. Ord a => a -> a -> Bool
<= Approx -> Extended Dyadic
upperBound Approx
e

-- |Turns a 'Dyadic' number into an exact approximation.
fromDyadic :: Dyadic -> Approx
fromDyadic :: Dyadic -> Approx
fromDyadic (Integer
m:^Int
s) = Integer -> Integer -> Int -> Approx
approxAutoMB Integer
m Integer
0 Int
s

-- |Turns a 'Dyadic' number into an exact approximation.
fromDyadicMB :: Int -> Dyadic -> Approx
fromDyadicMB :: Int -> Dyadic -> Approx
fromDyadicMB Int
mb (Integer
m:^Int
s) = Int -> Integer -> Integer -> Int -> Approx
approxMB Int
mb Integer
m Integer
0 Int
s

-- |Two approximations are equal if they encode the same interval.
instance Eq Approx where
    (Approx Int
_ Integer
m Integer
e Int
s) == :: Approx -> Approx -> Bool
== (Approx Int
_ Integer
n Integer
f Int
t)
        | Int
s Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
t = let k :: Int
k = Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
t
                   in Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
m Int
k Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
n Bool -> Bool -> Bool
&& Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
e Int
k Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
f
        | Int
s Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<  Int
t = let k :: Int
k = Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
s
                   in Integer
m Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
n Int
k Bool -> Bool -> Bool
&& Integer
e Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
f Int
k
    Approx
Bottom == Approx
Bottom = Bool
True
    Approx
_ == Approx
_ = Bool
False

-- |Not a sensible instance. Just used to allow to allow enumerating integers
-- using \'..\' notation.
instance Enum Approx where
    toEnum :: Int -> Approx
toEnum Int
n = Int -> Integer -> Integer -> Int -> Approx
Approx Int
64 (Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n) Integer
0 Int
0
    fromEnum :: Approx -> Int
fromEnum (Approx Int
_ Integer
m Integer
_ Int
s) = Integer -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Integer -> Int) -> Integer -> Int
forall a b. (a -> b) -> a -> b
$ Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
shift Integer
m Int
s
    fromEnum Approx
Bottom = Int
0

instance Num Approx where
    (Approx Int
mb1 Integer
m Integer
e Int
s) + :: Approx -> Approx -> Approx
+ (Approx Int
mb2 Integer
n Integer
f Int
t)
        | Int
s Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
t = let k :: Int
k = Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
t
                   in Int -> Int -> Integer -> Integer -> Int -> Approx
approxMB2 Int
mb1 Int
mb2 (Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
m Int
k Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
n) (Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
e Int
k Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
f) Int
t
        | Int
s Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<  Int
t = let k :: Int
k = Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
s
                   in Int -> Int -> Integer -> Integer -> Int -> Approx
approxMB2 Int
mb1 Int
mb2 (Integer
m Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
n Int
k) (Integer
e Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
f Int
k) Int
s
    Approx
_ + Approx
_ = Approx
Bottom
    (Approx Int
mb1 Integer
m Integer
e Int
s) * :: Approx -> Approx -> Approx
* (Approx Int
mb2 Integer
n Integer
f Int
t)
        | Integer
am Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>= Integer
e Bool -> Bool -> Bool
&& Integer
an Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>= Integer
f Bool -> Bool -> Bool
&& Integer
a Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
0           = Int -> Int -> Integer -> Integer -> Int -> Approx
approxMB2 Int
mb1 Int
mb2 (Integer
aInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
d) (Integer
abInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
ac) Int
u
        | Integer
am Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>= Integer
e Bool -> Bool -> Bool
&& Integer
an Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>= Integer
f Bool -> Bool -> Bool
&& Integer
a Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0           = Int -> Int -> Integer -> Integer -> Int -> Approx
approxMB2 Int
mb1 Int
mb2 (Integer
aInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
d) (Integer
abInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
ac) Int
u
        | Integer
am Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
e Bool -> Bool -> Bool
&& Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>= Integer
f                      = Int -> Int -> Integer -> Integer -> Int -> Approx
approxMB2 Int
mb1 Int
mb2 (Integer
aInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
b) (Integer
acInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
d) Int
u
        | Integer
am Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
e Bool -> Bool -> Bool
&& -Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>= Integer
f                     = Int -> Int -> Integer -> Integer -> Int -> Approx
approxMB2 Int
mb1 Int
mb2 (Integer
aInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
b) (Integer
acInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
d) Int
u
        | Integer
m Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>= Integer
e Bool -> Bool -> Bool
&& Integer
an Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
f                      = Int -> Int -> Integer -> Integer -> Int -> Approx
approxMB2 Int
mb1 Int
mb2 (Integer
aInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
c) (Integer
abInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
d) Int
u
        | -Integer
m Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>= Integer
e Bool -> Bool -> Bool
&& Integer
an Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
f                     = Int -> Int -> Integer -> Integer -> Int -> Approx
approxMB2 Int
mb1 Int
mb2 (Integer
aInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
c) (Integer
abInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
d) Int
u
        | Integer
a Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0                                = Int -> Int -> Integer -> Integer -> Int -> Approx
approxMB2 Int
mb1 Int
mb2 (Integer
0) (Integer
abInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
acInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
d) Int
u
        | Integer
am Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
e Bool -> Bool -> Bool
&& Integer
an Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
f Bool -> Bool -> Bool
&& Integer
a Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
0 Bool -> Bool -> Bool
&& Integer
ab Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
ac  = Int -> Int -> Integer -> Integer -> Int -> Approx
approxMB2 Int
mb1 Int
mb2 (Integer
aInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
ac) (Integer
abInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
d) Int
u
        | Integer
am Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
e Bool -> Bool -> Bool
&& Integer
an Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
f Bool -> Bool -> Bool
&& Integer
a Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
0 Bool -> Bool -> Bool
&& Integer
ab Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
ac = Int -> Int -> Integer -> Integer -> Int -> Approx
approxMB2 Int
mb1 Int
mb2 (Integer
aInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
ab) (Integer
acInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
d) Int
u
        | Integer
am Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
e Bool -> Bool -> Bool
&& Integer
an Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
f Bool -> Bool -> Bool
&& Integer
a Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 Bool -> Bool -> Bool
&& Integer
ab Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
ac  = Int -> Int -> Integer -> Integer -> Int -> Approx
approxMB2 Int
mb1 Int
mb2 (Integer
aInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
ac) (Integer
abInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
d) Int
u
        | Integer
am Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
e Bool -> Bool -> Bool
&& Integer
an Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
f Bool -> Bool -> Bool
&& Integer
a Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 Bool -> Bool -> Bool
&& Integer
ab Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
ac = Int -> Int -> Integer -> Integer -> Int -> Approx
approxMB2 Int
mb1 Int
mb2 (Integer
aInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
ab) (Integer
acInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
d) Int
u
      where am :: Integer
am = (Integer -> Integer
forall a. Num a => a -> a
abs Integer
m)
            an :: Integer
an = (Integer -> Integer
forall a. Num a => a -> a
abs Integer
n)
            a :: Integer
a = Integer
mInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
n
            b :: Integer
b = Integer
mInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
f
            c :: Integer
c = Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
e
            d :: Integer
d = Integer
eInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
f
            ab :: Integer
ab = (Integer -> Integer
forall a. Num a => a -> a
abs Integer
b)
            ac :: Integer
ac = (Integer -> Integer
forall a. Num a => a -> a
abs Integer
c)
            u :: Int
u = Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
t
    Approx
_ * Approx
_ = Approx
Bottom
    negate :: Approx -> Approx
negate (Approx Int
mb Integer
m Integer
e Int
s) = Int -> Integer -> Integer -> Int -> Approx
Approx Int
mb (-Integer
m) Integer
e Int
s
    negate Approx
Bottom = Approx
Bottom
    abs :: Approx -> Approx
abs (Approx Int
mb Integer
m Integer
e Int
s)
        | Integer
m' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
e    = let e' :: Integer
e' = Integer
m'Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
e
                      in Int -> Integer -> Integer -> Int -> Approx
Approx Int
mb Integer
e' Integer
e' (Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)
        | Bool
otherwise = Int -> Integer -> Integer -> Int -> Approx
Approx Int
mb Integer
m' Integer
e Int
s
      where m' :: Integer
m' = Integer -> Integer
forall a. Num a => a -> a
abs Integer
m
    abs Approx
Bottom = Approx
Bottom
    signum :: Approx -> Approx
signum (Approx Int
_ Integer
m Integer
e Int
_)
        | Integer
e Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 = Int -> Integer -> Integer -> Int -> Approx
Approx Int
64 (Integer -> Integer
forall a. Num a => a -> a
signum Integer
m) Integer
0 Int
0
        | Integer -> Integer
forall a. Num a => a -> a
abs Integer
m Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
e = Int -> Integer -> Integer -> Int -> Approx
Approx Int
64 Integer
0 Integer
1 Int
0
        | Integer -> Integer
forall a. Num a => a -> a
abs Integer
m Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
e = Int -> Integer -> Integer -> Int -> Approx
Approx Int
64 (Integer -> Integer
forall a. Num a => a -> a
signum Integer
m) Integer
1 (-Int
1)
        | Bool
otherwise = Int -> Integer -> Integer -> Int -> Approx
Approx Int
64 (Integer -> Integer
forall a. Num a => a -> a
signum Integer
m) Integer
0 Int
0
    signum Approx
Bottom = Int -> Integer -> Integer -> Int -> Approx
Approx Int
64 Integer
0 Integer
1 Int
0
    fromInteger :: Integer -> Approx
fromInteger Integer
i = (Int -> Int) -> Approx -> Approx
mapMB (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
64) (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Integer -> Integer -> Int -> Approx
approxAutoMB Integer
i Integer
0 Int
0

-- |Convert a rational number into an approximation of that number with
-- 'Precision' bits correct after the binary point.
toApprox :: Precision -> Rational -> Approx
toApprox :: Int -> Rational -> Approx
toApprox Int
t Rational
r = Integer -> Integer -> Int -> Approx
approxAutoMB Integer
m Integer
e (-Int
t Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)
    where
    m :: Integer
m = (Integer
2 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
r_scaled_rounded)
    r_scaled_rounded :: Integer
r_scaled_rounded = Rational -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
round Rational
r_scaled
    r_scaled :: Rational
r_scaled = Rational
rRational -> Rational -> Rational
forall a. Num a => a -> a -> a
*Rational
2Rational -> Int -> Rational
forall a b. (Fractional a, Integral b) => a -> b -> a
^^Int
t
    e :: Integer
e | Rational
r_scaled Rational -> Rational -> Bool
forall a. Eq a => a -> a -> Bool
== Integer -> Rational
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
r_scaled_rounded = Integer
0
      | Bool
otherwise = Integer
1

-- |Convert a rational number into an approximation of that number with
-- 'mBound' significant bits correct.
toApproxMB :: Int -> Rational -> Approx
toApproxMB :: Int -> Rational -> Approx
toApproxMB Int
mb Rational
r = 
    (Int -> Integer -> Integer -> Int -> Approx
Approx Int
mb (Rational -> Integer
forall a. Ratio a -> a
numerator Rational
r) Integer
0 Int
0) Approx -> Approx -> Approx
forall a. Fractional a => a -> a -> a
/ (Int -> Integer -> Integer -> Int -> Approx
Approx Int
mb (Rational -> Integer
forall a. Ratio a -> a
denominator Rational
r) Integer
0 Int
0)

-- |Not a proper Fractional type as Approx are intervals.
instance Fractional Approx where
    fromRational :: Rational -> Approx
fromRational = Int -> Rational -> Approx
toApprox Int
defaultPrecision
    recip :: Approx -> Approx
recip = Approx -> Approx
recipA
    Approx
a1 / :: Approx -> Approx -> Approx
/ Approx
a2 = Approx
a1 Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* (Approx -> Approx
recipA (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Int -> Approx -> Approx
setMB Int
mb Approx
a2)
        where
        mb :: Int
mb = Approx -> Int
mBound Approx
a1 Int -> Int -> Int
forall a. Ord a => a -> a -> a
`max` Approx -> Int
mBound Approx
a2

-- |Compute the reciprocal of an approximation. The number of bits after the
-- binary point is bounded by the 'midpoint bound' if the input is exact.
-- Otherwise, a good approximation with essentially the same significance as
-- the input will be computed.
recipA :: Approx -> Approx
recipA :: Approx -> Approx
recipA Approx
Bottom = Approx
Bottom
recipA (Approx Int
mb Integer
m Integer
e Int
s)
    | Integer
e Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 Bool -> Bool -> Bool
&& Integer
m Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= Integer
0
                  = let s' :: Int
s' = Integer -> Int
integerLog2 (Integer -> Integer
forall a. Num a => a -> a
abs Integer
m)
                    in if Integer -> Integer
forall a. Num a => a -> a
abs Integer
m Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Int -> Integer
forall a. Bits a => Int -> a
bit Int
s'
                       then
                            Int -> Integer -> Integer -> Int -> Approx
Approx Int
mb (Integer -> Integer
forall a. Num a => a -> a
signum Integer
m) Integer
0 (-Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
s')
                       else
                            Int -> Integer -> Integer -> Int -> Approx
Approx Int
mb
                            (Rational -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
round (Int -> Integer
forall a. Bits a => Int -> a
bit (Int
mbInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
s') Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% Integer
m))
                            Integer
1
                            (-Int
mbInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
s')
    | (Integer -> Integer
forall a. Num a => a -> a
abs Integer
m) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
e = let d :: Integer
d = Integer
mInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
mInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
eInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
e
                        d2 :: Integer
d2 = Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftR Integer
d Int
1
                        s' :: Int
s' = Integer -> Int
integerLog2 Integer
d Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
errorBits
                    in Approx -> Approx
boundErrorTerm (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Int -> Integer -> Integer -> Int -> Approx
approxMB Int
mb
                           ((Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
m Int
s' Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
d2) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div` Integer
d)
                           ((Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
e Int
s' Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
d2Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
1) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div` Integer
d Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1)
                           (-Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
s')
    --  (abs m) > e = let d = m*m-e*e
    --                     s' = 2 * (integerLog2 m + errorBits)
    --                 in boundErrorTerm $ Approx
    --                        (round (unsafeShiftL m s'%(d)))
    --                        (ceiling (1%2 + unsafeShiftL e s'%(d)))
    --                        (-s-s')
    | Bool
otherwise   = Approx
Bottom

-- |Divide an approximation by an integer.
divAInteger :: Approx -> Integer -> Approx
divAInteger :: Approx -> Integer -> Approx
divAInteger Approx
Bottom Integer
_ = Approx
Bottom
divAInteger (Approx Int
mb Integer
m Integer
e Int
s) Integer
n =
  let t :: Int
t = Integer -> Int
integerLog2 Integer
n
  in Int -> Integer -> Integer -> Int -> Approx
approxMB Int
mb
             (Rational -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
round (Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
m Int
t Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% Integer
n))
             (Rational -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
ceiling (Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
e Int
t Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% Integer
n))
             Int
s

-- |Compute the modulus of two approximations.
modA :: Approx -> Approx -> Approx
modA :: Approx -> Approx -> Approx
modA (Approx Int
mb1 Integer
m Integer
e Int
s) (Approx Int
mb2 Integer
n Integer
f Int
t) =
    let r :: Int
r = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
s Int
t
        (Integer
d,Integer
m') = Integer -> Integer -> (Integer, Integer)
forall a. Integral a => a -> a -> (a, a)
divMod (Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
m (Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
r)) (Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
n (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
r))
        e' :: Integer
e' = Integer -> Int -> Integer
forall a. Scalable a => a -> Int -> a
scale Integer
e (Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
r) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer -> Integer
forall a. Num a => a -> a
abs Integer
d Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer -> Int -> Integer
forall a. Scalable a => a -> Int -> a
scale Integer
f (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
r)
    in Int -> Int -> Integer -> Integer -> Int -> Approx
approxMB2 Int
mb1 Int
mb2 Integer
m' Integer
e' Int
r
modA Approx
_ Approx
_ = Approx
Bottom

-- |Compute the integer quotient (although returned as an 'Approx' since it
-- may be necessary to return 'Bottom' if the integer quotient can't be
-- determined) and the modulus as an approximation of two approximations.
divModA :: Approx -> Approx -> (Approx, Approx)
divModA :: Approx -> Approx -> (Approx, Approx)
divModA (Approx Int
mb1 Integer
m Integer
e Int
s) (Approx Int
mb2 Integer
n Integer
f Int
t) =
    let r :: Int
r = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
s Int
t
        (Integer
d,Integer
m') = Integer -> Integer -> (Integer, Integer)
forall a. Integral a => a -> a -> (a, a)
divMod (Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
m (Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
r)) (Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
n (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
r))
        e' :: Integer
e' = Integer
e Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer -> Integer
forall a. Num a => a -> a
abs Integer
d Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
f
    in (Integer -> Approx
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
d, Int -> Int -> Integer -> Integer -> Int -> Approx
approxMB2 Int
mb1 Int
mb2 Integer
m' Integer
e' Int
r)
divModA Approx
_ Approx
_ = (Approx
Bottom, Approx
Bottom)

-- |Not a proper Ord type as Approx are intervals.
instance Ord Approx where
    compare :: Approx -> Approx -> Ordering
compare (Approx Int
_ Integer
m Integer
e Int
s) (Approx Int
_ Integer
n Integer
f Int
t)
        | Integer
e Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 Bool -> Bool -> Bool
&& Integer
f Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 = Dyadic -> Dyadic -> Ordering
forall a. Ord a => a -> a -> Ordering
compare (Integer
mInteger -> Int -> Dyadic
:^Int
s) (Integer
nInteger -> Int -> Dyadic
:^Int
t)
        | Dyadic -> Dyadic
forall a. Num a => a -> a
abs ((Integer
mInteger -> Int -> Dyadic
:^Int
s)Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
-(Integer
nInteger -> Int -> Dyadic
:^Int
t)) Dyadic -> Dyadic -> Bool
forall a. Ord a => a -> a -> Bool
> (Integer
eInteger -> Int -> Dyadic
:^Int
s)Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
+(Integer
fInteger -> Int -> Dyadic
:^Int
t) = Dyadic -> Dyadic -> Ordering
forall a. Ord a => a -> a -> Ordering
compare (Integer
mInteger -> Int -> Dyadic
:^Int
s) (Integer
nInteger -> Int -> Dyadic
:^Int
t)
        | Bool
otherwise                           = String -> Ordering
forall a. HasCallStack => String -> a
error String
"compare: comparisons are partial on Approx"
    compare Approx
_ Approx
_ = String -> Ordering
forall a. HasCallStack => String -> a
error String
"compare: comparisons are partial on Approx"

-- |The 'toRational' function is partial since there is no good rational
-- number to return for the trivial approximation 'Bottom'.
--
-- Note that converting to a rational number will give only a single rational
-- point. Do not expect to be able to recover the interval from this value.
instance Real Approx where
    toRational :: Approx -> Rational
toRational (Approx Int
_ Integer
m Integer
e Int
s) = Rational -> Rational -> Rational
forall a. RealFrac a => a -> a -> Rational
approxRational
                                  (Dyadic -> Rational
forall a. Real a => a -> Rational
toRational (Integer
mInteger -> Int -> Dyadic
:^Int
s))
                                  (Dyadic -> Rational
forall a. Real a => a -> Rational
toRational (Integer
eInteger -> Int -> Dyadic
:^Int
s))
    toRational Approx
_ = Rational
forall a. HasCallStack => a
undefined

-- |Convert the centre of an approximation into a 'Maybe' 'Double'.
toDoubleA :: Approx -> Maybe Double
toDoubleA :: Approx -> Maybe Double
toDoubleA = (Dyadic -> Double) -> Maybe Dyadic -> Maybe Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (Rational -> Double
forall a. Fractional a => Rational -> a
fromRational (Rational -> Double) -> (Dyadic -> Rational) -> Dyadic -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Dyadic -> Rational
forall a. Real a => a -> Rational
toRational) (Maybe Dyadic -> Maybe Double)
-> (Approx -> Maybe Dyadic) -> Approx -> Maybe Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Approx -> Maybe Dyadic
centre


-- |Computes the precision of an approximation. This is roughly the number of
-- correct bits after the binary point.
precision :: Approx -> Extended Precision
precision :: Approx -> Extended Int
precision (Approx Int
_ Integer
_ Integer
0 Int
_) = Extended Int
forall a. Extended a
PosInf
precision (Approx Int
_ Integer
_ Integer
e Int
s) = Int -> Extended Int
forall a. a -> Extended a
Finite (Int -> Extended Int) -> Int -> Extended Int
forall a b. (a -> b) -> a -> b
$ - Int
s Int -> Int -> Int
forall a. Num a => a -> a -> a
- (Integer -> Int
integerLog2 Integer
e) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1
precision Approx
Bottom         = Extended Int
forall a. Extended a
NegInf

-- |Computes the significance of an approximation. This is roughly the number
-- of significant bits.
significance :: Approx -> Extended Int
significance :: Approx -> Extended Int
significance (Approx Int
_ Integer
_ Integer
0 Int
_) = Extended Int
forall a. Extended a
PosInf
significance (Approx Int
_ Integer
0 Integer
_ Int
_) = Extended Int
forall a. Extended a
NegInf
significance (Approx Int
_ Integer
m Integer
1 Int
_) =  Int -> Extended Int
forall a. a -> Extended a
Finite (Int -> Extended Int) -> Int -> Extended Int
forall a b. (a -> b) -> a -> b
$ Integer -> Int
integerLog2 (Integer -> Integer
forall a. Num a => a -> a
abs Integer
m) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1
significance (Approx Int
_ Integer
m Integer
e Int
_) =
    Int -> Extended Int
forall a. a -> Extended a
Finite (Int -> Extended Int) -> Int -> Extended Int
forall a b. (a -> b) -> a -> b
$ (Integer -> Int
integerLog2 (Integer -> Integer
forall a. Num a => a -> a
abs Integer
m)) Int -> Int -> Int
forall a. Num a => a -> a -> a
- (Integer -> Int
integerLog2 (Integer
eInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
1)) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1
significance Approx
Bottom         = Extended Int
forall a. Extended a
NegInf

{-|
This function bounds the error term of an 'Approx'.

It is always the case that @x `'better'` boundErrorTerm x@.

Consider an approximation @Approx mb m e s@. If @e@ has /k/ bits then that
essentially expresses that the last /k/ bits of @m@ are unknown or garbage. By
scaling both @m@ and @e@ so that @e@ has a small number of bits we save on
memory space and computational effort to compute operations. On the other
hand, if we remove too many bits in this way, the shift in the mid-point of the
interval becomes noticable and it may adversely affect convergence speed of
computations. The number of bits allowed for @e@ after the operation is
determined by the constant 'errorBits'.

== Domain interpretation and verification

For this implementation to be correct it is required that this function is
below the identity on the domain /D/ of 'Approx'. For efficiency it is
desirable to be as close to the identity as possible.

This function will map a converging sequence to a converging sequence.
-}
boundErrorTerm :: Approx -> Approx
boundErrorTerm :: Approx -> Approx
boundErrorTerm Approx
Bottom = Approx
Bottom
boundErrorTerm a :: Approx
a@(Approx Int
mb Integer
m Integer
e Int
s)
    | Integer
e Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
errorBound = Approx
a
    | Bool
otherwise =
        let k :: Int
k = Integer -> Int
integerLog2 Integer
e Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
errorBits
            t :: Bool
t = Integer -> Int -> Bool
forall a. Bits a => a -> Int -> Bool
testBit Integer
m (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)
            m' :: Integer
m' = Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftR Integer
m Int
k
            -- may overflow and use errorBits+1
            e' :: Integer
e' = Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftR (Integer
e Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Int -> Integer
forall a. Bits a => Int -> a
bit (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)) Int
k Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1
        in if Bool
t
           then Int -> Integer -> Integer -> Int -> Approx
Approx Int
mb (Integer
m'Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
1) Integer
e' (Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
k)
           else Int -> Integer -> Integer -> Int -> Approx
Approx Int
mb Integer
m'     Integer
e' (Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
k)

boundErrorTermMB :: Approx -> Approx
boundErrorTermMB :: Approx -> Approx
boundErrorTermMB Approx
Bottom = Approx
Bottom
boundErrorTermMB a :: Approx
a@(Approx Int
_ Integer
m Integer
e Int
s)
    | Integer
e Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
errorBoundMB = Approx
a
    | Bool
otherwise =
        let k :: Int
k = Integer -> Int
integerLog2 Integer
e Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
errorBitsMB
            t :: Bool
t = Integer -> Int -> Bool
forall a. Bits a => a -> Int -> Bool
testBit Integer
m (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)
            m' :: Integer
m' = Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftR Integer
m Int
k
            -- may overflow and use errorBits+1
            e' :: Integer
e' = Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftR (Integer
e Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Int -> Integer
forall a. Bits a => Int -> a
bit (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)) Int
k Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1
        in if Bool
t
           then Integer -> Integer -> Int -> Approx
approxAutoMB (Integer
m'Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
1) Integer
e' (Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
k)
           else Integer -> Integer -> Int -> Approx
approxAutoMB Integer
m'     Integer
e' (Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
k)

{-|
Limits the size of an approximation by restricting how much precision an
approximation can have.

It is always the case that @x `'better'` limitSize x@.

This is accomplished by restricting the exponent of the approximation from
below. In other words, we limit the precision possible.

It is conceivable to limit the significance of an approximation rather than
the precision. This would be an interesting research topic.

== Domain interpretation and verification

For this implementation to be correct it is required that this function is
below the identity on the domain /D/ of 'Approx'. For efficiency it is
desirable to be as close to the identity as possible.

This function will NOT map a converging sequence to a converging sequence for
a fixed precision argument. However, if the function is applied with
increasing precision for a converging sequence, then this will give a
converging sequence.
-}
limitSize :: Precision -> Approx -> Approx
limitSize :: Int -> Approx -> Approx
limitSize Int
_ Approx
Bottom = Approx
Bottom
limitSize Int
l a :: Approx
a@(Approx Int
mb Integer
m Integer
e Int
s)
    | Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0     = Int -> Integer -> Integer -> Int -> Approx
Approx (Int
mbInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
k)
                  ((if Integer -> Int -> Bool
forall a. Bits a => a -> Int -> Bool
testBit Integer
m (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) then (Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
1) else Integer -> Integer
forall a. a -> a
id) (Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftR Integer
m Int
k))
                  (Integer
1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ (Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftR (Integer
e Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Int -> Integer
forall a. Bits a => Int -> a
bit (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)) Int
k))
                  (-Int
l)
    | Bool
otherwise = Approx
a
    where k :: Int
k = (-Int
s)Int -> Int -> Int
forall a. Num a => a -> a -> a
-Int
l

-- |Throws an exception if the precision of an approximation is not larger
-- than the deafult minimum.
checkPrecisionLeft :: Approx -> Approx
checkPrecisionLeft :: Approx -> Approx
checkPrecisionLeft Approx
a
        | Approx -> Extended Int
precision Approx
a Extended Int -> Extended Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int -> Extended Int
forall (f :: * -> *) a. Applicative f => a -> f a
pure Int
defaultPrecision = Approx
a
        | Bool
otherwise = ArithException -> Approx
forall a e. Exception e => e -> a
throw (ArithException -> Approx) -> ArithException -> Approx
forall a b. (a -> b) -> a -> b
$ ArithException
LossOfPrecision

-- |Bounds the error term and limits the precision of an approximation.
--
-- It is always the case that @x `'better'` limitAndBound x@.
limitAndBound :: Precision -> Approx -> Approx
limitAndBound :: Int -> Approx -> Approx
limitAndBound Int
limit =
    Int -> Approx -> Approx
limitSize Int
limit (Approx -> Approx) -> (Approx -> Approx) -> Approx -> Approx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Approx -> Approx
boundErrorTerm

limitAndBoundMB :: Precision -> Approx -> Approx
limitAndBoundMB :: Int -> Approx -> Approx
limitAndBoundMB Int
limit =
    Int -> Approx -> Approx
limitSize Int
limit (Approx -> Approx) -> (Approx -> Approx) -> Approx -> Approx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Approx -> Approx
boundErrorTermMB

-- | Find the hull of two approximations.
unionA :: Approx -> Approx -> Approx
unionA :: Approx -> Approx -> Approx
unionA Approx
Bottom Approx
_ = Approx
Bottom
unionA Approx
_ Approx
Bottom = Approx
Bottom
unionA a :: Approx
a@(Approx Int
mb1 Integer
_ Integer
_ Int
_) b :: Approx
b@(Approx Int
mb2 Integer
_ Integer
_ Int
_) =
    Int -> Extended Dyadic -> Extended Dyadic -> Approx
endToApprox (Int
mb1 Int -> Int -> Int
forall a. Ord a => a -> a -> a
`max` Int
mb2) (Approx -> Extended Dyadic
lowerBound Approx
a Extended Dyadic -> Extended Dyadic -> Extended Dyadic
forall a. Ord a => a -> a -> a
`min` Approx -> Extended Dyadic
lowerBound Approx
b) (Approx -> Extended Dyadic
upperBound Approx
a Extended Dyadic -> Extended Dyadic -> Extended Dyadic
forall a. Ord a => a -> a -> a
`max` Approx -> Extended Dyadic
upperBound Approx
b)

-- | Find the intersection of two approximations.
intersectionA :: Approx -> Approx -> Approx
intersectionA :: Approx -> Approx -> Approx
intersectionA Approx
Bottom Approx
a = Approx
a
intersectionA Approx
a Approx
Bottom = Approx
a
intersectionA a :: Approx
a@(Approx Int
mb1 Integer
_ Integer
_ Int
_) b :: Approx
b@(Approx Int
mb2 Integer
_ Integer
_ Int
_) =
  if Extended Dyadic
l Extended Dyadic -> Extended Dyadic -> Bool
forall a. Ord a => a -> a -> Bool
<= Extended Dyadic
u
    then Int -> Extended Dyadic -> Extended Dyadic -> Approx
endToApprox (Int
mb1 Int -> Int -> Int
forall a. Ord a => a -> a -> a
`max` Int
mb2) Extended Dyadic
l Extended Dyadic
u
    else String -> Approx
forall a. HasCallStack => String -> a
error String
"Trying to take intersection of two non-overlapping Approx."
  where l :: Extended Dyadic
l = (Approx -> Extended Dyadic
lowerBound Approx
a Extended Dyadic -> Extended Dyadic -> Extended Dyadic
forall a. Ord a => a -> a -> a
`max` Approx -> Extended Dyadic
lowerBound Approx
b)
        u :: Extended Dyadic
u = (Approx -> Extended Dyadic
upperBound Approx
a Extended Dyadic -> Extended Dyadic -> Extended Dyadic
forall a. Ord a => a -> a -> a
`min` Approx -> Extended Dyadic
upperBound Approx
b)

-- | Determine if two approximations overlap.
consistentA :: Approx -> Approx -> Bool
consistentA :: Approx -> Approx -> Bool
consistentA Approx
Bottom Approx
_ = Bool
True
consistentA Approx
_ Approx
Bottom = Bool
True
consistentA Approx
a Approx
b = (Approx -> Extended Dyadic
lowerBound Approx
a Extended Dyadic -> Extended Dyadic -> Extended Dyadic
forall a. Ord a => a -> a -> a
`max` Approx -> Extended Dyadic
lowerBound Approx
b) Extended Dyadic -> Extended Dyadic -> Bool
forall a. Ord a => a -> a -> Bool
<= (Approx -> Extended Dyadic
upperBound Approx
a Extended Dyadic -> Extended Dyadic -> Extended Dyadic
forall a. Ord a => a -> a -> a
`min` Approx -> Extended Dyadic
upperBound Approx
b)

-- |Given a list of polynom coefficients and a value this evaluates the
-- polynomial at that value.
--
-- This works correctly only for exact coefficients.
--
-- Should give a tighter bound on the result since we reduce the dependency
-- problem.
poly :: [Approx] -> Approx -> Approx
poly :: [Approx] -> Approx -> Approx
poly [] Approx
_ = Approx
0
poly [Approx]
_ Approx
Bottom = Approx
Bottom
poly [Approx]
as x :: Approx
x@(Approx Int
mb Integer
_ Integer
_ Int
_) =
    let --poly' :: [Dyadic] -> Dyadic -> Dyadic
        poly' :: [c] -> c -> c
poly' [c]
as' c
x' = [c] -> c
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([c] -> c) -> ([c] -> [c]) -> [c] -> c
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (c -> c -> c) -> [c] -> [c] -> [c]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith c -> c -> c
forall a. Num a => a -> a -> a
(*) [c]
as' ([c] -> c) -> [c] -> c
forall a b. (a -> b) -> a -> b
$ c -> [c]
forall a. Num a => a -> [a]
pow c
x'
        ms :: [Dyadic]
ms = (Approx -> Dyadic) -> [Approx] -> [Dyadic]
forall a b. (a -> b) -> [a] -> [b]
map ((Dyadic -> (Dyadic -> Dyadic) -> Maybe Dyadic -> Dyadic
forall b a. b -> (a -> b) -> Maybe a -> b
maybe (String -> Dyadic
forall a. HasCallStack => String -> a
error String
"Can't compute poly with Bottom coefficients") Dyadic -> Dyadic
forall a. a -> a
id) (Maybe Dyadic -> Dyadic)
-> (Approx -> Maybe Dyadic) -> Approx -> Dyadic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Approx -> Maybe Dyadic
centre) [Approx]
as
        (Just Dyadic
c) = Approx -> Maybe Dyadic
centre Approx
x
        (Integer
m':^Int
s) = [Dyadic] -> Dyadic -> Dyadic
forall c. Num c => [c] -> c -> c
poly' [Dyadic]
ms Dyadic
c
        ds :: [Approx]
ds = (Approx -> Approx -> Approx) -> [Approx] -> [Approx] -> [Approx]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
(*) ([Approx] -> [Approx]
forall a. [a] -> [a]
tail [Approx]
as) ((Int -> Approx) -> [Int] -> [Approx]
forall a b. (a -> b) -> [a] -> [b]
map Int -> Approx
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Int
1,Int
2..] :: [Int]))
        (Finite Dyadic
b) = Approx -> Extended Dyadic
upperBound (Approx -> Extended Dyadic)
-> (Approx -> Approx) -> Approx -> Extended Dyadic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Approx -> Approx
forall a. Num a => a -> a
abs (Approx -> Extended Dyadic) -> Approx -> Extended Dyadic
forall a b. (a -> b) -> a -> b
$ [Approx] -> Approx -> Approx
forall c. Num c => [c] -> c -> c
poly' [Approx]
ds Approx
x
        (Finite (Integer
e':^Int
_)) = (Dyadic -> Dyadic) -> Extended Dyadic -> Extended Dyadic
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (Dyadic
bDyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
*) (Extended Dyadic -> Extended Dyadic)
-> Extended Dyadic -> Extended Dyadic
forall a b. (a -> b) -> a -> b
$ Approx -> Extended Dyadic
radius Approx
x
        -- exponent above will be same as s
    in Int -> Integer -> Integer -> Int -> Approx
approxMB Int
mb Integer
m' Integer
e' Int
s

-- |Gives a list of powers of a number, i.e., [1,x,x^2,...].
pow :: (Num a) => a -> [a]
pow :: a -> [a]
pow a
x = (a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate (a -> a -> a
forall a. Num a => a -> a -> a
* a
x) a
1

-- |Computes lists of binomial coefficients. [[1], [1,1], [1,2,1], [1,3,3,1], ...]
binomialCoefficients :: (Num a) => [[a]]
binomialCoefficients :: [[a]]
binomialCoefficients =
    let f :: [a] -> [a]
f [a]
ss = a
1 a -> [a] -> [a]
forall a. a -> [a] -> [a]
: (a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Num a => a -> a -> a
(+) [a]
ss ([a] -> [a]
forall a. [a] -> [a]
tail [a]
ss) [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a
1]
    in ([a] -> [a]) -> [a] -> [[a]]
forall a. (a -> a) -> a -> [a]
iterate [a] -> [a]
forall a. Num a => [a] -> [a]
f [a
1]

-- |Computes powers of approximations. Should give tighter intervals than the
-- general 'pow' since take the dependency problem into account. However, so
-- far benchmarking seems to indicate that the cost is too high, but this may
-- depend on the application.
powers :: Approx -> [Approx]
powers :: Approx -> [Approx]
powers (Approx Int
mb Integer
m Integer
e Int
s) =
    let ms :: [Integer]
ms = Integer -> [Integer]
forall a. Num a => a -> [a]
pow Integer
m
        es :: [Integer]
es = Integer -> [Integer]
forall a. Num a => a -> [a]
pow Integer
e
        f :: [Integer] -> [Integer]
f = [Integer] -> [Integer]
forall a. [a] -> [a]
reverse ([Integer] -> [Integer])
-> ([Integer] -> [Integer]) -> [Integer] -> [Integer]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer -> Integer -> Integer)
-> [Integer] -> [Integer] -> [Integer]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
(*) [Integer]
ms ([Integer] -> [Integer])
-> ([Integer] -> [Integer]) -> [Integer] -> [Integer]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Integer] -> [Integer]
forall a. [a] -> [a]
reverse ([Integer] -> [Integer])
-> ([Integer] -> [Integer]) -> [Integer] -> [Integer]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer -> Integer -> Integer)
-> [Integer] -> [Integer] -> [Integer]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
(*) [Integer]
es
        sumAlt :: [b] -> (b, b)
sumAlt [] = (b
0,b
0)
        sumAlt (b
x:[]) = (b
x,b
0)
        sumAlt (b
x:b
y:[b]
xs) = let (b
a,b
b) = [b] -> (b, b)
sumAlt [b]
xs in (b
ab -> b -> b
forall a. Num a => a -> a -> a
+b
x,b
bb -> b -> b
forall a. Num a => a -> a -> a
+b
y)
        g :: Int -> (Integer, Integer) -> Approx
g Int
s' (Integer
m', Integer
e') = Int -> Integer -> Integer -> Int -> Approx
approxMB Int
mb Integer
m' Integer
e' Int
s'
    in (Int -> (Integer, Integer) -> Approx)
-> [Int] -> [(Integer, Integer)] -> [Approx]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Int -> (Integer, Integer) -> Approx
g ((Int -> Int) -> Int -> [Int]
forall a. (a -> a) -> a -> [a]
iterate (Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
s) Int
0) ([(Integer, Integer)] -> [Approx])
-> [(Integer, Integer)] -> [Approx]
forall a b. (a -> b) -> a -> b
$ ([Integer] -> (Integer, Integer))
-> [[Integer]] -> [(Integer, Integer)]
forall a b. (a -> b) -> [a] -> [b]
map ([Integer] -> (Integer, Integer)
forall b. Num b => [b] -> (b, b)
sumAlt ([Integer] -> (Integer, Integer))
-> ([Integer] -> [Integer]) -> [Integer] -> (Integer, Integer)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Integer] -> [Integer]
f) [[Integer]]
forall a. Num a => [[a]]
binomialCoefficients
powers Approx
_ = Approx -> [Approx]
forall a. a -> [a]
repeat Approx
Bottom

{-|
Old implementation of sqrt using Heron's method. Using the current method
below proved to be more than twice as fast for small arguments (~50 bits) and
many times faster for larger arguments.
-}
sqrtHeronA :: Precision -> Approx -> Approx
sqrtHeronA :: Int -> Approx -> Approx
sqrtHeronA Int
_ Approx
Bottom = Approx
Bottom
sqrtHeronA Int
k a :: Approx
a@(Approx Int
mb Integer
m Integer
e Int
s)
    | -Integer
m Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
e    = String -> Approx
forall a. HasCallStack => String -> a
error String
"Attempting sqrt of Approx containing only negative numbers."
    | Integer
m Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
e     = Approx
Bottom
    | Integer
e Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0    = let (Integer
n:^Int
t) = Int -> Dyadic -> Dyadic
shiftD (-Int
k) (Dyadic -> Dyadic) -> Dyadic -> Dyadic
forall a b. (a -> b) -> a -> b
$ Int -> Dyadic -> Dyadic
sqrtD (-Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2) (Integer
mInteger -> Int -> Dyadic
:^Int
s)
                  in (Int -> Int) -> Approx -> Approx
mapMB (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
mb) (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Integer -> Integer -> Int -> Approx
approxAutoMB Integer
n Integer
1 Int
t
    | Integer
m Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
e    = let (Integer
n:^Int
t) = Int -> Dyadic -> Dyadic
sqrtD (Int
s Int -> Int -> Int
forall a. Integral a => a -> a -> a
`quot` Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
-Int
errorBits) ((Integer
mInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
e)Integer -> Int -> Dyadic
:^Int
s)
                      n' :: Integer
n' = (Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
2) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`quot` Integer
2
                  in Int -> Integer -> Integer -> Int -> Approx
approxMB Int
mb Integer
n' Integer
n' Int
t
    | Bool
otherwise = let (Finite Int
p) = Approx -> Extended Int
significance Approx
a
                      s' :: Int
s' = Int
s Int -> Int -> Int
forall a. Integral a => a -> a -> a
`quot` Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
errorBits
                      l :: Dyadic
l@(Integer
n:^Int
t) = Int -> Dyadic -> Dyadic
sqrtD Int
s' ((Integer
mInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
e)Integer -> Int -> Dyadic
:^Int
s)
                      (Integer
n':^Int
t') = Int -> Dyadic -> Dyadic -> Dyadic
sqrtD' Int
s' ((Integer
mInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
e)Integer -> Int -> Dyadic
:^Int
s) Dyadic
l
                  in Int -> Extended Dyadic -> Extended Dyadic -> Approx
endToApprox Int
mb (Dyadic -> Extended Dyadic
forall a. a -> Extended a
Finite ((Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
1)Integer -> Int -> Dyadic
:^Int
t)) (Dyadic -> Extended Dyadic
forall a. a -> Extended a
Finite ((Integer
n'Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
1)Integer -> Int -> Dyadic
:^Int
t'))

{-|
Compute the square root of an approximation.

This and many other operations on approximations is just a reimplementation of
interval arithmetic, with an extra argument limiting the effort put into the
computation. This is done via the precision argument.

The resulting approximation should approximate the image of every point in the
input approximation.
-}
sqrtA :: Approx -> Approx
sqrtA :: Approx -> Approx
sqrtA Approx
Bottom = Approx
Bottom
sqrtA x :: Approx
x@(Approx Int
_ Integer
0 Integer
0 Int
_) =  Approx
x
sqrtA x :: Approx
x@(Approx Int
mb Integer
_ Integer
_ Int
_) 
    | Approx -> Approx
upperA Approx
x Approx -> Approx -> Bool
forall a. Ord a => a -> a -> Bool
< Approx
1 = 
        Int -> Approx -> Approx
sqrtRecA Int
k (Approx -> Approx
recipA (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Int -> Approx -> Approx
setMB Int
k Approx
x)
    | Bool
otherwise =
        -- limitAndBoundMB mb $ 
        Approx
x Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Int -> Approx -> Approx
sqrtRecA Int
k Approx
x
    where
    k :: Int
k = Int
mb Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2

{-|
This uses Newton's method for computing the reciprocal of the square root.
-}
sqrtRecA :: Precision -> Approx -> Approx
sqrtRecA :: Int -> Approx -> Approx
sqrtRecA Int
_ Approx
Bottom = Approx
Bottom
sqrtRecA Int
k a :: Approx
a@(Approx Int
mb Integer
m Integer
e Int
s)
  | -Integer
m Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
e    = String -> Approx
forall a. HasCallStack => String -> a
error String
"Attempting sqrtRec of Approx containing only negative numbers."
  | Integer
m Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
e     = Approx
Bottom
  | Integer
e Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0    = let (Integer
n:^Int
t) = Int -> Dyadic -> Dyadic
shiftD (-Int
k) (Dyadic -> Dyadic) -> Dyadic -> Dyadic
forall a b. (a -> b) -> a -> b
$ Int -> Dyadic -> Dyadic
sqrtRecD (-Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2) (Integer
mInteger -> Int -> Dyadic
:^Int
s)
                in (Int -> Int) -> Approx -> Approx
mapMB (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
mb) (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Integer -> Integer -> Int -> Approx
approxAutoMB Integer
n Integer
1 Int
t
  | Integer
m Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
e    = let (Integer
n:^Int
t) = Int -> Dyadic -> Dyadic
sqrtRecD (Int
s Int -> Int -> Int
forall a. Integral a => a -> a -> a
`quot` Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
-Int
errorBits) ((Integer
mInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
e)Integer -> Int -> Dyadic
:^Int
s)
                    n' :: Integer
n' = (Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
2) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`quot` Integer
2
                in Int -> Integer -> Integer -> Int -> Approx
approxMB Int
mb Integer
n' Integer
n' Int
t
  | Bool
otherwise = let (Finite Int
p) = Approx -> Extended Int
significance Approx
a
                    s' :: Int
s' = Int
s Int -> Int -> Int
forall a. Integral a => a -> a -> a
`quot` Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
errorBits
                    (Integer
n:^Int
t) = Int -> Dyadic -> Dyadic
sqrtRecD Int
s' ((Integer
mInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
e)Integer -> Int -> Dyadic
:^Int
s) -- upper bound of result
                    -- We have tried to use sqrtRecD' with the above value as
                    -- a first approximation to the result. However, the low
                    -- endpoint may be too far away as a starting value to
                    -- ensure convergence to the right root. It's possible
                    -- that if we swap the order we would be fine. But as it
                    -- is, this computes a new first approximation.
                    (Integer
n':^Int
t') = Int -> Dyadic -> Dyadic
sqrtRecD Int
s' ((Integer
mInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
e)Integer -> Int -> Dyadic
:^Int
s) -- lower bound of result
                in Int -> Extended Dyadic -> Extended Dyadic -> Approx
endToApprox Int
mb (Dyadic -> Extended Dyadic
forall a. a -> Extended a
Finite ((Integer
n'Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
1)Integer -> Int -> Dyadic
:^Int
t')) (Dyadic -> Extended Dyadic
forall a. a -> Extended a
Finite ((Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
1)Integer -> Int -> Dyadic
:^Int
t))


{-|
The starting values for newton iterations can be found using the auxiliary function findStartingValues below.

For example, to generate the starting values for sqrtRecD above using three leading bits for both odd and even exponents the following was used:

> findStartingValues ((1/) . sqrt) [1,1.125..2]
[Approx 4172150648 1 (-32),Approx 3945434766 1 (-32),Approx 3752147976 1 (-32),Approx 3584793264 1 (-32),Approx 3438037830 1 (-32),Approx 3307969824 1 (-32),Approx 3191645366 1 (-32),Approx 3086800564 1 (-32)]
> mapM_ (putStrLn . showInBaseA 2 . limitSize 6) it
0.111110~
0.111011~
0.111000~
0.110101~
0.110011~
0.110001~
0.110000~
0.101110~
> findStartingValues ((1/) . sqrt) [2,2.25..4]
[Approx 2950156016 1 (-32),Approx 2789843678 1 (-32),Approx 2653169278 1 (-32),Approx 2534831626 1 (-32),Approx 2431059864 1 (-32),Approx 2339087894 1 (-32),Approx 2256834080 1 (-32),Approx 2182697612 1 (-32)]
> mapM_ (putStrLn . showInBaseA 2 . limitSize 6) it
0.101100~
0.101010~
0.101000~
0.100110~
0.100100~
0.100011~
0.100010~
0.100001~
-}
findStartingValues :: (Double -> Double) -> [Double] -> [Approx]
findStartingValues :: (Double -> Double) -> [Double] -> [Approx]
findStartingValues Double -> Double
f = (Double -> Approx) -> [Double] -> [Approx]
forall a b. (a -> b) -> [a] -> [b]
map (Rational -> Approx
forall a. Fractional a => Rational -> a
fromRational (Rational -> Approx) -> (Double -> Rational) -> Double -> Approx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Rational
forall a. Real a => a -> Rational
toRational (Double -> Rational) -> (Double -> Double) -> Double -> Rational
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2)) ([Double] -> [Approx])
-> ([Double] -> [Double]) -> [Double] -> [Approx]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (\[Double]
l -> (Double -> Double -> Double) -> [Double] -> [Double] -> [Double]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Double -> Double -> Double
forall a. Num a => a -> a -> a
(+) [Double]
l ([Double] -> [Double]
forall a. [a] -> [a]
tail [Double]
l)) ([Double] -> [Double])
-> ([Double] -> [Double]) -> [Double] -> [Double]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double -> Double) -> [Double] -> [Double]
forall a b. (a -> b) -> [a] -> [b]
map Double -> Double
f

-- |Square an approximation. Gives the exact image interval, as opposed to
-- multiplicating a number with itself which will give a slightly larger
-- interval due to the dependency problem.
sqrA :: Approx -> Approx
sqrA :: Approx -> Approx
sqrA Approx
Bottom = Approx
Bottom
sqrA (Approx Int
mb Integer
m Integer
e Int
s)
  | Integer
am Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
e = Int -> Integer -> Integer -> Int -> Approx
approxMB Int
mb (Integer
mInteger -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^(Int
2 :: Int) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
eInteger -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^(Int
2 :: Int)) (Integer
2Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
amInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
e) (Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
s)
  | Bool
otherwise = let m' :: Integer
m' = (Integer
am Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
e)Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^(Int
2 :: Int) in Int -> Integer -> Integer -> Int -> Approx
approxMB Int
mb Integer
m' Integer
m' (Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)
  where am :: Integer
am = Integer -> Integer
forall a. Num a => a -> a
abs Integer
m
-- Binary splitting

{-|
Binary splitting summation of linearly convergent series as described in
/'Fast multiprecision evaluation of series of rational numbers'/ by B Haible
and T Papanikolaou, ANTS-III Proceedings of the Third International Symposium
on Algorithmic Number Theory Pages 338-350, 1998.

The main idea is to balance the computations so that more operations are
performed with values of similar size. Using the underlying fast
multiplication algorithms this will give performance benefits.

The algorithm parallelises well. However, a final division is needed at the
end to compute /T\/BQ/ which amount to a substantial portion of the
computation time.
-}
abpq :: Num a => [Integer] -> [Integer] -> [a] -> [a] -> Int -> Int -> (a, a, Integer, a)
abpq :: [Integer]
-> [Integer] -> [a] -> [a] -> Int -> Int -> (a, a, Integer, a)
abpq [Integer]
as [Integer]
bs [a]
ps [a]
qs Int
n1 Int
n2
    | Int
n Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
1 = ([a]
ps [a] -> Int -> a
forall a. [a] -> Int -> a
!! Int
n1, [a]
qs [a] -> Int -> a
forall a. [a] -> Int -> a
!! Int
n1, [Integer]
bs [Integer] -> Int -> Integer
forall a. [a] -> Int -> a
!! Int
n1, Integer -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral ([Integer]
as [Integer] -> Int -> Integer
forall a. [a] -> Int -> a
!! Int
n1) a -> a -> a
forall a. Num a => a -> a -> a
* [a]
ps [a] -> Int -> a
forall a. [a] -> Int -> a
!! Int
n1)
    | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
6  = let as' :: [Integer]
as' = Int -> [Integer] -> [Integer]
forall a. Int -> [a] -> [a]
take Int
n ([Integer] -> [Integer]) -> [Integer] -> [Integer]
forall a b. (a -> b) -> a -> b
$ Int -> [Integer] -> [Integer]
forall a. Int -> [a] -> [a]
drop Int
n1 [Integer]
as
                   bs' :: [Integer]
bs' = Int -> [Integer] -> [Integer]
forall a. Int -> [a] -> [a]
take Int
n ([Integer] -> [Integer]) -> [Integer] -> [Integer]
forall a b. (a -> b) -> a -> b
$ Int -> [Integer] -> [Integer]
forall a. Int -> [a] -> [a]
drop Int
n1 [Integer]
bs
                   ps' :: [a]
ps' = Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
n ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop Int
n1 [a]
ps
                   qs' :: [a]
qs' = Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
take Int
n ([a] -> [a]) -> [a] -> [a]
forall a b. (a -> b) -> a -> b
$ Int -> [a] -> [a]
forall a. Int -> [a] -> [a]
drop Int
n1 [a]
qs
                   pbs :: Integer
pbs = [Integer] -> Integer
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [Integer]
bs'
                   bs'' :: [Integer]
bs'' = (Integer -> Integer) -> [Integer] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map (Integer
pbs Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div`) [Integer]
bs'
                   ps'' :: [a]
ps'' = (a -> a -> a) -> [a] -> [a]
forall a. (a -> a -> a) -> [a] -> [a]
scanl1 a -> a -> a
forall a. Num a => a -> a -> a
(*) [a]
ps'
                   qs'' :: [a]
qs'' = (a -> a -> a) -> [a] -> [a]
forall a. (a -> a -> a) -> [a] -> [a]
scanr1 a -> a -> a
forall a. Num a => a -> a -> a
(*) ([a] -> [a]
forall a. [a] -> [a]
tail [a]
qs' [a] -> [a] -> [a]
forall a. [a] -> [a] -> [a]
++ [a
1])
               in ([a]
ps'' [a] -> Int -> a
forall a. [a] -> Int -> a
!! (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1), [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [a]
qs', Integer
pbs
                  , [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([a] -> a) -> [a] -> a
forall a b. (a -> b) -> a -> b
$ (Integer -> Integer -> a -> a -> a)
-> [Integer] -> [Integer] -> [a] -> [a] -> [a]
forall a b c d e.
(a -> b -> c -> d -> e) -> [a] -> [b] -> [c] -> [d] -> [e]
zipWith4 (\Integer
a Integer
b a
p a
q -> Integer -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
a a -> a -> a
forall a. Num a => a -> a -> a
* Integer -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
b a -> a -> a
forall a. Num a => a -> a -> a
* a
p a -> a -> a
forall a. Num a => a -> a -> a
* a
q)
                                   [Integer]
as' [Integer]
bs'' [a]
ps'' [a]
qs'')
    | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
1  =
        let (a
pl, a
ql, Integer
bl, a
tl) = [Integer]
-> [Integer] -> [a] -> [a] -> Int -> Int -> (a, a, Integer, a)
forall a.
Num a =>
[Integer]
-> [Integer] -> [a] -> [a] -> Int -> Int -> (a, a, Integer, a)
abpq [Integer]
as [Integer]
bs [a]
ps [a]
qs Int
n1 Int
m
            (a
pr, a
qr, Integer
br, a
tr) = [Integer]
-> [Integer] -> [a] -> [a] -> Int -> Int -> (a, a, Integer, a)
forall a.
Num a =>
[Integer]
-> [Integer] -> [a] -> [a] -> Int -> Int -> (a, a, Integer, a)
abpq [Integer]
as [Integer]
bs [a]
ps [a]
qs Int
m Int
n2
        in (a
pl a -> a -> a
forall a. Num a => a -> a -> a
* a
pr, a
ql a -> a -> a
forall a. Num a => a -> a -> a
* a
qr, Integer
bl Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
br, Integer -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
br a -> a -> a
forall a. Num a => a -> a -> a
* a
qr a -> a -> a
forall a. Num a => a -> a -> a
* a
tl a -> a -> a
forall a. Num a => a -> a -> a
+ Integer -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
bl a -> a -> a
forall a. Num a => a -> a -> a
* a
pl a -> a -> a
forall a. Num a => a -> a -> a
* a
tr)
    | Bool
otherwise = String -> (a, a, Integer, a)
forall a. HasCallStack => String -> a
error String
"Non-expected case in binary splitting"
  where
    n :: Int
n = Int
n2 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
n1
    m :: Int
m = (Int
n1 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
n2 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2

ones :: Num a => [a]
ones :: [a]
ones = a -> [a]
forall a. a -> [a]
repeat a
1

{-|
Computes the list [lg 0!, lg 1!, lg 2!, ...].
-}
-- To be changed to Stirling formula if that is faster
log2Factorials :: [Int]
log2Factorials :: [Int]
log2Factorials = (Integer -> Int) -> [Integer] -> [Int]
forall a b. (a -> b) -> [a] -> [b]
map Integer -> Int
integerLog2 ([Integer] -> [Int])
-> ([Integer] -> [Integer]) -> [Integer] -> [Int]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Integer -> Integer -> Integer) -> [Integer] -> [Integer]
forall a. (a -> a -> a) -> [a] -> [a]
scanl1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
(*) ([Integer] -> [Int]) -> [Integer] -> [Int]
forall a b. (a -> b) -> a -> b
$ Integer
1Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:[Integer
1..]

-- Straighforward Taylor summation

{-|
Computes a sum of the form ∑ aₙ/qₙ where aₙ are approximations and qₙ are
integers. Terms are added as long as they are larger than the current
precision bound. The sum is adjusted for the tail of the series. For this to
be correct we need the the terms to converge geometrically to 0 by a factor of
at least 2.
-}
taylor :: Precision -> [Approx] -> [Integer] -> Approx
taylor :: Int -> [Approx] -> [Integer] -> Approx
taylor Int
res [Approx]
as [Integer]
qs =
  let res' :: Int
res' = Int
res Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
errorBits
      f :: Approx -> Integer -> Approx
f Approx
a Integer
q = Int -> Approx -> Approx
limitAndBound Int
res' (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Approx
a Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Approx -> Approx
recipA (Int -> Approx -> Approx
setMB (Approx -> Int
mBound Approx
a) (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Integer -> Approx
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
q)
      mb :: [Approx]
mb = (Approx -> Integer -> Approx) -> [Approx] -> [Integer] -> [Approx]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Approx -> Integer -> Approx
f [Approx]
as [Integer]
qs
      ([Approx]
cs,(Approx
d:[Approx]
_)) = (Approx -> Bool) -> [Approx] -> ([Approx], [Approx])
forall a. (a -> Bool) -> [a] -> ([a], [a])
span Approx -> Bool
nonZeroCentredA [Approx]
mb -- This span and the sum on the next line do probably not fuse!
  in Approx -> Approx -> Approx
fudge ([Approx] -> Approx
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Approx]
cs) Approx
d

-- | A list of factorial values.
fac :: Num a => [a]
fac :: [a]
fac = (Integer -> a) -> [Integer] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map Integer -> a
forall a. Num a => Integer -> a
fromInteger ([Integer] -> [a]) -> [Integer] -> [a]
forall a b. (a -> b) -> a -> b
$ Integer
1 Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
: (Integer -> Integer -> Integer) -> [Integer] -> [Integer]
forall a. (a -> a -> a) -> [a] -> [a]
scanl1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
(*) [Integer
1..]

-- | A list of the factorial values of odd numbers.
oddFac :: Num a => [a]
oddFac :: [a]
oddFac = let f :: [a] -> [a]
f (a
_:a
x:[a]
xs) = a
xa -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a] -> [a]
f [a]
xs
             f [a]
_ = String -> [a]
forall a. HasCallStack => String -> a
error String
"Impossible"
         in [a] -> [a]
forall a. [a] -> [a]
f [a]
forall a. Num a => [a]
fac

{-
evenFac :: Num a => [a]
evenFac = let f (x:_:xs) = x:f xs
              f _ = error "Impossible"
          in f fac
-}

-- | Checks if the centre of an approximation is not 0.
nonZeroCentredA :: Approx -> Bool
nonZeroCentredA :: Approx -> Bool
nonZeroCentredA Approx
Bottom = Bool
False
nonZeroCentredA (Approx Int
_ Integer
0 Integer
_ Int
_) = Bool
False
nonZeroCentredA Approx
_ = Bool
True

-- This version is faster especially far smaller precision.

{-|
Computes the sum of the form ∑ aₙxⁿ where aₙ and x are approximations.

Terms are added as long as they are larger than the current precision bound.
The sum is adjusted for the tail of the series. For this to be correct we need
the the terms to converge geometrically to 0 by a factor of at least 2.
-}
taylorA :: Precision -> [Approx] -> Approx -> Approx
taylorA :: Int -> [Approx] -> Approx -> Approx
taylorA Int
res [Approx]
as Approx
x =
  Approx -> Approx -> Approx
fudge Approx
sm Approx
d
  where
  (Approx
sm, Approx
d) = [(Approx, Approx)] -> (Approx, Approx)
forall b. [(Approx, b)] -> (Approx, b)
sumAndNext ([(Approx, Approx)] -> (Approx, Approx))
-> ([Approx] -> [(Approx, Approx)]) -> [Approx] -> (Approx, Approx)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Approx, Approx) -> Bool)
-> [(Approx, Approx)] -> [(Approx, Approx)]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Approx -> Bool
nonZeroCentredA (Approx -> Bool)
-> ((Approx, Approx) -> Approx) -> (Approx, Approx) -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Approx, Approx) -> Approx
forall a b. (a, b) -> a
fst) ([(Approx, Approx)] -> [(Approx, Approx)])
-> ([Approx] -> [(Approx, Approx)])
-> [Approx]
-> [(Approx, Approx)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Approx] -> [(Approx, Approx)]
forall a. [a] -> [(a, a)]
addNext ([Approx] -> [(Approx, Approx)])
-> ([Approx] -> [Approx]) -> [Approx] -> [(Approx, Approx)]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Approx -> Approx) -> [Approx] -> [Approx]
forall a b. (a -> b) -> [a] -> [b]
map (Int -> Approx -> Approx
limitAndBound Int
res) ([Approx] -> (Approx, Approx)) -> [Approx] -> (Approx, Approx)
forall a b. (a -> b) -> a -> b
$ (Approx -> Approx -> Approx) -> [Approx] -> [Approx] -> [Approx]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
(*) [Approx]
as (Approx -> [Approx]
forall a. Num a => a -> [a]
pow Approx
x)
  sumAndNext :: [(Approx, b)] -> (Approx, b)
sumAndNext = Approx -> [(Approx, b)] -> (Approx, b)
forall t b. Num t => t -> [(t, b)] -> (t, b)
aux Approx
0
    where
    aux :: t -> [(t, b)] -> (t, b)
aux t
a [(t
b,b
dd)] = (t
at -> t -> t
forall a. Num a => a -> a -> a
+t
b,b
dd)
    aux t
a ((t
b,b
_):[(t, b)]
rest) = t -> [(t, b)] -> (t, b)
aux (t
at -> t -> t
forall a. Num a => a -> a -> a
+t
b) [(t, b)]
rest
    aux t
_ [(t, b)]
_ = (t, b)
forall a. HasCallStack => a
undefined
  addNext :: [a] -> [(a, a)]
addNext (a
x1:a
x2:[a]
xs) = (a
x1,a
x2)(a, a) -> [(a, a)] -> [(a, a)]
forall a. a -> [a] -> [a]
:([a] -> [(a, a)]
addNext (a
x2a -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
xs))
  addNext [a]
_ = String -> [(a, a)]
forall a. HasCallStack => String -> a
error String
"taylorA: end of initite sequence" 

{- Exponential computed by standard Taylor expansion after range reduction.
-}

{-|
The exponential of an approximation. There are three implementation using
Taylor expansion here. This is just choosing between them.

More thorough benchmarking would be desirable.

Is faster for small approximations < ~2000 bits.
-}
expA :: Approx -> Approx
expA :: Approx -> Approx
expA = Approx -> Approx
expTaylorA'

-- | Exponential by binary splitting summation of Taylor series.
expBinarySplittingA :: Precision -> Approx -> Approx
expBinarySplittingA :: Int -> Approx -> Approx
expBinarySplittingA Int
_ Approx
Bottom = Approx
Bottom
expBinarySplittingA Int
res a :: Approx
a@(Approx Int
mb Integer
m Integer
e Int
s) =
  let s' :: Int
s' = Int
s Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Integer -> Int
integerLog2 Integer
m
      -- r' chosen so that a' below is smaller than 1/2
      r' :: Int
r' = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double -> Int) -> (Int -> Double) -> Int -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> (Int -> Double) -> Int -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Double) -> (Int -> Int) -> Int -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
5 (Int -> Int) -> Int -> Int
forall a b. (a -> b) -> a -> b
$ Int
res
      r :: Int
r = Int
s' Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
r'
      -- a' is a scaled by 2^k so that 2^(-r') <= a' < 2^(-r'+1)
      a' :: Approx
a' = Int -> Integer -> Integer -> Int -> Approx
Approx Int
mb Integer
m Integer
e (Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
r)
      mb' :: Int
mb' = Int
mbInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
res
      -- (Finite c) = min (significance a) (Finite res)
      (Just Int
n) = (Int -> Bool) -> [Int] -> Maybe Int
forall a. (a -> Bool) -> [a] -> Maybe Int
findIndex (Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
resInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
r) ([Int] -> Maybe Int) -> [Int] -> Maybe Int
forall a b. (a -> b) -> a -> b
$ (Int -> Int -> Int) -> [Int] -> [Int] -> [Int]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Int -> Int -> Int
forall a. Num a => a -> a -> a
(+) [Int]
log2Factorials [Int
0,Int
r'..]
      (Approx
p, Approx
q, Integer
b, Approx
t) = [Integer]
-> [Integer]
-> [Approx]
-> [Approx]
-> Int
-> Int
-> (Approx, Approx, Integer, Approx)
forall a.
Num a =>
[Integer]
-> [Integer] -> [a] -> [a] -> Int -> Int -> (a, a, Integer, a)
abpq [Integer]
forall a. Num a => [a]
ones
                          [Integer]
forall a. Num a => [a]
ones
                          ((Approx -> Approx) -> [Approx] -> [Approx]
forall a b. (a -> b) -> [a] -> [b]
map (Int -> Approx -> Approx
setMB Int
mb') ([Approx] -> [Approx]) -> [Approx] -> [Approx]
forall a b. (a -> b) -> a -> b
$ Approx
1Approx -> [Approx] -> [Approx]
forall a. a -> [a] -> [a]
:Approx -> [Approx]
forall a. a -> [a]
repeat Approx
a')
                          ((Approx -> Approx) -> [Approx] -> [Approx]
forall a b. (a -> b) -> [a] -> [b]
map (Int -> Approx -> Approx
setMB Int
mb') ([Approx] -> [Approx]) -> [Approx] -> [Approx]
forall a b. (a -> b) -> a -> b
$ Approx
1Approx -> [Approx] -> [Approx]
forall a. a -> [a] -> [a]
:[Approx
1..])
                          Int
0
                          Int
n
      nextTerm :: Approx
nextTerm = Approx
a Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Approx
p Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Approx -> Approx
recipA (Int -> Approx
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Approx
q)
      ss :: [Approx]
ss = (Approx -> Approx) -> Approx -> [Approx]
forall a. (a -> a) -> a -> [a]
iterate (Approx -> Approx
boundErrorTerm (Approx -> Approx) -> (Approx -> Approx) -> Approx -> Approx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Approx -> Approx
sqrA) (Approx -> [Approx]) -> Approx -> [Approx]
forall a b. (a -> b) -> a -> b
$ Approx -> Approx -> Approx
fudge (Approx
t Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Approx -> Approx
recipA (Integer -> Approx
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
bApprox -> Approx -> Approx
forall a. Num a => a -> a -> a
*Approx
q)) Approx
nextTerm
  in [Approx]
ss [Approx] -> Int -> Approx
forall a. [a] -> Int -> a
!! Int
r

-- | Exponential by summation of Taylor series.
expTaylorA :: Precision -> Approx -> Approx
expTaylorA :: Int -> Approx -> Approx
expTaylorA Int
_ Approx
Bottom = Approx
Bottom
expTaylorA Int
res (Approx Int
mb Integer
m Integer
e Int
s) =
  let s' :: Int
s' = Int
s Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Integer -> Int
integerLog2 Integer
m
      -- r' chosen so that a' below is smaller than 1/2
      r' :: Int
r' = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double -> Int) -> (Int -> Double) -> Int -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> (Int -> Double) -> Int -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Double) -> (Int -> Int) -> Int -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
5 (Int -> Int) -> Int -> Int
forall a b. (a -> b) -> a -> b
$ Int
res
      r :: Int
r = Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
0 (Int -> Int) -> Int -> Int
forall a b. (a -> b) -> a -> b
$ Int
s' Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
r'
      -- a' is a scaled by 2^k so that 2^(-r') <= a' < 2^(-r'+1)
      a' :: Approx
a' = (Int -> Integer -> Integer -> Int -> Approx
Approx (Int
mb Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
res) Integer
m Integer
e (Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
r))
      t :: Approx
t = Int -> [Approx] -> [Integer] -> Approx
taylor
            (Int
res Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
r)
            ((Approx -> Approx) -> Approx -> [Approx]
forall a. (a -> a) -> a -> [a]
iterate (Approx
a'Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
*) Approx
1)
            ((Integer -> Integer -> Integer) -> [Integer] -> [Integer]
forall a. (a -> a -> a) -> [a] -> [a]
scanl1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
(*) ([Integer] -> [Integer]) -> [Integer] -> [Integer]
forall a b. (a -> b) -> a -> b
$ Integer
1Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:[Integer
1..])
  in ([Approx] -> Int -> Approx
forall a. [a] -> Int -> a
!! Int
r) ([Approx] -> Approx) -> (Approx -> [Approx]) -> Approx -> Approx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Approx -> Approx) -> Approx -> [Approx]
forall a. (a -> a) -> a -> [a]
iterate (Approx -> Approx
boundErrorTermMB (Approx -> Approx) -> (Approx -> Approx) -> Approx -> Approx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Approx -> Approx
sqrA) (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Approx
t
   
-- | Exponential by summation of Taylor series.
expTaylorA' :: Approx -> Approx
expTaylorA' :: Approx -> Approx
expTaylorA' Approx
Bottom = Approx
Bottom
expTaylorA' Approx
a 
    | Approx -> Approx
upperA Approx
a Approx -> Approx -> Bool
forall a. Ord a => a -> a -> Bool
< Approx
0 = Approx -> Approx
recipA (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Approx -> Approx
aux (-Approx
a)
    | Bool
otherwise = Approx -> Approx
aux Approx
a
    where
    aux :: Approx -> Approx
aux Approx
Bottom = Approx
Bottom
    aux (Approx Int
mb Integer
0 Integer
0 Int
_) = Int -> Integer -> Integer -> Int -> Approx
Approx Int
mb Integer
1 Integer
0 Int
0
    aux (Approx Int
mb Integer
m Integer
0 Int
s) =
        let s' :: Int
s' = Int
s Int -> Int -> Int
forall a. Num a => a -> a -> a
+ (Integer -> Int
integerLog2 Integer
m)
            -- r' chosen so that a' below is smaller than 1/2
            r' :: Int
r' = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double -> Int) -> (Int -> Double) -> Int -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> (Int -> Double) -> Int -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Double) -> (Int -> Int) -> Int -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
5 (Int -> Int) -> Int -> Int
forall a b. (a -> b) -> a -> b
$ Int
mb
            r :: Int
r = Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
0 (Int -> Int) -> Int -> Int
forall a b. (a -> b) -> a -> b
$ Int
s' Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
r'
            mb'_ :: Int
mb'_ = Int
mb Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
r Int -> Int -> Int
forall a. Num a => a -> a -> a
+ (Integer -> Int
integerLog2 Integer
m) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1
            mb' :: Int
mb' = (Int
120Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
mb'_) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
100
            -- a' is a scaled by 2^k so that 2^(-r') <= a' < 2^(-r'+1)
            a' :: Approx
a' = (Int -> Integer -> Integer -> Int -> Approx
Approx Int
mb' Integer
m Integer
0 (Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
r))
            t :: Approx
t = Approx -> Approx
boundErrorTermMB (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Int -> [Approx] -> Approx -> Approx
taylorA Int
mb' ((Approx -> Approx) -> [Approx] -> [Approx]
forall a b. (a -> b) -> [a] -> [b]
map (Approx -> Approx
recipA (Approx -> Approx) -> (Approx -> Approx) -> Approx -> Approx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Approx -> Approx
setMB Int
mb') [Approx]
forall a. Num a => [a]
fac) Approx
a'
        in ([Approx] -> Int -> Approx
forall a. [a] -> Int -> a
!! Int
r) ([Approx] -> Approx) -> (Approx -> [Approx]) -> Approx -> Approx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Approx -> Approx) -> Approx -> [Approx]
forall a. (a -> a) -> a -> [a]
iterate (Approx -> Approx
boundErrorTermMB (Approx -> Approx) -> (Approx -> Approx) -> Approx -> Approx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Approx -> Approx
sqrA) (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Approx
t
    aux Approx
a2 = Approx -> Approx
aux (Approx -> Approx
lowerA Approx
a2) Approx -> Approx -> Approx
`unionA` Approx -> Approx
aux (Approx -> Approx
upperA Approx
a2)
   
{- Logarithms computed by ln x = 2*atanh ((x-1)/(x+1)) after range reduction.
-}

{-|

Computing the logarithm of an approximation. This chooses the fastest implementation.

More thorough benchmarking is desirable.

Binary splitting is faster than Taylor. AGM should be used over ~1000 bits.
-}
logA :: Approx -> Approx
-- This implementation asks for the dyadic approximation of the endpoints, we
-- should instead use that, after the first range reduction, the derivative is
-- less than 3/2 on the interval, so it easy to just compute one expensive
-- computation. We could even make use of the fact that the derivative on the
-- interval x is bounded by 1/x to get a tighter bound on the error.
logA :: Approx -> Approx
logA Approx
Bottom = Approx
Bottom
logA x :: Approx
x@(Approx Int
_ Integer
m Integer
e Int
_)
  | Integer
m Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
e Bool -> Bool -> Bool
&& Approx -> Approx
upperA Approx
x Approx -> Approx -> Bool
forall a. Ord a => a -> a -> Bool
< Approx
1 = -(Approx -> Approx
logInternal (Approx -> Approx
recipA Approx
x))
  | Integer
m Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
e = Approx -> Approx
logInternal Approx
x
--    let (n :^ t) = logD (negate p) $ (m-e) :^ s
--        (n' :^ t') = logD (negate p) $ (m+e) :^ s
--    in endToApprox (Finite ((n-1):^t)) (Finite ((n'+1):^t'))
  | Bool
otherwise = Approx
Bottom

logInternal :: Approx -> Approx
logInternal :: Approx -> Approx
logInternal Approx
Bottom = String -> Approx
forall a. HasCallStack => String -> a
error String
"LogInternal: impossible"
logInternal (Approx Int
mb Integer
m Integer
e Int
s) =
  let t' :: Int
t' = (Int -> Int
forall a. Num a => a -> a
negate Int
mb) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
10 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
0 (Integer -> Int
integerLog2 Integer
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
s) -- (5 + size of argument) guard digits
      r :: Int
r = Int
s Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Integer -> Int
integerLog2 (Integer
3Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
m) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1
      x :: Dyadic
x = Dyadic -> Int -> Dyadic
forall a. Scalable a => a -> Int -> a
scale (Integer
m Integer -> Int -> Dyadic
:^ Int
s) (-Int
r) -- 2/3 <= x' <= 4/3
      y :: Dyadic
y = Int -> Dyadic -> Dyadic -> Dyadic
divD' Int
t' (Dyadic
x Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
- Dyadic
1) (Dyadic
x Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
+ Dyadic
1) -- so |y| <= 1/5
      (Integer
n :^ Int
s') = (Dyadic -> Int -> Dyadic) -> Int -> Dyadic -> Dyadic
forall a b c. (a -> b -> c) -> b -> a -> c
flip Dyadic -> Int -> Dyadic
forall a. Scalable a => a -> Int -> a
scale Int
1 (Dyadic -> Dyadic) -> Dyadic -> Dyadic
forall a b. (a -> b) -> a -> b
$ Int -> Dyadic -> Dyadic
atanhD Int
t' Dyadic
y
      (Integer
e' :^ Int
s'') = Int -> Dyadic -> Dyadic -> Dyadic
divD' Int
t' (Integer
eInteger -> Int -> Dyadic
:^(Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
r)) Dyadic
x -- Estimate error term.
      res :: Approx
res = Int -> Integer -> Integer -> Int -> Approx
approxMB Int
mb Integer
n (Integer -> Int -> Integer
forall a. Scalable a => a -> Int -> a
scale (Integer
e' Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1) (Int
s'' Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
s')) Int
s'
  in Approx -> Approx
boundErrorTerm (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Approx
res Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
+ Int -> Approx
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
r Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Int -> Approx
log2A Int
t'

-- | Logarithm by binary splitting summation of Taylor series.
logBinarySplittingA :: Precision -> Approx -> Approx
logBinarySplittingA :: Int -> Approx -> Approx
logBinarySplittingA Int
_ Approx
Bottom = Approx
Bottom
logBinarySplittingA Int
res a :: Approx
a@(Approx Int
mb Integer
m Integer
e Int
s) =
    if Integer
m Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
e then Approx
Bottom -- only defined for strictly positive arguments
    else
        let r :: Int
r = Int
s Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Integer -> Int
integerLog2 (Integer
3Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
m) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1
            a' :: Approx
a' = Int -> Integer -> Integer -> Int -> Approx
Approx (Int
mbInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
res) Integer
m Integer
e (Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
r)  -- a' is a scaled by a power of 2 so that 2/3 <= |a'| <= 4/3
            u :: Approx
u = Approx
a' Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
- Approx
1
            v :: Approx
v = Approx
a' Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
+ Approx
1
            u2 :: Approx
u2 = Approx -> Approx
sqrA Approx
u
            v2 :: Approx
v2 = Approx -> Approx
sqrA Approx
v
            Finite Int
res' = Extended Int -> Extended Int -> Extended Int
forall a. Ord a => a -> a -> a
min (Approx -> Extended Int
significance Approx
a) (Int -> Extended Int
forall a. a -> Extended a
Finite Int
res)
            n :: Int
n = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
ceiling (Double -> Int) -> (Double -> Double) -> Double -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2) (Double -> Int) -> Double -> Int
forall a b. (a -> b) -> a -> b
$ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (-Int
res')Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double -> Double
forall a. Floating a => a -> a
log Double
0.2Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double -> Double
forall a. Floating a => a -> a
log Double
2) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1
            (Approx
_, Approx
q, Integer
b, Approx
t) = [Integer]
-> [Integer]
-> [Approx]
-> [Approx]
-> Int
-> Int
-> (Approx, Approx, Integer, Approx)
forall a.
Num a =>
[Integer]
-> [Integer] -> [a] -> [a] -> Int -> Int -> (a, a, Integer, a)
abpq (Integer -> [Integer]
forall a. a -> [a]
repeat Integer
2)
                                [Integer
1,Integer
3..]
                                (Approx
uApprox -> [Approx] -> [Approx]
forall a. a -> [a] -> [a]
:Approx -> [Approx]
forall a. a -> [a]
repeat Approx
u2)
                                (Approx
vApprox -> [Approx] -> [Approx]
forall a. a -> [a] -> [a]
:Approx -> [Approx]
forall a. a -> [a]
repeat Approx
v2)
                                Int
0
                                Int
n
            nextTerm :: Approx
nextTerm = Approx -> Approx
recipA (Int -> Approx -> Approx
setMB (Int
mbInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
res) Approx
5) Approx -> Int -> Approx
forall a b. (Fractional a, Integral b) => a -> b -> a
^^ (Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)
        in Approx -> Approx
boundErrorTerm (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Approx -> Approx -> Approx
fudge (Approx
t Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Approx -> Approx
recipA (Integer -> Approx
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
bApprox -> Approx -> Approx
forall a. Num a => a -> a -> a
*Approx
q) Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
+ Int -> Approx
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
r Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Int -> Approx
log2A (-Int
res)) Approx
nextTerm

-- | Logarithm by summation of Taylor series.
logTaylorA :: Precision -> Approx -> Approx
logTaylorA :: Int -> Approx -> Approx
logTaylorA Int
_ Approx
Bottom = Approx
Bottom
logTaylorA Int
res (Approx Int
mb Integer
m Integer
e Int
s) =
    if Integer
m Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
e then Approx
Bottom -- only defined for strictly positive arguments
    else
        let res' :: Int
res' = Int
res Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
errorBits
            r :: Int
r = Int
s Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Integer -> Int
integerLog2 (Integer
3Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
m) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1
            a' :: Approx
a' = Int -> Integer -> Integer -> Int -> Approx
approxMB Int
mb Integer
m Integer
e (Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
r)  -- a' is a scaled by a power of 2 so that 2/3 <= a' <= 4/3
            u :: Approx
u = Approx
a' Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
- Approx
1
            v :: Approx
v = Approx
a' Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
+ Approx
1
            x :: Approx
x = Approx
u Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Approx -> Approx
recipA Approx
v  -- so |u/v| <= 1/5
            x2 :: Approx
x2 = Approx -> Approx
boundErrorTerm (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Approx -> Approx
sqrA Approx
x
            t :: Approx
t = Int -> [Approx] -> [Integer] -> Approx
taylor
                  Int
res'
                  ((Approx -> Approx) -> Approx -> [Approx]
forall a. (a -> a) -> a -> [a]
iterate (Approx
x2Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
*) Approx
x)
                  [Integer
1,Integer
3..]
        in Approx -> Approx
boundErrorTerm (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Approx
2 Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Approx
t Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
+ Int -> Approx
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
r Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Int -> Approx
log2A (-Int
res')

-- Sine computed by Taylor expansion after 2 stage range reduction.

-- | Computes sine by summation of Taylor series after two levels of range reductions.
sinTaylorA :: Approx -> Approx
sinTaylorA :: Approx -> Approx
sinTaylorA Approx
Bottom = Int -> Integer -> Integer -> Int -> Approx
Approx Int
64 Integer
0 Integer
1 Int
0 -- [-1,1]
sinTaylorA a :: Approx
a@(Approx Int
mb Integer
_ Integer
e Int
_) 
    | Integer
e Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 = Approx -> Approx
sinTaylorRed2A Approx
aRed
    | Bool
otherwise = Approx
sL Approx -> Approx -> Approx
`unionA` Approx
sR -- aRed is in the interval [-π/2,π/2] where sine is monotone
    where
    (Approx
aRed, (Maybe Approx
maRedL, Maybe Approx
maRedR)) = Approx -> (Approx, (Maybe Approx, Maybe Approx))
sinTaylorRed1A Approx
a
    sL :: Approx
sL =
        case Maybe Approx
maRedL of
            Maybe Approx
Nothing -> Int -> Integer -> Integer -> Int -> Approx
Approx Int
mb (-Integer
1) Integer
0 Int
0 -- aRed probably contains -pi/2
            Just Approx
aRedL -> Approx -> Approx
sinTaylorRed2A Approx
aRedL
    sR :: Approx
sR =
        case Maybe Approx
maRedR of
            Maybe Approx
Nothing -> Int -> Integer -> Integer -> Int -> Approx
Approx Int
mb Integer
1 Integer
0 Int
0 -- aRed probably contains +pi/2
            Just Approx
aRedR -> Approx -> Approx
sinTaylorRed2A Approx
aRedR

-- | First level of range reduction for sine. Brings it into the interval [-π/2,π/2].
sinTaylorRed1A :: Approx -> (Approx, (Maybe Approx, Maybe Approx))
sinTaylorRed1A :: Approx -> (Approx, (Maybe Approx, Maybe Approx))
sinTaylorRed1A Approx
Bottom = (Approx
Bottom, (Maybe Approx
forall a. Maybe a
Nothing, Maybe Approx
forall a. Maybe a
Nothing))
sinTaylorRed1A a :: Approx
a@(Approx Int
mb Integer
_ Integer
_ Int
_) = 
  let _pi :: Approx
_pi = Int -> Approx
piA (Int
mbInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
10)
      _halfPi :: Approx
_halfPi = Approx -> Int -> Approx
forall a. Scalable a => a -> Int -> a
scale Approx
_pi (-Int
1)
      x :: Approx
x = Int -> Approx -> Approx
setMB Int
mb (Approx -> Approx) -> (Approx -> Approx) -> Approx -> Approx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
subtract Approx
_halfPi) (Approx -> Approx) -> (Approx -> Approx) -> Approx -> Approx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Approx -> Approx
forall a. Num a => a -> a
abs (Approx -> Approx) -> (Approx -> Approx) -> Approx -> Approx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Approx
_pi Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
-) (Approx -> Approx) -> (Approx -> Approx) -> Approx -> Approx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Approx -> Approx
forall a. Num a => a -> a
abs (Approx -> Approx) -> (Approx -> Approx) -> Approx -> Approx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
subtract Approx
_halfPi) (Approx -> Approx) -> (Approx -> Approx) -> Approx -> Approx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Approx -> Approx -> Approx
modA Approx
a (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Approx
2Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
*Approx
_pi
      xL :: Approx
xL = Approx -> Approx
lowerA Approx
x
      xR :: Approx
xR = Approx -> Approx
upperA Approx
x
      _halfPiL :: Approx
_halfPiL = Approx -> Approx
lowerA Approx
_halfPi
  in (Approx
x, 
        (if (- Approx
_halfPiL) Approx -> Approx -> Bool
forall a. Ord a => a -> a -> Bool
<= Approx
xL       then Approx -> Maybe Approx
forall a. a -> Maybe a
Just Approx
xL else Maybe Approx
forall a. Maybe a
Nothing, -- guarantee -π/2 <= xL
         if           Approx
xR Approx -> Approx -> Bool
forall a. Ord a => a -> a -> Bool
<= Approx
_halfPiL then Approx -> Maybe Approx
forall a. a -> Maybe a
Just Approx
xR else Maybe Approx
forall a. Maybe a
Nothing)) -- guarantee xR <= π/2
  
-- | Second level of range reduction for sine.
sinTaylorRed2A :: Approx -> Approx
sinTaylorRed2A :: Approx -> Approx
sinTaylorRed2A Approx
Bottom = Int -> Integer -> Integer -> Int -> Approx
Approx Int
64 Integer
0 Integer
1 Int
0 -- [-1,1]
sinTaylorRed2A a :: Approx
a@(Approx Int
mb Integer
m Integer
_ Int
s) = 
  let k :: Int
k = Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
0 (Integer -> Int
integerLog2 Integer
m Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
s Int -> Int -> Int
forall a. Num a => a -> a -> a
+ (Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double -> Int) -> (Double -> Double) -> Double -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Int) -> Double -> Int
forall a b. (a -> b) -> a -> b
$ Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
mb))
      a' :: Approx
a' = Approx
a Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Approx -> Approx
recipA ((Int -> Approx -> Approx
setMB Int
mb Approx
3)Approx -> Int -> Approx
forall a b. (Num a, Integral b) => a -> b -> a
^Int
k)
      a2 :: Approx
a2 = Approx -> Approx
forall a. Num a => a -> a
negate (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Approx -> Approx
sqrA Approx
a'
      t :: Approx
t = Int -> [Approx] -> Approx -> Approx
taylorA Int
mb ((Approx -> Approx) -> [Approx] -> [Approx]
forall a b. (a -> b) -> [a] -> [b]
map (Approx -> Approx
recipA (Approx -> Approx) -> (Approx -> Approx) -> Approx -> Approx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Approx -> Approx
setMB Int
mb) [Approx]
forall a. Num a => [a]
oddFac) Approx
a2
      step :: Approx -> Approx
step Approx
x = Approx -> Approx
boundErrorTerm (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Approx
x Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* (Approx
3 Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
- Approx
4 Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Approx -> Approx
sqrA Approx
x)
  in ([Approx] -> Int -> Approx
forall a. [a] -> Int -> a
!! Int
k) ([Approx] -> Approx) -> (Approx -> [Approx]) -> Approx -> Approx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Approx -> Approx) -> Approx -> [Approx]
forall a. (a -> a) -> a -> [a]
iterate Approx -> Approx
step (Approx -> [Approx]) -> (Approx -> Approx) -> Approx -> [Approx]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Approx -> Approx
boundErrorTerm (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Approx
t Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Approx
a'

-- | Computes the sine of an approximation. Chooses the best implementation.
sinA :: Approx -> Approx
sinA :: Approx -> Approx
sinA = Approx -> Approx
sinTaylorA

-- | Computes the cosine of an approximation. Chooses the best implementation.
cosA :: Approx -> Approx
cosA :: Approx -> Approx
cosA Approx
Bottom = Int -> Integer -> Integer -> Int -> Approx
Approx Int
64 Integer
0 Integer
1 Int
0 -- [-1,1]
cosA x :: Approx
x@(Approx Int
mb Integer
_ Integer
_ Int
_) = Approx -> Approx
sinA ((Int -> Integer -> Integer -> Int -> Approx
Approx Int
1 Integer
1 Integer
0 (-Int
1)) Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Int -> Approx
piA (Int
mbInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
2) Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
- Approx
x)

-- | Computes the sine of an approximation by binary splitting summation of Taylor series.
--
-- Begins by reducing the interval to [0,π/4].
sinBinarySplittingA :: Precision -> Approx -> Approx
sinBinarySplittingA :: Int -> Approx -> Approx
sinBinarySplittingA Int
_ Approx
Bottom = Approx
Bottom
sinBinarySplittingA Int
res Approx
a =
    let _pi :: Approx
_pi = Int -> Approx
piBorweinA Int
res
        (Approx Int
mb' Integer
m' Integer
e' Int
s') = Approx
4 Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Approx
a Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Approx -> Approx
recipA Approx
_pi
        (Integer
k,Integer
m1) = Integer
m' Integer -> Integer -> (Integer, Integer)
forall a. Integral a => a -> a -> (a, a)
`divMod` Int -> Integer
forall a. Bits a => Int -> a
bit (-Int
s')
        a2 :: Approx
a2 = Approx
_pi Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Int -> Dyadic -> Approx
fromDyadicMB Int
mb' (Integer
1Integer -> Int -> Dyadic
:^(-Int
2)) Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* (Int -> Integer -> Integer -> Int -> Approx
Approx Int
mb' Integer
m1 Integer
e' Int
s')
    in case Integer
k Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
8 of
         Integer
0 -> Int -> Approx -> Approx
sinInRangeA Int
res Approx
a2
         Integer
1 -> Int -> Approx -> Approx
cosInRangeA Int
res (Approx
_pi Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Int -> Dyadic -> Approx
fromDyadicMB Int
mb' (Integer
1Integer -> Int -> Dyadic
:^(-Int
2)) Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
- Approx
a2)
         Integer
2 -> Int -> Approx -> Approx
cosInRangeA Int
res Approx
a2
         Integer
3 -> Int -> Approx -> Approx
sinInRangeA Int
res (Approx
_pi Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Int -> Dyadic -> Approx
fromDyadicMB Int
mb' (Integer
1Integer -> Int -> Dyadic
:^(-Int
2)) Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
- Approx
a2)
         Integer
4 -> - Int -> Approx -> Approx
sinInRangeA Int
res Approx
a2
         Integer
5 -> - Int -> Approx -> Approx
cosInRangeA Int
res (Approx
_pi Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Int -> Dyadic -> Approx
fromDyadicMB Int
mb' (Integer
1Integer -> Int -> Dyadic
:^(-Int
2)) Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
- Approx
a2)
         Integer
6 -> - Int -> Approx -> Approx
cosInRangeA Int
res Approx
a2
         Integer
7 -> - Int -> Approx -> Approx
sinInRangeA Int
res (Approx
_pi Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Int -> Dyadic -> Approx
fromDyadicMB Int
mb' (Integer
1Integer -> Int -> Dyadic
:^(-Int
2)) Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
- Approx
a2)
         Integer
_ -> String -> Approx
forall a. HasCallStack => String -> a
error String
"Impossible"

-- | Computes the cosine of an approximation by binary splitting summation of Taylor series.
--
-- Begins by reducing the interval to [0,π/4].
cosBinarySplittingA :: Precision -> Approx -> Approx
cosBinarySplittingA :: Int -> Approx -> Approx
cosBinarySplittingA Int
_ Approx
Bottom = Approx
Bottom
cosBinarySplittingA Int
res Approx
a =
    let _pi :: Approx
_pi = Int -> Approx
piBorweinA Int
res
        (Approx Int
mb' Integer
m' Integer
e' Int
s') = Approx
4 Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Approx
a Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Approx -> Approx
recipA Approx
_pi
        (Integer
k,Integer
m1) = Integer
m' Integer -> Integer -> (Integer, Integer)
forall a. Integral a => a -> a -> (a, a)
`divMod` Int -> Integer
forall a. Bits a => Int -> a
bit (-Int
s')
        a2 :: Approx
a2 = Approx
_pi Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Int -> Dyadic -> Approx
fromDyadicMB Int
mb' (Integer
1Integer -> Int -> Dyadic
:^(-Int
2)) Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* (Int -> Integer -> Integer -> Int -> Approx
Approx Int
mb' Integer
m1 Integer
e' Int
s')
    in case Integer
k Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` Integer
8 of
         Integer
0 -> Int -> Approx -> Approx
cosInRangeA Int
res Approx
a2
         Integer
1 -> Int -> Approx -> Approx
sinInRangeA Int
res (Approx
_pi Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Int -> Dyadic -> Approx
fromDyadicMB Int
mb' (Integer
1Integer -> Int -> Dyadic
:^(-Int
2)) Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
- Approx
a2)
         Integer
2 -> - Int -> Approx -> Approx
sinInRangeA Int
res Approx
a2
         Integer
3 -> - Int -> Approx -> Approx
cosInRangeA Int
res (Approx
_pi Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Int -> Dyadic -> Approx
fromDyadicMB Int
mb' (Integer
1Integer -> Int -> Dyadic
:^(-Int
2)) Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
- Approx
a2)
         Integer
4 -> - Int -> Approx -> Approx
cosInRangeA Int
res Approx
a2
         Integer
5 -> - Int -> Approx -> Approx
sinInRangeA Int
res (Approx
_pi Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Int -> Dyadic -> Approx
fromDyadicMB Int
mb' (Integer
1Integer -> Int -> Dyadic
:^(-Int
2)) Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
- Approx
a2)
         Integer
6 -> Int -> Approx -> Approx
sinInRangeA Int
res Approx
a2
         Integer
7 -> Int -> Approx -> Approx
cosInRangeA Int
res (Approx
_pi Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Int -> Dyadic -> Approx
fromDyadicMB Int
mb' (Integer
1Integer -> Int -> Dyadic
:^(-Int
2)) Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
- Approx
a2)
         Integer
_ -> String -> Approx
forall a. HasCallStack => String -> a
error String
"Impossible"


-- | Computes the arc tangent of an approximation. Chooses the best implementation.
atanA :: Precision -> Approx -> Approx
atanA :: Int -> Approx -> Approx
atanA = Int -> Approx -> Approx
atanTaylorA

-- | Computes the arc tangent of an approximation by binary splitting summation of Taylor series.
atanBinarySplittingA :: Precision -> Approx -> Approx
atanBinarySplittingA :: Int -> Approx -> Approx
atanBinarySplittingA Int
_ Approx
Bottom = Approx
Bottom
atanBinarySplittingA Int
res Approx
a =
  let rr :: Approx -> Approx
rr Approx
x = Approx
x Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Approx -> Approx
recipA (Approx
1 Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
+ Approx -> Approx
sqrtA (Approx
1 Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
+ Approx -> Approx
sqrA Approx
x))
      a' :: Approx
a' = Approx -> Approx
rr (Approx -> Approx) -> (Approx -> Approx) -> Approx -> Approx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Approx -> Approx
rr (Approx -> Approx) -> (Approx -> Approx) -> Approx -> Approx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Approx -> Approx
rr (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Approx
a -- range reduction so that |a'| < 1/4
      a2 :: Approx
a2 = - Approx -> Approx
sqrA Approx
a'
      res' :: Int
res' = case (Approx -> Extended Int
significance Approx
a) of
               (Finite Int
_r) -> Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
res Int
_r
               Extended Int
_ -> Int
res
--      Finite res' = min (significance a) (Finite res)
      n :: Int
n = (Int
res' Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2
      (Approx
_, Approx
q, Integer
b, Approx
t) = [Integer]
-> [Integer]
-> [Approx]
-> [Approx]
-> Int
-> Int
-> (Approx, Approx, Integer, Approx)
forall a.
Num a =>
[Integer]
-> [Integer] -> [a] -> [a] -> Int -> Int -> (a, a, Integer, a)
abpq [Integer]
forall a. Num a => [a]
ones
                          [Integer
1,Integer
3..]
                          (Approx
a'Approx -> [Approx] -> [Approx]
forall a. a -> [a] -> [a]
:Approx -> [Approx]
forall a. a -> [a]
repeat Approx
a2)
                          (Approx -> [Approx]
forall a. a -> [a]
repeat Approx
1)
                          Int
0
                          Int
n
      nextTerm :: Approx
nextTerm = Int -> Integer -> Integer -> Int -> Approx
Approx (Approx -> Int
mBound Approx
a Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
res) Integer
1 Integer
0 (-Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
n)
  in Approx -> Approx
boundErrorTerm (Approx -> Approx) -> (Approx -> Approx) -> Approx -> Approx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Approx
8Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
*) (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Approx -> Approx -> Approx
fudge (Approx
t Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Approx -> Approx
recipA (Integer -> Approx
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
bApprox -> Approx -> Approx
forall a. Num a => a -> a -> a
*Approx
q)) Approx
nextTerm

-- + Bottom
-- + Deal with sign -- Because of next line, not worthwhile
-- + if lowerbound(abs a) > 2 then pi/2 - atan (1/a) -- Don't want to do this, what if 0 \in a?
-- + else
--   - r = min res (significance a)
--   - k = floor (sqrt r) / 2 `min` 2 (guarantee |x| < 1/2)
--   - x = rr^k(a)
--   - taylor (r + k + 5) (-x^2) [1,3..]
--   - (x*)
--   - same error as x
--   - (2^k *)

atanTaylorA :: Precision -> Approx -> Approx
atanTaylorA :: Int -> Approx -> Approx
atanTaylorA Int
_ Approx
Bottom = Approx
Bottom
atanTaylorA Int
res a :: Approx
a@(Approx Int
mb Integer
_ Integer
_ Int
_) =
  let (Finite Int
r) = Extended Int -> Extended Int -> Extended Int
forall a. Ord a => a -> a -> a
min (Int -> Extended Int
forall (f :: * -> *) a. Applicative f => a -> f a
pure Int
res) (Approx -> Extended Int
significance Approx
a)
      k :: Int
k = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min (Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
floor (Double -> Double
forall a. Floating a => a -> a
sqrt (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
r)) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2) Int
2
      res' :: Int
res' = (Int
mb Int -> Int -> Int
forall a. Ord a => a -> a -> a
`max` Int
res) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
k Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
5
      rr :: Approx -> Approx
rr Approx
_x = Approx
_x Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Approx -> Approx
recipA (Approx
1 Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
+ Approx -> Approx
sqrtA (Approx
1 Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
+ Approx -> Approx
sqrA Approx
_x))
      x :: Approx
x = Approx -> Approx
boundErrorTerm (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ (Approx -> Approx) -> Approx -> [Approx]
forall a. (a -> a) -> a -> [a]
iterate Approx -> Approx
rr (Int -> Approx -> Approx
setMB Int
res' Approx
a) [Approx] -> Int -> Approx
forall a. [a] -> Int -> a
!! Int
k
      x2 :: Approx
x2 = Approx -> Approx
forall a. Num a => a -> a
negate (Approx -> Approx
sqrA Approx
x)
      t :: Approx
t = Approx -> Approx
boundErrorTerm (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Approx
x Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Int -> [Approx] -> Approx -> Approx
taylorA Int
res' ((Approx -> Approx) -> [Approx] -> [Approx]
forall a b. (a -> b) -> [a] -> [b]
map (Approx -> Approx
recipA (Approx -> Approx) -> (Approx -> Approx) -> Approx -> Approx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Approx -> Approx
setMB Int
res') [Approx
1,Approx
3..]) Approx
x2
  in Approx -> Int -> Approx
forall a. Scalable a => a -> Int -> a
scale Approx
t Int
k

-- > let x = fromDouble (-0.2939788524332769)
-- > require 10 $ x
-- Approx (-5295852201093248) 1 (-54)
-- > require 10 . tan $ atan x
-- Approx (-10845905307838971904) 907 (-65)
-- > scale (-5295852201093248) 11
-- -10845905307838971904
--
-- problemet är att 1 måste skalas till 2^11, men blev bara 907
--
-- Men problemet verkar vara i tan, inte i atan.


-- | Computes the arc tangent of an approximation by summation of Taylor series.
-- atanTaylorA :: Precision -> Approx -> Approx
-- atanTaylorA _ Bottom = Bottom
-- atanTaylorA res a =
--   let rr x = x * recipA res (1 + sqrtA res (1 + sqrA x))
--       a' = rr . rr . rr $ a -- range reduction so that |a'| < 1/4
--       a2 = - sqrA a'
--       res' = case (significance a) of
--                (Finite _r) -> min res _r
--                _ -> res
-- --      Finite res' = min (significance a) (Finite res)
--       t = taylorA res' (map (recipA res') [1,3..]) a2
--   in boundErrorTerm . (8*) $ t

{-
swapSinCos :: Precision -> Approx -> Approx
swapSinCos res a = sqrtA res $ 1 - sqrA a
-}

-- Computes sine if second argument is in the range [0,pi/4]
sinInRangeA :: Precision -> Approx -> Approx
sinInRangeA :: Int -> Approx -> Approx
sinInRangeA Int
_ Approx
Bottom = Approx
Bottom
sinInRangeA Int
res Approx
a =
    let n :: Int
n = Int
res Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2        -- need to improve this estimate (is valid from res>=80)
        (Approx
_, Approx
q, Integer
b, Approx
t) = [Integer]
-> [Integer]
-> [Approx]
-> [Approx]
-> Int
-> Int
-> (Approx, Approx, Integer, Approx)
forall a.
Num a =>
[Integer]
-> [Integer] -> [a] -> [a] -> Int -> Int -> (a, a, Integer, a)
abpq [Integer]
forall a. Num a => [a]
ones
                            [Integer]
forall a. Num a => [a]
ones
                            (Approx
aApprox -> [Approx] -> [Approx]
forall a. a -> [a] -> [a]
:Approx -> [Approx]
forall a. a -> [a]
repeat (- Approx -> Approx
sqrA Approx
a))
                            (Approx
1Approx -> [Approx] -> [Approx]
forall a. a -> [a] -> [a]
:[Approx
2Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
*Approx
iApprox -> Approx -> Approx
forall a. Num a => a -> a -> a
*(Approx
2Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
*Approx
iApprox -> Approx -> Approx
forall a. Num a => a -> a -> a
+Approx
1) | Approx
i <- [Approx
1..]] :: [Approx])
                            Int
0
                            Int
n
        nextTerm :: Approx
nextTerm = Int -> Dyadic -> Approx
fromDyadicMB (Approx -> Int
mBound Approx
a Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
res) (Integer
1Integer -> Int -> Dyadic
:^(-Int
res))
    in Approx -> Approx
boundErrorTerm (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Approx -> Approx -> Approx
fudge (Approx
t Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Approx -> Approx
recipA (Integer -> Approx
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
bApprox -> Approx -> Approx
forall a. Num a => a -> a -> a
*Approx
q)) Approx
nextTerm

-- Computes cosine if second argument is in the range [0,pi/4]
cosInRangeA :: Precision -> Approx -> Approx
cosInRangeA :: Int -> Approx -> Approx
cosInRangeA Int
_ Approx
Bottom = Approx
Bottom
cosInRangeA Int
res Approx
a =
    let n :: Int
n = Int
res Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2        -- need to improve this estimate (is valid from res>=80)
        (Approx
_, Approx
q, Integer
b, Approx
t) = [Integer]
-> [Integer]
-> [Approx]
-> [Approx]
-> Int
-> Int
-> (Approx, Approx, Integer, Approx)
forall a.
Num a =>
[Integer]
-> [Integer] -> [a] -> [a] -> Int -> Int -> (a, a, Integer, a)
abpq [Integer]
forall a. Num a => [a]
ones
                            [Integer]
forall a. Num a => [a]
ones
                            (Approx
1Approx -> [Approx] -> [Approx]
forall a. a -> [a] -> [a]
:Approx -> [Approx]
forall a. a -> [a]
repeat (- Approx -> Approx
sqrA Approx
a))
                            (Approx
1Approx -> [Approx] -> [Approx]
forall a. a -> [a] -> [a]
:[Approx
2Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
*Approx
iApprox -> Approx -> Approx
forall a. Num a => a -> a -> a
*(Approx
2Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
*Approx
iApprox -> Approx -> Approx
forall a. Num a => a -> a -> a
-Approx
1) | Approx
i <- [Approx
1..]] :: [Approx])
                            Int
0
                            Int
n
        nextTerm :: Approx
nextTerm = Int -> Dyadic -> Approx
fromDyadicMB (Approx -> Int
mBound Approx
a Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
res) (Integer
1Integer -> Int -> Dyadic
:^(-Int
res))
    in Approx -> Approx
boundErrorTerm (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Approx -> Approx -> Approx
fudge (Approx
t Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Approx -> Approx
recipA (Integer -> Approx
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
bApprox -> Approx -> Approx
forall a. Num a => a -> a -> a
*Approx
q)) Approx
nextTerm

{-|
Computes a sequence of approximations of π using binary splitting summation of
Ramanujan's series. See Haible and Papanikolaou 1998.
-}
piRaw :: [Approx]
piRaw :: [Approx]
piRaw = ((Int, (Integer, Integer, Integer, Integer))
 -> Maybe (Approx, (Int, (Integer, Integer, Integer, Integer))))
-> (Int, (Integer, Integer, Integer, Integer)) -> [Approx]
forall b a. (b -> Maybe (a, b)) -> b -> [a]
unfoldr (Int, (Integer, Integer, Integer, Integer))
-> Maybe (Approx, (Int, (Integer, Integer, Integer, Integer)))
f (Int
1, (Integer
1, Integer
1, Integer
1, Integer
13591409))
    where as :: [Integer]
as = [Integer
13591409,Integer
13591409Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
545140134..]
          bs :: [Integer]
bs = [Integer]
forall a. Num a => [a]
ones
          ps :: [Integer]
ps = (Integer
1Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:[-(Integer
6Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
5)Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*(Integer
2Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
1)Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*(Integer
6Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
1) | Integer
n <- [Integer
1,Integer
2..]])
          qs :: [Integer]
qs = (Integer
1Integer -> [Integer] -> [Integer]
forall a. a -> [a] -> [a]
:[Integer
nInteger -> Integer -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
3Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
640320Integer -> Integer -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
26680 | Integer
n <- [Integer
1,Integer
2..]])
          f :: (Int, (Integer, Integer, Integer, Integer))
-> Maybe (Approx, (Int, (Integer, Integer, Integer, Integer)))
f (Int
i, (Integer
pl, Integer
ql, Integer
bl, Integer
tl)) = 
            let i2 :: Int
i2 = Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
*Int
2
                (Integer
pr, Integer
qr, Integer
br, Integer
tr) = [Integer]
-> [Integer]
-> [Integer]
-> [Integer]
-> Int
-> Int
-> (Integer, Integer, Integer, Integer)
forall a.
Num a =>
[Integer]
-> [Integer] -> [a] -> [a] -> Int -> Int -> (a, a, Integer, a)
abpq [Integer]
as [Integer]
bs [Integer]
ps [Integer]
qs Int
i Int
i2
                n :: Int
n = Int
21Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
47Int -> Int -> Int
forall a. Num a => a -> a -> a
*(Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)
                x :: Approx
x = Integer -> Approx
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
tl Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Approx -> Approx
recipA (Int -> Approx -> Approx
setMB Int
n (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Integer -> Approx
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Integer
blInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
ql))
                x1 :: Approx
x1 = Approx -> Approx -> Approx
fudge Approx
x (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Int -> Dyadic -> Approx
fromDyadicMB Int
n (Integer
1Integer -> Int -> Dyadic
:^(-Int
n))
                x2 :: Approx
x2 = Approx -> Approx
boundErrorTermMB (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ Approx -> Approx
sqrtA (Int -> Approx -> Approx
setMB Int
n Approx
1823176476672000) Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
* Approx -> Approx
recipA Approx
x1
            in (Approx, (Int, (Integer, Integer, Integer, Integer)))
-> Maybe (Approx, (Int, (Integer, Integer, Integer, Integer)))
forall a. a -> Maybe a
Just ( Approx
x2
                    , (Int
i2, (Integer
pl Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
pr, Integer
ql Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
qr, Integer
bl Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
br, Integer -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
br Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
qr Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
tl Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
bl Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
pl Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
tr))
                    )

-- | Computes an 'Approx' of π of the given precision.
piA :: Precision -> Approx
piA :: Int -> Approx
piA Int
res = Int -> Approx -> Approx
limitAndBound Int
res (Approx -> Approx) -> ([Approx] -> Approx) -> [Approx] -> Approx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Approx] -> Approx
forall a. [a] -> a
head ([Approx] -> Approx) -> [Approx] -> Approx
forall a b. (a -> b) -> a -> b
$ (Approx -> Bool) -> [Approx] -> [Approx]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile ((Extended Int -> Extended Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int -> Extended Int
forall (f :: * -> *) a. Applicative f => a -> f a
pure Int
res) (Extended Int -> Bool)
-> (Approx -> Extended Int) -> Approx -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Approx -> Extended Int
precision) [Approx]
piRaw

{-|
Second argument is noise to be added to first argument. Used to allow for the
error term when truncating a series.
-}
fudge :: Approx -> Approx -> Approx
fudge :: Approx -> Approx -> Approx
fudge Approx
a (Approx Int
_ Integer
0 Integer
0 Int
_) = Approx
a
fudge (Approx Int
mb Integer
m Integer
0 Int
s) (Approx Int
mb' Integer
m' Integer
e' Int
s') =
  Int -> Int -> Integer -> Integer -> Int -> Approx
approxMB2 Int
mb Int
mb' (Integer
m Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
`shift` (Int
s Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
s')) (Integer -> Integer
forall a. Num a => a -> a
abs Integer
m' Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
e' Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1) Int
s'
fudge (Approx Int
mb Integer
m Integer
e Int
s) (Approx Int
mb' Integer
m' Integer
e' Int
s') =
  let m'' :: Integer
m'' = Integer
1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ (Integer -> Integer
forall a. Num a => a -> a
abs Integer
m' Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
e') Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
`shift` (Int
s' Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
s Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
  in Int -> Int -> Integer -> Integer -> Int -> Approx
approxMB2 Int
mb Int
mb' Integer
m (Integer
eInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
m'') Int
s
fudge Approx
_ Approx
_  = Approx
Bottom

--

-- | Compute π using Machin's formula. Lifted from computation on dyadic numbers.
piMachinA :: Precision -> Approx
piMachinA :: Int -> Approx
piMachinA Int
t = let (Integer
m:^Int
s) = Int -> Dyadic
piMachinD (-Int
t) in Integer -> Integer -> Int -> Approx
approxAutoMB Integer
m Integer
1 Int
s

-- | Compute π using AGM as described in Borwein and Borwein's book 'Pi and
-- the AGM'. Lifted from computation on dyadic numbers.
piBorweinA :: Precision -> Approx
piBorweinA :: Int -> Approx
piBorweinA Int
t = let (Integer
m:^Int
s) = Int -> Dyadic
piBorweinD (-Int
t) in Integer -> Integer -> Int -> Approx
approxAutoMB Integer
m Integer
1 Int
s


-- | Compute π using AGM as described in Borwein and Borwein's book 'Pi and
-- the AGM'.
piAgmA :: Precision -> Approx -> Approx
piAgmA :: Int -> Approx -> Approx
piAgmA Int
t x_ :: Approx
x_@(Approx Int
mb_ Integer
_ Integer
_ Int
_) = 
             let -- t' = t - 10
                 mb :: Int
mb = Int
mb_ Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
t
                 a :: Approx
a = Int -> Approx -> Approx
setMB Int
mb Approx
1
                 x :: Approx
x = Int -> Approx -> Approx
setMB Int
mb Approx
x_
                 b :: Approx
b = Approx -> Approx
boundErrorTerm (Approx -> Approx) -> Approx -> Approx
forall a b. (a -> b) -> a -> b
$ (Approx
2Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
*Approx
xApprox -> Approx -> Approx
forall a. Num a => a -> a -> a
*Approx -> Approx
recipA (Approx
xApprox -> Integer -> Approx
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
-Approx
1))Approx -> Integer -> Approx
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2
                 ss :: [(Approx, Approx)]
ss = Int -> Approx -> Approx -> [(Approx, Approx)]
agmA Int
t Approx
a Approx
b
                 c :: Approx
c = Approx -> Approx
boundErrorTerm (Approx -> Approx)
-> ([(Approx, Approx)] -> Approx) -> [(Approx, Approx)] -> Approx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Approx
1Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
-) (Approx -> Approx)
-> ([(Approx, Approx)] -> Approx) -> [(Approx, Approx)] -> Approx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
*Approx -> Approx
recipA (Approx
1Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
-Approx
bApprox -> Integer -> Approx
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2)) (Approx -> Approx)
-> ([(Approx, Approx)] -> Approx) -> [(Approx, Approx)] -> Approx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Approx] -> Approx
agm2 ([Approx] -> Approx)
-> ([(Approx, Approx)] -> [Approx]) -> [(Approx, Approx)] -> Approx
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [(Approx, Approx)] -> [Approx]
agm1 ([(Approx, Approx)] -> Approx) -> [(Approx, Approx)] -> Approx
forall a b. (a -> b) -> a -> b
$ [(Approx, Approx)]
ss
                 d :: Approx
d = Approx -> Approx
sqrtA (Approx
1Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
+Approx
b)
                 b2 :: Approx
b2 = Approx
bApprox -> Integer -> Approx
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2
                 b3 :: Approx
b3 = Approx
b2Approx -> Approx -> Approx
forall a. Num a => a -> a -> a
*Approx
b
                 b4 :: Approx
b4 = Approx
b2Approx -> Integer -> Approx
forall a b. (Num a, Integral b) => a -> b -> a
^Integer
2
                 l :: Approx
l = Approx -> Approx
boundErrorTerm (Approx -> Approx) -> Approx -> Approx
forall a b