-- ERA: Exact Real Arithmetic (version 1.0)

--

-- A tolerably efficient and possibly correct implementation of the computable

-- reals using Haskell 1.2.

--

-- David Lester, Department of Computer Science, Manchester University, M13 9PL.

--           (2000-2001)


module Data.Number.CReal(CReal, showCReal) where
import Data.Ratio
import Numeric(readFloat, readSigned)

-- |The 'CReal' type implements (constructive) real numbers.

--

-- Note that the comparison operations on 'CReal' may diverge

-- since it is (by necessity) impossible to implementent them

-- correctly and always terminating.

--

-- This implementation is really David Lester's ERA package.

data CReal = CR (Int -> Integer)

instance Eq  CReal where
  x :: CReal
x == :: CReal -> CReal -> Bool
== y :: CReal
y = Int -> Integer
s' (Int -> Int
digitsToBits Int
digits) Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 0 where (CR s' :: Int -> Integer
s') = CReal
xCReal -> CReal -> CReal
forall a. Num a => a -> a -> a
-CReal
y

instance Ord CReal where
  x :: CReal
x <= :: CReal -> CReal -> Bool
<= y :: CReal
y = Int -> Integer
s' (Int -> Int
digitsToBits Int
digits) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= 0 where (CR s' :: Int -> Integer
s') = CReal
xCReal -> CReal -> CReal
forall a. Num a => a -> a -> a
-CReal
y
  x :: CReal
x < :: CReal -> CReal -> Bool
<  y :: CReal
y = Int -> Integer
s' (Int -> Int
digitsToBits Int
digits) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<  0 where (CR s' :: Int -> Integer
s') = CReal
xCReal -> CReal -> CReal
forall a. Num a => a -> a -> a
-CReal
y
  x :: CReal
x >= :: CReal -> CReal -> Bool
>= y :: CReal
y = Int -> Integer
s' (Int -> Int
digitsToBits Int
digits) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>= 0 where (CR s' :: Int -> Integer
s') = CReal
xCReal -> CReal -> CReal
forall a. Num a => a -> a -> a
-CReal
y
  x :: CReal
x > :: CReal -> CReal -> Bool
>  y :: CReal
y = Int -> Integer
s' (Int -> Int
digitsToBits Int
digits) Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>  0 where (CR s' :: Int -> Integer
s') = CReal
xCReal -> CReal -> CReal
forall a. Num a => a -> a -> a
-CReal
y
  max :: CReal -> CReal -> CReal
max (CR x' :: Int -> Integer
x') (CR y' :: Int -> Integer
y') = (Int -> Integer) -> CReal
CR (\p :: Int
p -> Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
max (Int -> Integer
x' Int
p) (Int -> Integer
y' Int
p))
  min :: CReal -> CReal -> CReal
min (CR x' :: Int -> Integer
x') (CR y' :: Int -> Integer
y') = (Int -> Integer) -> CReal
CR (\p :: Int
p -> Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
min (Int -> Integer
x' Int
p) (Int -> Integer
y' Int
p))

instance Num CReal where
  (CR x' :: Int -> Integer
x') + :: CReal -> CReal -> CReal
+ (CR y' :: Int -> Integer
y') = (Int -> Integer) -> CReal
CR (\p :: Int
p -> Rational -> Integer
round_uk ((Int -> Integer
x' (Int
pInt -> Int -> Int
forall a. Num a => a -> a -> a
+2) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Int -> Integer
y' (Int
pInt -> Int -> Int
forall a. Num a => a -> a -> a
+2)) Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% 4))
  (CR x' :: Int -> Integer
x') * :: CReal -> CReal -> CReal
* (CR y' :: Int -> Integer
y') = (Int -> Integer) -> CReal
CR (\p :: Int
p -> Rational -> Integer
round_uk ((Int -> Integer
x' (Int
pInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
sy)Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Int -> Integer
y' (Int
pInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
sx)) Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% 2Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^(Int
pInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
sxInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
sy)))
                        where x0 :: Integer
x0 = Integer -> Integer
forall a. Num a => a -> a
abs (Int -> Integer
x' 0)Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+2; y0 :: Integer
y0 = Integer -> Integer
forall a. Num a => a -> a
abs (Int -> Integer
y' 0)Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+2
                              sx :: Int
sx = Integer -> Int -> Int
sizeinbase Integer
x0 2Int -> Int -> Int
forall a. Num a => a -> a -> a
+3; sy :: Int
sy = Integer -> Int -> Int
sizeinbase Integer
y0 2Int -> Int -> Int
forall a. Num a => a -> a -> a
+3
  negate :: CReal -> CReal
negate (CR x' :: Int -> Integer
x')    = (Int -> Integer) -> CReal
CR (\p :: Int
p -> Integer -> Integer
forall a. Num a => a -> a
negate (Int -> Integer
x' Int
p))
  abs :: CReal -> CReal
abs x :: CReal
x             = CReal -> CReal -> CReal
forall a. Ord a => a -> a -> a
max CReal
x (CReal -> CReal
forall a. Num a => a -> a
negate CReal
x)
  signum :: CReal -> CReal
signum (CR x' :: Int -> Integer
x')    = Integer -> CReal
forall a. Num a => Integer -> a
fromInteger (Integer -> Integer
forall a. Num a => a -> a
signum (Int -> Integer
x' (Int -> Int
digitsToBits Int
digits)))
  fromInteger :: Integer -> CReal
fromInteger n :: Integer
n     = (Int -> Integer) -> CReal
CR (\p :: Int
p -> Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*2Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^Int
p)

instance Fractional CReal where
  recip :: CReal -> CReal
recip (CR x' :: Int -> Integer
x') = (Int -> Integer) -> CReal
CR (\p :: Int
p -> let s :: Int
s = [Int] -> Int
forall a. [a] -> a
head [Int
n | Int
n <- [0..], 3 Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer -> Integer
forall a. Num a => a -> a
abs (Int -> Integer
x' Int
n)]
                              in Rational -> Integer
round_uk (2Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^(2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
pInt -> Int -> Int
forall a. Num a => a -> a -> a
+2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
+2) Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% (Int -> Integer
x' (Int
pInt -> Int -> Int
forall a. Num a => a -> a -> a
+2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
sInt -> Int -> Int
forall a. Num a => a -> a -> a
+2))))
  fromRational :: Rational -> CReal
fromRational x :: Rational
x = Integer -> CReal
forall a. Num a => Integer -> a
fromInteger (Rational -> Integer
forall a. Ratio a -> a
numerator Rational
x) CReal -> CReal -> CReal
forall a. Fractional a => a -> a -> a
/ Integer -> CReal
forall a. Num a => Integer -> a
fromInteger (Rational -> Integer
forall a. Ratio a -> a
denominator Rational
x)

-- two useful scaling functions:


div2n :: CReal -> Int -> CReal
div2n :: CReal -> Int -> CReal
div2n (CR x' :: Int -> Integer
x') n :: Int
n = (Int -> Integer) -> CReal
CR (\p :: Int
p -> if Int
p Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
n then Int -> Integer
x' (Int
pInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
n) else Rational -> Integer
round_uk (Int -> Integer
x' Int
p Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% 2Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^Int
n))

mul2n :: CReal -> Int -> CReal
mul2n :: CReal -> Int -> CReal
mul2n (CR x' :: Int -> Integer
x') n :: Int
n = (Int -> Integer) -> CReal
CR (\p :: Int
p -> Int -> Integer
x' (Int
pInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
n))

-- transcendental functions (mostly range reductions):


instance Floating CReal where
  pi :: CReal
pi = 16 CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
* CReal -> CReal
forall a. Floating a => a -> a
atan (Rational -> CReal
forall a. Fractional a => Rational -> a
fromRational (1 Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% 5)) 
                CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
- 4 CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
* CReal -> CReal
forall a. Floating a => a -> a
atan (Rational -> CReal
forall a. Fractional a => Rational -> a
fromRational (1 Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% 239))
  sqrt :: CReal -> CReal
sqrt x :: CReal
x  = (Int -> Integer) -> CReal
CR (\p :: Int
p -> Integer -> Integer
floorsqrt (Int -> Integer
x' (2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Int
p))) where (CR x' :: Int -> Integer
x') = CReal
x

  log :: CReal -> CReal
log x :: CReal
x   = if Integer
t Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< 0 then [Char] -> CReal
forall a. HasCallStack => [Char] -> a
error "log of negative number\n" else
            if Integer
t Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< 4 then - CReal -> CReal
forall a. Floating a => a -> a
log (CReal -> CReal
forall a. Fractional a => a -> a
recip CReal
x)                  else
            if Integer
t Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< 8 then CReal -> CReal
log_dr CReal
x                         else
            {- 7 < t -}   CReal -> CReal
log_dr (CReal -> Int -> CReal
div2n CReal
x Int
n) CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
+ Int -> CReal
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
* CReal
log2
            where (CR x' :: Int -> Integer
x') = CReal
x; t :: Integer
t = Int -> Integer
x' 2; n :: Int
n = Integer -> Int -> Int
sizeinbase Integer
t 2 Int -> Int -> Int
forall a. Num a => a -> a -> a
- 3
  exp :: CReal -> CReal
exp x :: CReal
x   = if Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< 0 then CReal -> Int -> CReal
div2n (CReal -> CReal
exp_dr CReal
s) (Integer -> Int
forall a. Num a => Integer -> a
fromInteger (-Integer
n)) else
            if Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> 0 then CReal -> Int -> CReal
mul2n (CReal -> CReal
exp_dr CReal
s) (Integer -> Int
forall a. Num a => Integer -> a
fromInteger Integer
n) else CReal -> CReal
exp_dr CReal
s
            where (CR u' :: Int -> Integer
u') = CReal
xCReal -> CReal -> CReal
forall a. Fractional a => a -> a -> a
/CReal
log2; n :: Integer
n = Int -> Integer
u' 0; s :: CReal
s = CReal
xCReal -> CReal -> CReal
forall a. Num a => a -> a -> a
-Integer -> CReal
forall a. Num a => Integer -> a
fromInteger Integer
nCReal -> CReal -> CReal
forall a. Num a => a -> a -> a
*CReal
log2
  sin :: CReal -> CReal
sin x :: CReal
x   = if Integer
n Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 0 then CReal -> CReal
sin_dr CReal
y                           else
            if Integer
n Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 1 then CReal
sqrt1By2 CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
* (CReal -> CReal
cos_dr CReal
y CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
+ CReal -> CReal
sin_dr CReal
y)   else
            if Integer
n Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 2 then CReal -> CReal
cos_dr CReal
y                           else
            if Integer
n Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 3 then CReal
sqrt1By2 CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
* (CReal -> CReal
cos_dr CReal
y CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
- CReal -> CReal
sin_dr CReal
y)   else
            if Integer
n Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 4 then - CReal -> CReal
sin_dr CReal
y                         else
            if Integer
n Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 5 then - CReal
sqrt1By2 CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
* (CReal -> CReal
cos_dr CReal
y CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
+ CReal -> CReal
sin_dr CReal
y) else
            if Integer
n Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 6 then - CReal -> CReal
cos_dr CReal
y                         else
            {- n == 7 -}   - CReal
sqrt1By2 CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
* (CReal -> CReal
cos_dr CReal
y CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
- CReal -> CReal
sin_dr CReal
y)
            where (CR z' :: Int -> Integer
z') = CReal
xCReal -> CReal -> CReal
forall a. Fractional a => a -> a -> a
/CReal
piBy4; s :: Integer
s = Rational -> Integer
round_uk (Int -> Integer
z' 2 Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% 4); n :: Integer
n = Integer
s Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` 8
                  y :: CReal
y = CReal
x CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
- CReal
piBy4 CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
* Integer -> CReal
forall a. Num a => Integer -> a
fromInteger Integer
s
  cos :: CReal -> CReal
cos x :: CReal
x   = if Integer
n Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 0 then CReal -> CReal
cos_dr CReal
y                           else
            if Integer
n Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 1 then CReal
sqrt1By2 CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
* (CReal -> CReal
cos_dr CReal
y CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
- CReal -> CReal
sin_dr CReal
y)   else
            if Integer
n Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 2 then - CReal -> CReal
sin_dr CReal
y                           else
            if Integer
n Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 3 then - CReal
sqrt1By2 CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
* (CReal -> CReal
cos_dr CReal
y CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
+ CReal -> CReal
sin_dr CReal
y)   else
            if Integer
n Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 4 then - CReal -> CReal
cos_dr CReal
y                         else
            if Integer
n Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 5 then - CReal
sqrt1By2 CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
* (CReal -> CReal
cos_dr CReal
y CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
- CReal -> CReal
sin_dr CReal
y) else
            if Integer
n Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 6 then CReal -> CReal
sin_dr CReal
y                         else
            {- n == 7 -}   CReal
sqrt1By2 CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
* (CReal -> CReal
cos_dr CReal
y CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
+ CReal -> CReal
sin_dr CReal
y)
            where (CR z' :: Int -> Integer
z') = CReal
xCReal -> CReal -> CReal
forall a. Fractional a => a -> a -> a
/CReal
piBy4; s :: Integer
s = Rational -> Integer
round_uk (Int -> Integer
z' 2 Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% 4); n :: Integer
n = Integer
s Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`mod` 8
                  y :: CReal
y = CReal
x CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
- CReal
piBy4 CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
* Integer -> CReal
forall a. Num a => Integer -> a
fromInteger Integer
s
  atan :: CReal -> CReal
atan x :: CReal
x  = if Integer
t Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<  -5 then CReal -> CReal
atan_dr (CReal -> CReal
forall a. Num a => a -> a
negate (CReal -> CReal
forall a. Fractional a => a -> a
recip CReal
x)) CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
- CReal
piBy2 else
            if Integer
t Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== -4 then -CReal
piBy4 CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
- CReal -> CReal
atan_dr (CReal
xp1CReal -> CReal -> CReal
forall a. Fractional a => a -> a -> a
/CReal
xm1)         else
            if Integer
t Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<   4 then CReal -> CReal
atan_dr CReal
x                          else
            if Integer
t Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
==  4 then CReal
piBy4 CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
+ CReal -> CReal
atan_dr (CReal
xm1CReal -> CReal -> CReal
forall a. Fractional a => a -> a -> a
/CReal
xp1)          else
            {- t >   4 -}   CReal
piBy2 CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
- CReal -> CReal
atan_dr (CReal -> CReal
forall a. Fractional a => a -> a
recip CReal
x)
            where (CR x' :: Int -> Integer
x') = CReal
x; t :: Integer
t = Int -> Integer
x' 2
                  xp1 :: CReal
xp1 = CReal
xCReal -> CReal -> CReal
forall a. Num a => a -> a -> a
+1; xm1 :: CReal
xm1 = CReal
xCReal -> CReal -> CReal
forall a. Num a => a -> a -> a
-1
  asin :: CReal -> CReal
asin x :: CReal
x  = if Integer
x0 Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>  0 then CReal
forall a. Floating a => a
pi CReal -> CReal -> CReal
forall a. Fractional a => a -> a -> a
/ 2 CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
- CReal -> CReal
forall a. Floating a => a -> a
atan (CReal
sCReal -> CReal -> CReal
forall a. Fractional a => a -> a -> a
/CReal
x) else
            if Integer
x0 Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 0 then CReal -> CReal
forall a. Floating a => a -> a
atan (CReal
xCReal -> CReal -> CReal
forall a. Fractional a => a -> a -> a
/CReal
s)                      else
            {- x0 <  0 -}   - CReal -> CReal
forall a. Floating a => a -> a
atan (CReal
sCReal -> CReal -> CReal
forall a. Fractional a => a -> a -> a
/CReal
x) CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
- CReal
forall a. Floating a => a
pi CReal -> CReal -> CReal
forall a. Fractional a => a -> a -> a
/ 2
            where (CR x' :: Int -> Integer
x') = CReal
x; x0 :: Integer
x0 = Int -> Integer
x' 0; s :: CReal
s = CReal -> CReal
forall a. Floating a => a -> a
sqrt (1 CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
- CReal
xCReal -> CReal -> CReal
forall a. Num a => a -> a -> a
*CReal
x)
  acos :: CReal -> CReal
acos x :: CReal
x  = CReal
forall a. Floating a => a
pi CReal -> CReal -> CReal
forall a. Fractional a => a -> a -> a
/ 2 CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
- CReal -> CReal
forall a. Floating a => a -> a
asin CReal
x
  sinh :: CReal -> CReal
sinh x :: CReal
x  = (CReal
y CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
- CReal -> CReal
forall a. Fractional a => a -> a
recip CReal
y) CReal -> CReal -> CReal
forall a. Fractional a => a -> a -> a
/ 2 where y :: CReal
y = CReal -> CReal
forall a. Floating a => a -> a
exp CReal
x
  cosh :: CReal -> CReal
cosh x :: CReal
x  = (CReal
y CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
+ CReal -> CReal
forall a. Fractional a => a -> a
recip CReal
y) CReal -> CReal -> CReal
forall a. Fractional a => a -> a -> a
/ 2 where y :: CReal
y = CReal -> CReal
forall a. Floating a => a -> a
exp CReal
x
  tanh :: CReal -> CReal
tanh x :: CReal
x  = (CReal
y CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
- CReal
y') CReal -> CReal -> CReal
forall a. Fractional a => a -> a -> a
/ (CReal
y CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
+ CReal
y') where y :: CReal
y = CReal -> CReal
forall a. Floating a => a -> a
exp CReal
x; y' :: CReal
y' = CReal -> CReal
forall a. Fractional a => a -> a
recip CReal
y
  asinh :: CReal -> CReal
asinh x :: CReal
x = CReal -> CReal
forall a. Floating a => a -> a
log (CReal
x CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
+ CReal -> CReal
forall a. Floating a => a -> a
sqrt (CReal
xCReal -> CReal -> CReal
forall a. Num a => a -> a -> a
*CReal
x CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
+ 1))
  acosh :: CReal -> CReal
acosh x :: CReal
x = CReal -> CReal
forall a. Floating a => a -> a
log (CReal
x CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
+ CReal -> CReal
forall a. Floating a => a -> a
sqrt (CReal
xCReal -> CReal -> CReal
forall a. Num a => a -> a -> a
*CReal
x CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
- 1))
  atanh :: CReal -> CReal
atanh x :: CReal
x = CReal -> CReal
forall a. Floating a => a -> a
log ((1 CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
+ CReal
x) CReal -> CReal -> CReal
forall a. Fractional a => a -> a -> a
/ (1 CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
- CReal
x)) CReal -> CReal -> CReal
forall a. Fractional a => a -> a -> a
/ 2


acc_seq :: (Rational -> Integer -> Rational) -> [Rational]
acc_seq :: (Rational -> Integer -> Rational) -> [Rational]
acc_seq f :: Rational -> Integer -> Rational
f = (Rational -> Integer -> Rational)
-> Rational -> [Integer] -> [Rational]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl Rational -> Integer -> Rational
f (1 Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% 1) [1..]

exp_dr :: CReal -> CReal
exp_dr :: CReal -> CReal
exp_dr = [Rational] -> (Int -> Int) -> CReal -> CReal
power_series ((Rational -> Integer -> Rational) -> [Rational]
acc_seq (\a :: Rational
a n :: Integer
n -> Rational
aRational -> Rational -> Rational
forall a. Num a => a -> a -> a
*(1 Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% Integer
n))) Int -> Int
forall a. a -> a
id

log_dr :: CReal -> CReal
log_dr :: CReal -> CReal
log_dr x :: CReal
x = CReal
y CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
* CReal -> CReal
log_drx CReal
y where y :: CReal
y = (CReal
x CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
- 1) CReal -> CReal -> CReal
forall a. Fractional a => a -> a -> a
/ CReal
x

log_drx :: CReal -> CReal
log_drx :: CReal -> CReal
log_drx = [Rational] -> (Int -> Int) -> CReal -> CReal
power_series [1 Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% Integer
n | Integer
n <- [1..]] (Int -> Int -> Int
forall a. Num a => a -> a -> a
+1)

sin_dr :: CReal -> CReal
sin_dr :: CReal -> CReal
sin_dr x :: CReal
x = CReal
xCReal -> CReal -> CReal
forall a. Num a => a -> a -> a
*[Rational] -> (Int -> Int) -> CReal -> CReal
power_series ((Rational -> Integer -> Rational) -> [Rational]
acc_seq (\a :: Rational
a n :: Integer
n -> -Rational
aRational -> Rational -> Rational
forall a. Num a => a -> a -> a
*(1 Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% (2Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*(2Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+1))))) Int -> Int
forall a. a -> a
id (CReal
xCReal -> CReal -> CReal
forall a. Num a => a -> a -> a
*CReal
x)

cos_dr :: CReal -> CReal
cos_dr :: CReal -> CReal
cos_dr x :: CReal
x = [Rational] -> (Int -> Int) -> CReal -> CReal
power_series ((Rational -> Integer -> Rational) -> [Rational]
acc_seq (\a :: Rational
a n :: Integer
n -> -Rational
aRational -> Rational -> Rational
forall a. Num a => a -> a -> a
*(1 Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% (2Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*(2Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-1))))) Int -> Int
forall a. a -> a
id (CReal
xCReal -> CReal -> CReal
forall a. Num a => a -> a -> a
*CReal
x)

atan_dr :: CReal -> CReal
atan_dr :: CReal -> CReal
atan_dr x :: CReal
x = (CReal
xCReal -> CReal -> CReal
forall a. Fractional a => a -> a -> a
/CReal
y) CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
* CReal -> CReal
atan_drx ((CReal
xCReal -> CReal -> CReal
forall a. Num a => a -> a -> a
*CReal
x)CReal -> CReal -> CReal
forall a. Fractional a => a -> a -> a
/CReal
y) where y :: CReal
y = CReal
xCReal -> CReal -> CReal
forall a. Num a => a -> a -> a
*CReal
xCReal -> CReal -> CReal
forall a. Num a => a -> a -> a
+1

atan_drx :: CReal -> CReal
atan_drx :: CReal -> CReal
atan_drx = [Rational] -> (Int -> Int) -> CReal -> CReal
power_series ((Rational -> Integer -> Rational) -> [Rational]
acc_seq (\a :: Rational
a n :: Integer
n -> Rational
aRational -> Rational -> Rational
forall a. Num a => a -> a -> a
*((2Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
n) Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% (2Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+1)))) (Int -> Int -> Int
forall a. Num a => a -> a -> a
+1)

-- power_series takes as arguments:

--   a (rational) list of the coefficients of the power series

--   a function from the desired accuracy to the number of terms needed

--   the argument x


power_series :: [Rational] -> (Int -> Int) -> CReal -> CReal
power_series :: [Rational] -> (Int -> Int) -> CReal -> CReal
power_series ps :: [Rational]
ps terms :: Int -> Int
terms (CR x' :: Int -> Integer
x')
  = (Int -> Integer) -> CReal
CR (\p :: Int
p -> let t :: Int
t = Int -> Int
terms Int
p; l2t :: Int
l2t = 2Int -> Int -> Int
forall a. Num a => a -> a -> a
*Integer -> Int -> Int
sizeinbase (Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
tInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+1) 2Int -> Int -> Int
forall a. Num a => a -> a -> a
+6; p' :: Int
p' = Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
l2t
                  xr :: Integer
xr = Int -> Integer
x' Int
p'; xn :: Integer
xn = 2Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^Int
p'; g :: Integer -> Integer
g yn :: Integer
yn = Rational -> Integer
round_uk ((Integer
ynInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
xr) Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% (2Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^Int
p'))
               in Rational -> Integer
round_uk ([Integer] -> [Rational] -> Integer
accumulate ((Integer -> Integer) -> Integer -> [Integer]
forall a. (a -> a) -> a -> [a]
iterate Integer -> Integer
g Integer
xn) (Int -> [Rational] -> [Rational]
forall a. Int -> [a] -> [a]
take Int
t [Rational]
ps) Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% (2Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^Int
l2t)))
    where accumulate :: [Integer] -> [Rational] -> Integer
accumulate _      []     = 0
          accumulate []     _      = [Char] -> Integer
forall a. HasCallStack => [Char] -> a
error "CReal.power_series.accumulate"
          accumulate (x :: Integer
x:xs :: [Integer]
xs) (c :: Rational
c:cs :: [Rational]
cs) = let t :: Integer
t = Rational -> Integer
round_uk (Rational
cRational -> Rational -> Rational
forall a. Num a => a -> a -> a
*(Integer
x Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% 1)) in
                                     if Integer
t Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 0 then 0 else Integer
t Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ [Integer] -> [Rational] -> Integer
accumulate [Integer]
xs [Rational]
cs

-- Some useful constants:


piBy2 :: CReal
piBy2 :: CReal
piBy2 = CReal -> Int -> CReal
div2n CReal
forall a. Floating a => a
pi 1

piBy4 :: CReal
piBy4 :: CReal
piBy4 = CReal -> Int -> CReal
div2n CReal
forall a. Floating a => a
pi 2

log2 :: CReal
log2 :: CReal
log2 = CReal -> Int -> CReal
div2n (CReal -> CReal
log_drx (CReal -> CReal
forall a. Fractional a => a -> a
recip 2)) 1

sqrt1By2 :: CReal
sqrt1By2 :: CReal
sqrt1By2 = CReal -> CReal
forall a. Floating a => a -> a
sqrt (CReal -> CReal
forall a. Fractional a => a -> a
recip 2)

instance Enum CReal where
  toEnum :: Int -> CReal
toEnum i :: Int
i         = Int -> CReal
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i
  fromEnum :: CReal -> Int
fromEnum _       = [Char] -> Int
forall a. HasCallStack => [Char] -> a
error "Cannot fromEnum CReal"
  enumFrom :: CReal -> [CReal]
enumFrom         = (CReal -> CReal) -> CReal -> [CReal]
forall a. (a -> a) -> a -> [a]
iterate (CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
+ 1)
  enumFromTo :: CReal -> CReal -> [CReal]
enumFromTo n :: CReal
n e :: CReal
e   = (CReal -> Bool) -> [CReal] -> [CReal]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (CReal -> CReal -> Bool
forall a. Ord a => a -> a -> Bool
<= CReal
e) ([CReal] -> [CReal]) -> [CReal] -> [CReal]
forall a b. (a -> b) -> a -> b
$ (CReal -> CReal) -> CReal -> [CReal]
forall a. (a -> a) -> a -> [a]
iterate (CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
+ 1)CReal
n
  enumFromThen :: CReal -> CReal -> [CReal]
enumFromThen n :: CReal
n m :: CReal
m = (CReal -> CReal) -> CReal -> [CReal]
forall a. (a -> a) -> a -> [a]
iterate (CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
+(CReal
mCReal -> CReal -> CReal
forall a. Num a => a -> a -> a
-CReal
n)) CReal
n
  enumFromThenTo :: CReal -> CReal -> CReal -> [CReal]
enumFromThenTo n :: CReal
n m :: CReal
m e :: CReal
e = if CReal
m CReal -> CReal -> Bool
forall a. Ord a => a -> a -> Bool
>= CReal
n then (CReal -> Bool) -> [CReal] -> [CReal]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (CReal -> CReal -> Bool
forall a. Ord a => a -> a -> Bool
<= CReal
e) ([CReal] -> [CReal]) -> [CReal] -> [CReal]
forall a b. (a -> b) -> a -> b
$ (CReal -> CReal) -> CReal -> [CReal]
forall a. (a -> a) -> a -> [a]
iterate (CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
+(CReal
mCReal -> CReal -> CReal
forall a. Num a => a -> a -> a
-CReal
n)) CReal
n
                          else (CReal -> Bool) -> [CReal] -> [CReal]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (CReal -> CReal -> Bool
forall a. Ord a => a -> a -> Bool
>= CReal
e) ([CReal] -> [CReal]) -> [CReal] -> [CReal]
forall a b. (a -> b) -> a -> b
$ (CReal -> CReal) -> CReal -> [CReal]
forall a. (a -> a) -> a -> [a]
iterate (CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
+(CReal
mCReal -> CReal -> CReal
forall a. Num a => a -> a -> a
-CReal
n)) CReal
n

instance Real CReal where
 -- toRational x@(CR x') = x' n % 2^n where n = digitsToBits digits

  toRational :: CReal -> Rational
toRational _ = [Char] -> Rational
forall a. HasCallStack => [Char] -> a
error "CReal.toRational"
instance RealFrac CReal where
  properFraction :: CReal -> (b, CReal)
properFraction x :: CReal
x@(CR x' :: Int -> Integer
x') = (Integer -> b
forall a. Num a => Integer -> a
fromInteger Integer
n, CReal
x CReal -> CReal -> CReal
forall a. Num a => a -> a -> a
- Integer -> CReal
forall a. Num a => Integer -> a
fromInteger Integer
n) where n :: Integer
n = Int -> Integer
x' 0

instance RealFloat CReal where
  floatRadix :: CReal -> Integer
floatRadix _ = [Char] -> Integer
forall a. HasCallStack => [Char] -> a
error "CReal.floatRadix"
  floatDigits :: CReal -> Int
floatDigits _ = [Char] -> Int
forall a. HasCallStack => [Char] -> a
error "CReal.floatDigits"
  floatRange :: CReal -> (Int, Int)
floatRange _ = [Char] -> (Int, Int)
forall a. HasCallStack => [Char] -> a
error "CReal.floatRange"
  decodeFloat :: CReal -> (Integer, Int)
decodeFloat _ = [Char] -> (Integer, Int)
forall a. HasCallStack => [Char] -> a
error "CReal.decodeFloat"
  encodeFloat :: Integer -> Int -> CReal
encodeFloat _ _ = [Char] -> CReal
forall a. HasCallStack => [Char] -> a
error "CReal.encodeFloat"
  exponent :: CReal -> Int
exponent _ = 0
  scaleFloat :: Int -> CReal -> CReal
scaleFloat 0 x :: CReal
x = CReal
x
  significand :: CReal -> CReal
significand x :: CReal
x = CReal
x
  isNaN :: CReal -> Bool
isNaN _ = Bool
False
  isInfinite :: CReal -> Bool
isInfinite _ = Bool
False
  isDenormalized :: CReal -> Bool
isDenormalized _ = Bool
False
  isNegativeZero :: CReal -> Bool
isNegativeZero _ = Bool
False
  isIEEE :: CReal -> Bool
isIEEE _ = Bool
False

-- printing and reading the reals:


-- |The 'showCReal' function connverts a 'CReal' to a 'String'.

showCReal :: Int                -- ^ The number of decimals

          -> CReal              -- ^ The real number

          -> String             -- ^ The resulting string

showCReal :: Int -> CReal -> [Char]
showCReal d :: Int
d (CR x' :: Int -> Integer
x')
  = (if Bool
s then "-" else "") [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
zs [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ (if Int
d Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= 0 then '.'Char -> [Char] -> [Char]
forall a. a -> [a] -> [a]
:[Char]
fs' else "")
    where b :: Int
b  = Int -> Int
digitsToBits Int
d
          n :: Integer
n  = Int -> Integer
x' Int
b
          ds :: [Char]
ds = Integer -> [Char]
forall a. Show a => a -> [Char]
show (Rational -> Integer
round_uk ((Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*10Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^Int
d) Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% 2Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^Int
b))
          (s :: Bool
s,ds' :: [Char]
ds') = let sgn :: Bool
sgn = [Char] -> Char
forall a. [a] -> a
head [Char]
ds Char -> Char -> Bool
forall a. Eq a => a -> a -> Bool
== '-' in (Bool
sgn, if Bool
sgn then [Char] -> [Char]
forall a. [a] -> [a]
tail [Char]
ds else [Char]
ds)
          ds'' :: [Char]
ds'' = Int -> [Char] -> [Char]
forall a. Int -> [a] -> [a]
take (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max (Int
dInt -> Int -> Int
forall a. Num a => a -> a -> a
+1Int -> Int -> Int
forall a. Num a => a -> a -> a
-[Char] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Char]
ds') 0) (Char -> [Char]
forall a. a -> [a]
repeat '0') [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
ds'
          (zs :: [Char]
zs,fs :: [Char]
fs) = Int -> [Char] -> ([Char], [Char])
forall a. Int -> [a] -> ([a], [a])
splitAt ([Char] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [Char]
ds'' Int -> Int -> Int
forall a. Num a => a -> a -> a
-Int
d) [Char]
ds''
          fs' :: [Char]
fs' = case [Char] -> [Char]
forall a. [a] -> [a]
reverse ([Char] -> [Char]) -> [Char] -> [Char]
forall a b. (a -> b) -> a -> b
$ (Char -> Bool) -> [Char] -> [Char]
forall a. (a -> Bool) -> [a] -> [a]
dropWhile (Char -> Char -> Bool
forall a. Eq a => a -> a -> Bool
== '0') ([Char] -> [Char]) -> [Char] -> [Char]
forall a b. (a -> b) -> a -> b
$ [Char] -> [Char]
forall a. [a] -> [a]
reverse [Char]
fs of
                "" -> "0"
                xs :: [Char]
xs -> [Char]
xs

digitsToBits :: Int -> Int
digitsToBits :: Int -> Int
digitsToBits d :: Int
d = Double -> Int
forall a b. (RealFrac a, Integral b) => a -> b
ceiling (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
d Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double -> Double -> Double
forall a. Floating a => a -> a -> a
logBase 2.0 10.0 :: Double)) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ 4

digits :: Int
digits :: Int
digits = 40

instance Read CReal where
  readsPrec :: Int -> ReadS CReal
readsPrec _p :: Int
_p = ReadS CReal -> ReadS CReal
forall a. Real a => ReadS a -> ReadS a
readSigned ReadS CReal
forall a. RealFrac a => ReadS a
readFloat

instance Show CReal where
  showsPrec :: Int -> CReal -> [Char] -> [Char]
showsPrec p :: Int
p x :: CReal
x = let xs :: [Char]
xs = Int -> CReal -> [Char]
showCReal Int
digits CReal
x in
                  if [Char] -> Char
forall a. [a] -> a
head [Char]
xs Char -> Char -> Bool
forall a. Eq a => a -> a -> Bool
== '-' then Bool -> ([Char] -> [Char]) -> [Char] -> [Char]
showParen (Int
p Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> 6) ([Char] -> [Char] -> [Char]
showString [Char]
xs)
                                    else [Char] -> [Char] -> [Char]
showString [Char]
xs

-- GMP functions not provided by Haskell


sizeinbase :: Integer -> Int -> Int
sizeinbase :: Integer -> Int -> Int
sizeinbase i :: Integer
i b :: Int
b = Integer -> Int
forall a. Num a => Integer -> a
f (Integer -> Integer
forall a. Num a => a -> a
abs Integer
i)
                 where f :: Integer -> p
f n :: Integer
n = if Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= 1 then 1 else 1 p -> p -> p
forall a. Num a => a -> a -> a
+ Integer -> p
f (Integer
n Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div` Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
b)

floorsqrt :: Integer -> Integer
floorsqrt :: Integer -> Integer
floorsqrt x :: Integer
x = (Integer -> Bool) -> (Integer -> Integer) -> Integer -> Integer
forall a. (a -> Bool) -> (a -> a) -> a -> a
until Integer -> Bool
satisfy Integer -> Integer
forall b. Integral b => Integer -> b
improve Integer
x
              where improve :: Integer -> b
improve y :: Integer
y = Rational -> b
forall a b. (RealFrac a, Integral b) => a -> b
floor ((Integer
yInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
yInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
x) Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% (2Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
y))
                    satisfy :: Integer -> Bool
satisfy y :: Integer
y = Integer
yInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
y Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
x Bool -> Bool -> Bool
&& Integer
x Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= (Integer
yInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+1)Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*(Integer
yInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+1)

round_uk :: Rational -> Integer
round_uk :: Rational -> Integer
round_uk x :: Rational
x = Rational -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
floor (Rational
xRational -> Rational -> Rational
forall a. Num a => a -> a -> a
+1 Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% 2)