{- |The Dyadic module provides dyadic numbers. Syntax for dyadic number
   'a * 2 ^ s' is 'a :^ s'. The exponent 's' is an 'Int', but the 'a' is an
   arbitrary 'Integer'.
-}
module Data.CDAR.Dyadic (Dyadic(..),shiftD,sqrtD,sqrtD',sqrtRecD,sqrtRecD',initSqrtRecD,initSqrtRecDoubleD,piMachinD,piBorweinD,ln2D,logD,atanhD,divD,divD') where

import Data.Ratio
import Data.Bits
import Data.CDAR.Classes
import Data.CDAR.IntegerLog

-- Dyadic numbers

infix 8 :^

-- |The Dyadic data type.
data Dyadic = Integer :^ Int deriving (ReadPrec [Dyadic]
ReadPrec Dyadic
Int -> ReadS Dyadic
ReadS [Dyadic]
(Int -> ReadS Dyadic)
-> ReadS [Dyadic]
-> ReadPrec Dyadic
-> ReadPrec [Dyadic]
-> Read Dyadic
forall a.
(Int -> ReadS a)
-> ReadS [a] -> ReadPrec a -> ReadPrec [a] -> Read a
readListPrec :: ReadPrec [Dyadic]
$creadListPrec :: ReadPrec [Dyadic]
readPrec :: ReadPrec Dyadic
$creadPrec :: ReadPrec Dyadic
readList :: ReadS [Dyadic]
$creadList :: ReadS [Dyadic]
readsPrec :: Int -> ReadS Dyadic
$creadsPrec :: Int -> ReadS Dyadic
Read,Int -> Dyadic -> ShowS
[Dyadic] -> ShowS
Dyadic -> String
(Int -> Dyadic -> ShowS)
-> (Dyadic -> String) -> ([Dyadic] -> ShowS) -> Show Dyadic
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Dyadic] -> ShowS
$cshowList :: [Dyadic] -> ShowS
show :: Dyadic -> String
$cshow :: Dyadic -> String
showsPrec :: Int -> Dyadic -> ShowS
$cshowsPrec :: Int -> Dyadic -> ShowS
Show)

instance Eq Dyadic where
    Integer
a :^ Int
s == :: Dyadic -> Dyadic -> Bool
== Integer
b :^ Int
t | Int
s Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
t    = Integer
a Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
b (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
s)
                     | Bool
otherwise = Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
a (Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
t) Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
b

instance Ord Dyadic where
    compare :: Dyadic -> Dyadic -> Ordering
compare (Integer
a :^ Int
s) (Integer
b :^ Int
t) | Int
s Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
t    = Integer -> Integer -> Ordering
forall a. Ord a => a -> a -> Ordering
compare Integer
a (Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
b (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
s))
                              | Bool
otherwise = Integer -> Integer -> Ordering
forall a. Ord a => a -> a -> Ordering
compare (Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
a (Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
t)) Integer
b

instance Scalable Dyadic where
  scale :: Dyadic -> Int -> Dyadic
scale (Integer
a :^ Int
s) Int
n = Integer
a Integer -> Int -> Dyadic
:^ (Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
n)

instance Num Dyadic where
    Integer
a :^ Int
s + :: Dyadic -> Dyadic -> Dyadic
+ Integer
b :^ Int
t 
        | Int
s Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
t    = (Integer
a Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer -> Int -> Integer
forall a. Scalable a => a -> Int -> a
scale Integer
b (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
s)) Integer -> Int -> Dyadic
:^ Int
s
        | Bool
otherwise = (Integer -> Int -> Integer
forall a. Scalable a => a -> Int -> a
scale Integer
a (Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
t) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
b) Integer -> Int -> Dyadic
:^ Int
t
    Integer
a :^ Int
s * :: Dyadic -> Dyadic -> Dyadic
* Integer
b :^ Int
t = (Integer
a Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
b) Integer -> Int -> Dyadic
:^ (Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
t)
    negate :: Dyadic -> Dyadic
negate (Integer
a :^ Int
s) = (Integer -> Integer
forall a. Num a => a -> a
negate Integer
a) Integer -> Int -> Dyadic
:^ Int
s
    abs :: Dyadic -> Dyadic
abs (Integer
a :^ Int
s) = (Integer -> Integer
forall a. Num a => a -> a
abs Integer
a) Integer -> Int -> Dyadic
:^ Int
s
    signum :: Dyadic -> Dyadic
signum (Integer
a :^ Int
_) = (Integer -> Integer
forall a. Num a => a -> a
signum Integer
a)Integer -> Int -> Dyadic
:^ Int
0
    fromInteger :: Integer -> Dyadic
fromInteger Integer
i = Integer
i Integer -> Int -> Dyadic
:^ Int
0

instance Real Dyadic where
    toRational :: Dyadic -> Rational
toRational (Integer
a :^ Int
s) = (Integer -> Rational
forall a. Real a => a -> Rational
toRational Integer
a)Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
*Rational
2Rational -> Int -> Rational
forall a b. (Fractional a, Integral b) => a -> b -> a
^^Int
s

-- | Shift a dyadic number to a given base and round in case of right shifts.
shiftD :: Int -> Dyadic -> Dyadic
shiftD :: Int -> Dyadic -> Dyadic
shiftD Int
t (Integer
m:^Int
s) =
    if Int
t Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
s
    then 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
t) Integer -> Int -> Dyadic
:^ Int
t
    else Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftR (Integer
m Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Int -> Integer
forall a. Bits a => Int -> a
bit (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)) (Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
s) Integer -> Int -> Dyadic
:^ Int
t

-- Square root

divD :: Dyadic -> Dyadic -> Dyadic
divD :: Dyadic -> Dyadic -> Dyadic
divD Dyadic
a (Integer
n:^Int
t) = let (Integer
m:^Int
_) = Int -> Dyadic -> Dyadic
shiftD (Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
t) Dyadic
a
                 in Rational -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
round (Integer
m Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% Integer
n) Integer -> Int -> Dyadic
:^ Int
t

{- |Computes the square root of a Dyadic number to the specified base. The
   Newton-Raphson method may overestimates the square root, but the
   overestimate is bounded by 1 ulp. For example, sqrtD 0 2 will give 2,
   whereas the closest integer to the square root is 1. Need double precision
   to guarantee correct rounding, which was not considered worthwhile.

   This is actually Heron's method and is no longer used in Approx as it is
   faster to use sqrtRecD.
-}
sqrtD :: Int -> Dyadic -> Dyadic
sqrtD :: Int -> Dyadic -> Dyadic
sqrtD Int
t Dyadic
x = Int -> Dyadic -> Dyadic -> Dyadic
sqrtD' Int
t Dyadic
x (Dyadic -> Dyadic) -> Dyadic -> Dyadic
forall a b. (a -> b) -> a -> b
$ Dyadic -> Dyadic
initSqrtD Dyadic
x
    where
      -- Initial approximation of square root to about 3 bits
      initSqrtD :: Dyadic -> Dyadic
      initSqrtD :: Dyadic -> Dyadic
initSqrtD (Integer
0:^Int
_) = (Integer
0Integer -> Int -> Dyadic
:^Int
0)
      initSqrtD (Integer
m:^Int
s) = let i :: Int
i = Integer -> Int
integerLog2 Integer
m
                             n :: Integer
n = Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
shift Integer
m (Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
-Int
i)
                             s' :: Int
s' = (Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
i) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
3
                         in if Int -> Bool
forall a. Integral a => a -> Bool
odd (Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
i)
                            then (Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
8)Integer -> Int -> Dyadic
:^Int
s'
                            else (Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
4)Integer -> Int -> Dyadic
:^Int
s'

-- |Square root with initial approximation as third argument.
sqrtD' :: Int -> Dyadic -> Dyadic -> Dyadic
sqrtD' :: Int -> Dyadic -> Dyadic -> Dyadic
sqrtD' Int
t x :: Dyadic
x@(Integer
m:^Int
_) Dyadic
y
    | Integer
m Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 = Integer
0Integer -> Int -> Dyadic
:^Int
0
    | Integer
m Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
0  = [Dyadic] -> Dyadic
converge ([Dyadic] -> Dyadic) -> (Dyadic -> [Dyadic]) -> Dyadic -> Dyadic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Dyadic -> Dyadic) -> Dyadic -> [Dyadic]
forall a. (a -> a) -> a -> [a]
iterate (Dyadic -> Int -> Dyadic -> Dyadic
newtonStep Dyadic
x Int
t) (Dyadic -> Dyadic) -> Dyadic -> Dyadic
forall a b. (a -> b) -> a -> b
$ Int -> Dyadic -> Dyadic
shiftD Int
t Dyadic
y
    | Bool
otherwise = String -> Dyadic
forall a. HasCallStack => String -> a
error String
"Attempting sqrt of negative dyadic number."
    where
      -- One step of Newton iteration to find square root.
      -- The base of a need to be the same as t', and the result will have base t'.
      newtonStep :: Dyadic -> Int -> Dyadic -> Dyadic
      newtonStep :: Dyadic -> Int -> Dyadic -> Dyadic
newtonStep Dyadic
d Int
t' Dyadic
a = Int -> Dyadic -> Dyadic
shiftD Int
t' (Dyadic -> Dyadic) -> Dyadic -> Dyadic
forall a b. (a -> b) -> a -> b
$ (Integer
1Integer -> Int -> Dyadic
:^(-Int
1)) Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
* (Dyadic
a Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
+ Dyadic -> Dyadic -> Dyadic
divD Dyadic
d Dyadic
a)

-- |Reciprocal of square root using Newton's iteration.
sqrtRecD :: Int -> Dyadic -> Dyadic
sqrtRecD :: Int -> Dyadic -> Dyadic
sqrtRecD Int
t Dyadic
a = Int -> Dyadic -> Dyadic -> Dyadic
sqrtRecD' Int
t Dyadic
a (Dyadic -> Dyadic) -> Dyadic -> Dyadic
forall a b. (a -> b) -> a -> b
$ Dyadic -> Dyadic
initSqrtRecDoubleD Dyadic
a

-- |Gives initial values for the iteration. Based on the three most
-- significant bits of the argument. Should consider to use machine floating
-- point to give initial value.
initSqrtRecD :: Dyadic -> Dyadic
initSqrtRecD :: Dyadic -> Dyadic
initSqrtRecD (Integer
0 :^ Int
_) = String -> Dyadic
forall a. HasCallStack => String -> a
error String
"Divide by zero in initSqrtRecD"
initSqrtRecD (Integer
m :^ Int
s) =
  let i :: Int
i = Integer -> Int
integerLog2 Integer
m
      n :: Integer
n = Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
shift Integer
m (Int
3Int -> Int -> Int
forall a. Num a => a -> a -> a
-Int
i)
      s' :: Int
s' = (-Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
5
  in if Int -> Bool
forall a. Integral a => a -> Bool
even (Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
i)
     then ([Integer
62,Integer
59,Integer
56,Integer
53,Integer
51,Integer
49,Integer
48,Integer
46][Integer] -> Int -> Integer
forall a. [a] -> Int -> a
!!(Integer -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
8)) Integer -> Int -> Dyadic
:^ Int
s'
     else ([Integer
44,Integer
42,Integer
40,Integer
38,Integer
36,Integer
35,Integer
34,Integer
33][Integer] -> Int -> Integer
forall a. [a] -> Int -> a
!!(Integer -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
8)) Integer -> Int -> Dyadic
:^ Int
s'

-- Using Double to find a good first approximation to SqrtRec. Extracts 52
-- bits from the dyadic number and encodes a float in the range [1/2,2).
-- Decodes the result of the Double computation as a dyadic taking care of the
-- exponent explicitly (as the range of Double exponents is too small).
initSqrtRecDoubleD :: Dyadic -> Dyadic
initSqrtRecDoubleD :: Dyadic -> Dyadic
initSqrtRecDoubleD (Integer
0 :^ Int
_) = String -> Dyadic
forall a. HasCallStack => String -> a
error String
"Divide by zero in initSqrtRecDoubleD"
initSqrtRecDoubleD (Integer
m :^ Int
s) =
  let i :: Int
i = Integer -> Int
integerLog2 Integer
m
      n :: Integer
n = Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
shift Integer
m (Int
52Int -> Int -> Int
forall a. Num a => a -> a -> a
-Int
i)
      s' :: Int
s' = (-Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
s) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
52
      (Integer
m',Int
_) = Double -> (Integer, Int)
forall a. RealFloat a => a -> (Integer, Int)
decodeFloat (Double -> (Integer, Int))
-> (Double -> Double) -> Double -> (Integer, Int)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/) (Double -> Double) -> (Double -> Double) -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> (Integer, Int)) -> Double -> (Integer, Int)
forall a b. (a -> b) -> a -> b
$ Integer -> Int -> Double
forall a. RealFloat a => Integer -> Int -> a
encodeFloat Integer
n (if Int -> Bool
forall a. Integral a => a -> Bool
even (Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
i) then (-Int
52) else (-Int
51))
      t :: Int
t = if Integer
m' Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= Int -> Integer
forall a. Bits a => Int -> a
bit Int
52 Bool -> Bool -> Bool
&& Int -> Bool
forall a. Integral a => a -> Bool
even (Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
i) then Int
s'Int -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1 else Int
s'
  in Integer
m' Integer -> Int -> Dyadic
:^ Int
t

-- |Reciprocal of square root using Newton's iteration with inital value as third argument.
sqrtRecD' :: Int -> Dyadic -> Dyadic -> Dyadic
sqrtRecD' :: Int -> Dyadic -> Dyadic -> Dyadic
sqrtRecD' Int
t Dyadic
a Dyadic
x0 =
  let step :: Dyadic -> Dyadic
step Dyadic
x = Int -> Dyadic -> Dyadic
shiftD Int
t (Dyadic -> Dyadic) -> Dyadic -> Dyadic
forall a b. (a -> b) -> a -> b
$ Dyadic
x Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
+ Int -> Dyadic -> Dyadic
shiftD Int
t (Dyadic
x Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
* (Dyadic
1 Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
- Int -> Dyadic -> Dyadic
shiftD Int
t (Dyadic
x Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
* Dyadic
x Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
* Dyadic
a)) Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
* (Integer
1 Integer -> Int -> Dyadic
:^ (-Int
1)))
      xs :: [Dyadic]
xs = (Dyadic -> Dyadic) -> Dyadic -> [Dyadic]
forall a. (a -> a) -> a -> [a]
iterate Dyadic -> Dyadic
step Dyadic
x0
  in [Dyadic] -> Dyadic
converge [Dyadic]
xs

-- isqrt :: Integral t => t -> t
-- isqrt 0 = 0
-- isqrt 1 = 1
-- isqrt n = head $ dropWhile (\x -> x*x > n) $ iterate (\x -> (x + n `div` x) `div` 2) (n `div` 2)

-- Assuming a seqence of dyadic numbers with the same exponent converges
-- quadratically. Then this returns the first element where it is known that
-- the sequence will become essentially constant (differing by at most 1). If
-- the numbers contain sufficiently many bits we need only wait until the
-- difference has less than half as many bits (with a small margin).
converge :: [Dyadic] -> Dyadic
converge :: [Dyadic] -> Dyadic
converge ((Integer
0:^Int
s):[Dyadic]
_) = (Integer
0Integer -> Int -> Dyadic
:^Int
s)
converge ((Integer
m:^Int
_):ds :: [Dyadic]
ds@(d :: Dyadic
d@(Integer
n:^Int
_):[Dyadic]
_)) =
  let a :: Int
a = Integer -> Int
integerLog2 Integer
m
      b :: Int
b = 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
n))
  in if Int
a Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
6
     then if Integer -> Integer
forall a. Num a => a -> a
abs (Integer
mInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
n) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
1 then Dyadic
d else [Dyadic] -> Dyadic
converge [Dyadic]
ds
     else if Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
b Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
a Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
3 then Dyadic
d else [Dyadic] -> Dyadic
converge [Dyadic]
ds
converge [Dyadic]
_ = String -> Dyadic
forall a. HasCallStack => String -> a
error String
"List terminating in converge."

divD' :: Int -> Dyadic -> Dyadic -> Dyadic
divD' :: Int -> Dyadic -> Dyadic -> Dyadic
divD' Int
p Dyadic
a Dyadic
b = let (Integer
m:^Int
_) = Int -> Dyadic -> Dyadic
shiftD (Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
p) Dyadic
a
                  (Integer
n:^Int
_) = Int -> Dyadic -> Dyadic
shiftD Int
p Dyadic
b
              in Rational -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
round (Integer
mInteger -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%Integer
n) Integer -> Int -> Dyadic
:^ Int
p

-- |Compute dyadic values close to pi by Machin's formula.
piMachinD :: Int -> Dyadic
piMachinD :: Int -> Dyadic
piMachinD Int
t = let t' :: Int
t' = Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
10Int -> Int -> Int
forall a. Num a => a -> a -> a
-Integer -> Int
integerLog2 (Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Int
forall a. Num a => a -> a
abs Int
t))
                  a :: [Dyadic]
a = (Integer -> Dyadic) -> [Integer] -> [Dyadic]
forall a b. (a -> b) -> [a] -> [b]
map ((Integer -> Int -> Dyadic
:^Int
t') (Integer -> Dyadic) -> (Integer -> Integer) -> Integer -> Dyadic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Rational -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
round (Rational -> Integer)
-> (Integer -> Rational) -> Integer -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int -> Integer
forall a. Bits a => Int -> a
bit (-Int
t') Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%)) ([Integer] -> [Dyadic]) -> [Integer] -> [Dyadic]
forall a b. (a -> b) -> a -> b
$ (Integer -> Integer) -> Integer -> [Integer]
forall a. (a -> a) -> a -> [a]
iterate ((-Integer
25)Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*) (Integer
5 :: Integer)
                  b :: [Dyadic]
b = (Integer -> Dyadic) -> [Integer] -> [Dyadic]
forall a b. (a -> b) -> [a] -> [b]
map ((Integer -> Int -> Dyadic
:^Int
t') (Integer -> Dyadic) -> (Integer -> Integer) -> Integer -> Dyadic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Rational -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
round (Rational -> Integer)
-> (Integer -> Rational) -> Integer -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int -> Integer
forall a. Bits a => Int -> a
bit (-Int
t') Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%)) ([Integer] -> [Dyadic]) -> [Integer] -> [Dyadic]
forall a b. (a -> b) -> a -> b
$ (Integer -> Integer) -> Integer -> [Integer]
forall a. (a -> a) -> a -> [a]
iterate ((-Integer
57121)Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*) (Integer
239 :: Integer)
                  c :: [Dyadic]
c = (Integer -> Dyadic) -> [Integer] -> [Dyadic]
forall a b. (a -> b) -> [a] -> [b]
map ((Integer -> Int -> Dyadic
:^Int
t') (Integer -> Dyadic) -> (Integer -> Integer) -> Integer -> Dyadic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Rational -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
round (Rational -> Integer)
-> (Integer -> Rational) -> Integer -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int -> Integer
forall a. Bits a => Int -> a
bit (-Int
t') Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%)) ([Integer
1,Integer
3..] :: [Integer])
                  d :: [Dyadic]
d = (Dyadic -> Bool) -> [Dyadic] -> [Dyadic]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Dyadic -> Dyadic -> Bool
forall a. Eq a => a -> a -> Bool
/= Dyadic
0) ([Dyadic] -> [Dyadic])
-> ([Dyadic] -> [Dyadic]) -> [Dyadic] -> [Dyadic]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Dyadic -> Dyadic) -> [Dyadic] -> [Dyadic]
forall a b. (a -> b) -> [a] -> [b]
map (Int -> Dyadic -> Dyadic
shiftD Int
t') ([Dyadic] -> [Dyadic]) -> [Dyadic] -> [Dyadic]
forall a b. (a -> b) -> a -> b
$ (Dyadic -> Dyadic -> Dyadic) -> [Dyadic] -> [Dyadic] -> [Dyadic]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
(*) [Dyadic]
a [Dyadic]
c
                  e :: [Dyadic]
e = (Dyadic -> Bool) -> [Dyadic] -> [Dyadic]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Dyadic -> Dyadic -> Bool
forall a. Eq a => a -> a -> Bool
/= Dyadic
0) ([Dyadic] -> [Dyadic])
-> ([Dyadic] -> [Dyadic]) -> [Dyadic] -> [Dyadic]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Dyadic -> Dyadic) -> [Dyadic] -> [Dyadic]
forall a b. (a -> b) -> [a] -> [b]
map (Int -> Dyadic -> Dyadic
shiftD Int
t') ([Dyadic] -> [Dyadic]) -> [Dyadic] -> [Dyadic]
forall a b. (a -> b) -> a -> b
$ (Dyadic -> Dyadic -> Dyadic) -> [Dyadic] -> [Dyadic] -> [Dyadic]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
(*) [Dyadic]
b [Dyadic]
c
              in Int -> Dyadic -> Dyadic
shiftD Int
t (Dyadic -> Dyadic) -> Dyadic -> Dyadic
forall a b. (a -> b) -> a -> b
$ Dyadic
4 Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
* (Dyadic
4 Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
* [Dyadic] -> Dyadic
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Dyadic]
d Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
- [Dyadic] -> Dyadic
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Dyadic]
e)

-- |Compute dyadic values close to pi by Borwein's formula.
piBorweinD :: Int -> Dyadic
piBorweinD :: Int -> Dyadic
piBorweinD Int
t = let t' :: Int
t' = Int
tInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
10Int -> Int -> Int
forall a. Num a => a -> a -> a
-Integer -> Int
integerLog2 (Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Int
forall a. Num a => a -> a
abs Int
t))
                   s :: Dyadic
s = Int -> Dyadic -> Dyadic
sqrtD Int
t' Dyadic
2
                   d :: (Dyadic, Dyadic, Int)
d = [(Dyadic, Dyadic, Int)] -> (Dyadic, Dyadic, Int)
forall a. [a] -> a
head ([(Dyadic, Dyadic, Int)] -> (Dyadic, Dyadic, Int))
-> ([(Dyadic, Dyadic, Int)] -> [(Dyadic, Dyadic, Int)])
-> [(Dyadic, Dyadic, Int)]
-> (Dyadic, Dyadic, Int)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. ((Dyadic, Dyadic, Int) -> Bool)
-> [(Dyadic, Dyadic, Int)] -> [(Dyadic, Dyadic, Int)]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile (\(Dyadic
_,Dyadic
y,Int
_) -> Dyadic
y Dyadic -> Dyadic -> Bool
forall a. Eq a => a -> a -> Bool
/= Dyadic
0) ([(Dyadic, Dyadic, Int)] -> (Dyadic, Dyadic, Int))
-> [(Dyadic, Dyadic, Int)] -> (Dyadic, Dyadic, Int)
forall a b. (a -> b) -> a -> b
$ ((Dyadic, Dyadic, Int) -> (Dyadic, Dyadic, Int))
-> (Dyadic, Dyadic, Int) -> [(Dyadic, Dyadic, Int)]
forall a. (a -> a) -> a -> [a]
iterate (Int -> (Dyadic, Dyadic, Int) -> (Dyadic, Dyadic, Int)
f Int
t') (Dyadic
6Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
-Dyadic
4Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
*Dyadic
s,Dyadic
sDyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
-Dyadic
1,Int
1 :: Int)
               in Int -> Dyadic -> Dyadic
shiftD Int
t (Dyadic -> Dyadic)
-> ((Dyadic, Dyadic, Int) -> Dyadic)
-> (Dyadic, Dyadic, Int)
-> Dyadic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (\(Dyadic
a,Dyadic
_,Int
_) -> Int -> Dyadic -> Dyadic -> Dyadic
divD' Int
t' Dyadic
1 Dyadic
a) ((Dyadic, Dyadic, Int) -> Dyadic)
-> (Dyadic, Dyadic, Int) -> Dyadic
forall a b. (a -> b) -> a -> b
$ (Dyadic, Dyadic, Int)
d
    where
      f :: Int -> (Dyadic, Dyadic, Int) -> (Dyadic, Dyadic, Int)
      f :: Int -> (Dyadic, Dyadic, Int) -> (Dyadic, Dyadic, Int)
f Int
l (Dyadic
a,Dyadic
y,Int
k) = let u :: Dyadic
u = Int -> Dyadic -> Dyadic
sqrtD Int
l (Dyadic -> Dyadic) -> (Dyadic -> Dyadic) -> Dyadic -> Dyadic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Dyadic -> Dyadic
sqrtD Int
l (Dyadic -> Dyadic) -> Dyadic -> Dyadic
forall a b. (a -> b) -> a -> b
$ Dyadic
1Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
-Dyadic
yDyadic -> Int -> Dyadic
forall a b. (Num a, Integral b) => a -> b -> a
^(Int
4::Int)
                        y' :: Dyadic
y' = Int -> Dyadic -> Dyadic
shiftD Int
l (Dyadic -> Dyadic) -> Dyadic -> Dyadic
forall a b. (a -> b) -> a -> b
$ Dyadic -> Dyadic -> Dyadic
divD (Dyadic
1Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
-Dyadic
u) (Dyadic
1Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
+Dyadic
u)
                        a' :: Dyadic
a' = Int -> Dyadic -> Dyadic
shiftD Int
l (Dyadic -> Dyadic) -> Dyadic -> Dyadic
forall a b. (a -> b) -> a -> b
$ Dyadic
aDyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
*(Dyadic
1Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
+Dyadic
y')Dyadic -> Int -> Dyadic
forall a b. (Num a, Integral b) => a -> b -> a
^(Int
4::Int)Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
-(Int -> Integer
forall a. Bits a => Int -> a
bit (Int
2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)Integer -> Int -> Dyadic
:^Int
0)Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
*Dyadic
y'Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
*(Dyadic
1Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
+Dyadic
y'Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
+Dyadic
y'Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
*Dyadic
y')
                    in (Dyadic
a',Dyadic
y',Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1)

-- |Compute dyadic values close to ln 2.
ln2D :: Int -> Dyadic
ln2D :: Int -> Dyadic
ln2D Int
t = let t' :: Int
t' = Int
t Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
10 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Integer -> Int
integerLog2 (Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int -> Int
forall a. Num a => a -> a
abs Int
t))
             a :: [Dyadic]
a = (Integer -> Dyadic) -> [Integer] -> [Dyadic]
forall a b. (a -> b) -> [a] -> [b]
map ((Integer -> Int -> Dyadic
:^Int
t') (Integer -> Dyadic) -> (Integer -> Integer) -> Integer -> Dyadic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Rational -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
round (Rational -> Integer)
-> (Integer -> Rational) -> Integer -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int -> Integer
forall a. Bits a => Int -> a
bit (-Int
t') Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%)) ([Integer] -> [Dyadic]) -> [Integer] -> [Dyadic]
forall a b. (a -> b) -> a -> b
$ (Integer -> Integer) -> Integer -> [Integer]
forall a. (a -> a) -> a -> [a]
iterate (Integer
3Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*) (Integer
3 :: Integer)
             b :: [Dyadic]
b = (Integer -> Dyadic) -> [Integer] -> [Dyadic]
forall a b. (a -> b) -> [a] -> [b]
map ((Integer -> Int -> Dyadic
:^Int
t') (Integer -> Dyadic) -> (Integer -> Integer) -> Integer -> Dyadic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Rational -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
round (Rational -> Integer)
-> (Integer -> Rational) -> Integer -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int -> Integer
forall a. Bits a => Int -> a
bit (-Int
t') Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%)) ([Integer] -> [Dyadic]) -> [Integer] -> [Dyadic]
forall a b. (a -> b) -> a -> b
$ (Integer -> Integer) -> Integer -> [Integer]
forall a. (a -> a) -> a -> [a]
iterate (Integer
4Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*) (Integer
4 :: Integer)
             c :: [Dyadic]
c = (Integer -> Dyadic) -> [Integer] -> [Dyadic]
forall a b. (a -> b) -> [a] -> [b]
map ((Integer -> Int -> Dyadic
:^Int
t') (Integer -> Dyadic) -> (Integer -> Integer) -> Integer -> Dyadic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Rational -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
round (Rational -> Integer)
-> (Integer -> Rational) -> Integer -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Int -> Integer
forall a. Bits a => Int -> a
bit (-Int
t') Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%)) ([Integer
1..] :: [Integer])
             d :: [Dyadic]
d = (Dyadic -> Dyadic -> Dyadic) -> [Dyadic] -> [Dyadic] -> [Dyadic]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
(+) [Dyadic]
a [Dyadic]
b
             e :: [Dyadic]
e = (Dyadic -> Bool) -> [Dyadic] -> [Dyadic]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Dyadic -> Dyadic -> Bool
forall a. Eq a => a -> a -> Bool
/= Dyadic
0) ([Dyadic] -> [Dyadic])
-> ([Dyadic] -> [Dyadic]) -> [Dyadic] -> [Dyadic]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Dyadic -> Dyadic) -> [Dyadic] -> [Dyadic]
forall a b. (a -> b) -> [a] -> [b]
map (Int -> Dyadic -> Dyadic
shiftD Int
t') ([Dyadic] -> [Dyadic]) -> [Dyadic] -> [Dyadic]
forall a b. (a -> b) -> a -> b
$ (Dyadic -> Dyadic -> Dyadic) -> [Dyadic] -> [Dyadic] -> [Dyadic]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
(*) [Dyadic]
d [Dyadic]
c
         in Int -> Dyadic -> Dyadic
shiftD Int
t (Dyadic -> Dyadic) -> Dyadic -> Dyadic
forall a b. (a -> b) -> a -> b
$ [Dyadic] -> Dyadic
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [Dyadic]
e

logD :: Int -> Dyadic -> Dyadic
logD :: Int -> Dyadic -> Dyadic
logD Int
t x :: Dyadic
x@(Integer
a :^ Int
s) =
  if Integer
a Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
0 then String -> Dyadic
forall a. HasCallStack => String -> a
error String
"logD: Non-positive argument"
  else let t' :: Int
t' = Int
t 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
a 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
a) 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 Dyadic
x (-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
           z :: Dyadic
z = Int -> Dyadic
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
r Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
* Int -> Dyadic
ln2D Int
t'
       in  Int -> Dyadic -> Dyadic
shiftD Int
t (Dyadic -> Dyadic) -> (Dyadic -> Dyadic) -> Dyadic -> Dyadic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
+Dyadic
z) (Dyadic -> Dyadic) -> (Dyadic -> Dyadic) -> Dyadic -> Dyadic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (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

atanhD :: Int -> Dyadic -> Dyadic
atanhD :: Int -> Dyadic -> Dyadic
atanhD Int
t x :: Dyadic
x@(Integer
a :^ Int
s) =
  if Integer -> Int
integerLog2 (Integer -> Integer
forall a. Num a => a -> a
abs Integer
a) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
s Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0 then String -> Dyadic
forall a. HasCallStack => String -> a
error String
"atanhD: Argument outside domain, (-1,1)"
  else let 
    -- number of guard digits is 5+k where k depends on the precision [how do we know if this is enough?]
    t' :: Int
t' = Int
t Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
5 Int -> Int -> Int
forall a. Num a => a -> a -> a
- Integer -> Int
integerLog2 (Integer -> Integer
forall a. Num a => a -> a
abs (Integer -> Integer) -> Integer -> Integer
forall a b. (a -> b) -> a -> b
$ Int -> Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
t) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2
    g :: Dyadic -> Dyadic -> Dyadic
g Dyadic
_x Dyadic
_y = Int -> Dyadic -> Dyadic
shiftD Int
t' (Dyadic
_x Dyadic -> Dyadic -> Dyadic
forall a. Num a => a -> a -> a
* Dyadic
_y)
    x2 :: Dyadic
x2 = Dyadic -> Dyadic -> Dyadic
g Dyadic
x Dyadic
x
    b :: [Dyadic]
b = (Dyadic -> Dyadic) -> Dyadic -> [Dyadic]
forall a. (a -> a) -> a -> [a]
iterate (Dyadic -> Dyadic -> Dyadic
g Dyadic
x2) Dyadic
1
    c :: [Dyadic]
c = ((Integer -> Dyadic) -> [Integer] -> [Dyadic]
forall a b. (a -> b) -> [a] -> [b]
map (Int -> Dyadic -> Dyadic -> Dyadic
divD' Int
t' Dyadic
1 (Dyadic -> Dyadic) -> (Integer -> Dyadic) -> Integer -> Dyadic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Dyadic
forall a b. (Integral a, Num b) => a -> b
fromIntegral) [Integer
1,Integer
3..])
    in Int -> Dyadic -> Dyadic
shiftD Int
t (Dyadic -> Dyadic) -> ([Dyadic] -> Dyadic) -> [Dyadic] -> Dyadic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Dyadic -> Dyadic -> Dyadic
g Dyadic
x (Dyadic -> Dyadic) -> ([Dyadic] -> Dyadic) -> [Dyadic] -> Dyadic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [Dyadic] -> Dyadic
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([Dyadic] -> Dyadic)
-> ([Dyadic] -> [Dyadic]) -> [Dyadic] -> Dyadic
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Dyadic -> Bool) -> [Dyadic] -> [Dyadic]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Dyadic -> Dyadic -> Bool
forall a. Eq a => a -> a -> Bool
/= Dyadic
0) ([Dyadic] -> Dyadic) -> [Dyadic] -> Dyadic
forall a b. (a -> b) -> a -> b
$ (Dyadic -> Dyadic -> Dyadic) -> [Dyadic] -> [Dyadic] -> [Dyadic]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Dyadic -> Dyadic -> Dyadic
g [Dyadic]
b [Dyadic]
c

{-
agmD :: Int -> Dyadic -> Dyadic -> Dyadic
agmD t a b = let t' = t - 5
                 agmStep (c,d) = ((1:^(-1)) * (c+d), sqrtD t' (c*d))
                 close (c,d) = abs (c-d) < 1:^t
             in fst . last . takeWhile (not . close) $ iterate agmStep (shiftD t' a, shiftD t' b)

theta2D :: Int -> Dyadic -> Dyadic
theta2D t q = let q2 = shiftD t $ q*q
                  q4 = sqrtD t . sqrtD t $ q
                  a = iterate (shiftD t . (q2*)) 1
                  b = scanl1 (\c d -> shiftD t (c*d)) a
              in (2*) . shiftD t . (q4*) . sum $ takeWhile (/= 0) b

theta3D :: Int -> Dyadic -> Dyadic
theta3D t q = let q2 = shiftD t $ q*q
                  -- q4 = sqrtD t . sqrtD t $ q
                  a = iterate (shiftD t . (q2*)) q
                  b = scanl1 (\c d -> shiftD t (c*d)) a
              in (1+) . (2*) . shiftD t . sum $ takeWhile (/= 0) b

lnBigD :: Int -> Dyadic -> Dyadic
lnBigD t x = let t' = t-10-integerLog2 (fromIntegral (abs t))
                 t2 = theta2D t' x
                 t3 = theta3D t' x
                 a = agmD t' (t2*t2) (t3*t3)
             in shiftD t $ divD' t' (piBorweinD t') a
-}
{-
agm a b = let step (a,b) = (0.5 * (a+b), sqrt (a*b))
              close (a,b) = abs (a-b) < 1e-6
          in fst . last . takeWhile (not . close) $ iterate step (a, b)

theta2 q = let q2 = q^2
               q4 = sqrt . sqrt $ q
               a = iterate (q2*) 1
               b = scanl1 (*) a
           in (2*) . (q4*) . sum $ takeWhile (> 1e-6) b

theta3 q = let q2 = q^2
               q4 = sqrt . sqrt $ q
               a = iterate (q2*) q
               b = scanl1 (*) a
           in (1+) . (2*) . sum $ takeWhile (> 1e-6) b
-}

{-
-- checking that converging dyadic numbers only differ in last bit
check ((1:^_):(1:^_):_) = False
check ((1:^_):xs) = check xs
check ((-1:^_):(-1:^_):_) = False
check ((-1:^_):xs) = check xs
check ((0:^_):xs) = check xs
check ((_:^_):_) = False
check [] = True
check [_] = True
-}