{-# LANGUAGE DataKinds #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE MagicHash #-}
{-# LANGUAGE MultiWayIf #-}
{-# LANGUAGE PostfixOperators #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE BangPatterns #-}

-----------------------------------------------------------------------------
-- | This module exports a bunch of utilities for working inside the CReal
-- datatype. One should be careful to maintain the CReal invariant when using
-- these functions
----------------------------------------------------------------------------
module Data.CReal.Internal
  (
    -- * The CReal type
    CReal(..)
    -- ** Memoization
  , Cache(..)
    -- ** Simple utilities
  , atPrecision
  , crealPrecision

    -- * More efficient variants of common functions
    -- Note that the preconditions to these functions are not checked
    -- ** Additive
  , plusInteger
    -- ** Multiplicative
  , mulBounded
  , (.*.)
  , mulBoundedL
  , (.*)
  , (*.)
  , recipBounded
  , shiftL
  , shiftR
  , square
  , squareBounded

    -- ** Exponential
  , expBounded
  , expPosNeg
  , logBounded

    -- ** Trigonometric
  , atanBounded
  , sinBounded
  , cosBounded

    -- * Utilities for operating inside CReals
  , crMemoize
  , powerSeries
  , alternateSign

    -- ** Integer operations
  , (/.)
  , (/^)
  , log2
  , log10
  , isqrt

    -- * Utilities for converting CReals to Strings
  , showAtPrecision
  , decimalDigitsAtPrecision
  , rationalToDecimal
  ) where

import Data.List (scanl')
import qualified Data.Bits as B
import Data.Bits hiding (shiftL, shiftR)
import GHC.Base (Int(..))
import GHC.Integer.Logarithms (integerLog2#, integerLogBase#)
import GHC.Real (Ratio(..), (%))
import GHC.TypeLits
import Text.Read
import qualified Text.Read.Lex as L
import System.Random (Random(..))
import Control.Concurrent.MVar
import Control.Exception
import System.IO.Unsafe (unsafePerformIO)

{-# ANN module ("HLint: ignore Reduce duplication" :: String) #-}

-- $setup
-- >>> :set -XDataKinds
-- >>> :set -XPostfixOperators

default ()

-- | The Cache type represents a way to memoize a `CReal`. It holds the largest
-- precision the number has been evaluated that, as well as the value. Rounding
-- it down gives the value for lower numbers.
data Cache
  = Never
  | Current {-# UNPACK #-} !Int !Integer
  deriving (Int -> Cache -> ShowS
[Cache] -> ShowS
Cache -> String
(Int -> Cache -> ShowS)
-> (Cache -> String) -> ([Cache] -> ShowS) -> Show Cache
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [Cache] -> ShowS
$cshowList :: [Cache] -> ShowS
show :: Cache -> String
$cshow :: Cache -> String
showsPrec :: Int -> Cache -> ShowS
$cshowsPrec :: Int -> Cache -> ShowS
Show)

-- | The type CReal represents a fast binary Cauchy sequence. This is a Cauchy
-- sequence with the invariant that the pth element divided by 2^p will be
-- within 2^-p of the true value. Internally this sequence is represented as a
-- function from Ints to Integers, as well as an `MVar` to hold the highest
-- precision cached value.
data CReal (n :: Nat) = CR {-# UNPACK #-} !(MVar Cache) (Int -> Integer)

-- | 'crMemoize' takes a fast binary Cauchy sequence and returns a CReal
-- represented by that sequence which will memoize the values at each
-- precision. This is essential for getting good performance.
crMemoize :: (Int -> Integer) -> CReal n
crMemoize :: (Int -> Integer) -> CReal n
crMemoize Int -> Integer
fn = IO (CReal n) -> CReal n
forall a. IO a -> a
unsafePerformIO (IO (CReal n) -> CReal n) -> IO (CReal n) -> CReal n
forall a b. (a -> b) -> a -> b
$ do
  MVar Cache
mvc <- Cache -> IO (MVar Cache)
forall a. a -> IO (MVar a)
newMVar Cache
Never
  CReal n -> IO (CReal n)
forall (m :: * -> *) a. Monad m => a -> m a
return (CReal n -> IO (CReal n)) -> CReal n -> IO (CReal n)
forall a b. (a -> b) -> a -> b
$ MVar Cache -> (Int -> Integer) -> CReal n
forall (n :: Nat). MVar Cache -> (Int -> Integer) -> CReal n
CR MVar Cache
mvc Int -> Integer
fn

-- | crealPrecision x returns the type level parameter representing x's default
-- precision.
--
-- >>> crealPrecision (1 :: CReal 10)
-- 10
crealPrecision :: KnownNat n => CReal n -> Int
crealPrecision :: CReal n -> Int
crealPrecision = Integer -> Int
forall a. Num a => Integer -> a
fromInteger (Integer -> Int) -> (CReal n -> Integer) -> CReal n -> Int
forall b c a. (b -> c) -> (a -> b) -> a -> c
. CReal n -> Integer
forall (n :: Nat) (proxy :: Nat -> *).
KnownNat n =>
proxy n -> Integer
natVal

-- | @x \`atPrecision\` p@ returns the numerator of the pth element in the
-- Cauchy sequence represented by x. The denominator is 2^p.
--
-- >>> 10 `atPrecision` 10
-- 10240
atPrecision :: CReal n -> Int -> Integer
(CR MVar Cache
mvc Int -> Integer
f) atPrecision :: CReal n -> Int -> Integer
`atPrecision` (!Int
p) = IO Integer -> Integer
forall a. IO a -> a
unsafePerformIO (IO Integer -> Integer) -> IO Integer -> Integer
forall a b. (a -> b) -> a -> b
$ MVar Cache -> (Cache -> IO (Cache, Integer)) -> IO Integer
forall a b. MVar a -> (a -> IO (a, b)) -> IO b
modifyMVar MVar Cache
mvc ((Cache -> IO (Cache, Integer)) -> IO Integer)
-> (Cache -> IO (Cache, Integer)) -> IO Integer
forall a b. (a -> b) -> a -> b
$ \Cache
vc -> do
  Cache
vc' <- Cache -> IO Cache
forall a. a -> IO a
evaluate Cache
vc
  case Cache
vc' of
    Current Int
j Integer
v | Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
p ->
      (Cache, Integer) -> IO (Cache, Integer)
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Cache
vc', Integer
v Integer -> Int -> Integer
/^ (Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
p))
    Cache
_ -> do
      Integer
v <- Integer -> IO Integer
forall a. a -> IO a
evaluate (Integer -> IO Integer) -> Integer -> IO Integer
forall a b. (a -> b) -> a -> b
$ Int -> Integer
f Int
p
      let !vcn :: Cache
vcn = Int -> Integer -> Cache
Current Int
p Integer
v
      (Cache, Integer) -> IO (Cache, Integer)
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Cache
vcn, Integer
v)
{-# INLINABLE atPrecision #-}

-- | A CReal with precision p is shown as a decimal number d such that d is
-- within 2^-p of the true value.
--
-- >>> show (47176870 :: CReal 0)
-- "47176870"
--
-- >>> show (pi :: CReal 230)
-- "3.1415926535897932384626433832795028841971693993751058209749445923078164"
instance KnownNat n => Show (CReal n) where
  show :: CReal n -> String
show CReal n
x = Int -> CReal n -> String
forall (n :: Nat). Int -> CReal n -> String
showAtPrecision (CReal n -> Int
forall (n :: Nat). KnownNat n => CReal n -> Int
crealPrecision CReal n
x) CReal n
x

-- | The instance of Read will read an optionally signed number expressed in
-- decimal scientific notation
instance Read (CReal n) where
  readPrec :: ReadPrec (CReal n)
readPrec = ReadPrec (CReal n) -> ReadPrec (CReal n)
forall a. ReadPrec a -> ReadPrec a
parens (ReadPrec (CReal n) -> ReadPrec (CReal n))
-> ReadPrec (CReal n) -> ReadPrec (CReal n)
forall a b. (a -> b) -> a -> b
$ do
    Lexeme
lit <- ReadPrec Lexeme
lexP
    case Lexeme
lit of
      Number Number
n -> CReal n -> ReadPrec (CReal n)
forall (m :: * -> *) a. Monad m => a -> m a
return (CReal n -> ReadPrec (CReal n)) -> CReal n -> ReadPrec (CReal n)
forall a b. (a -> b) -> a -> b
$ Rational -> CReal n
forall a. Fractional a => Rational -> a
fromRational (Rational -> CReal n) -> Rational -> CReal n
forall a b. (a -> b) -> a -> b
$ Number -> Rational
L.numberToRational Number
n
      Symbol String
"-" -> Int -> ReadPrec (CReal n) -> ReadPrec (CReal n)
forall a. Int -> ReadPrec a -> ReadPrec a
prec Int
6 (ReadPrec (CReal n) -> ReadPrec (CReal n))
-> ReadPrec (CReal n) -> ReadPrec (CReal n)
forall a b. (a -> b) -> a -> b
$ do
        Lexeme
lit' <- ReadPrec Lexeme
lexP
        case Lexeme
lit' of
          Number Number
n -> CReal n -> ReadPrec (CReal n)
forall (m :: * -> *) a. Monad m => a -> m a
return (CReal n -> ReadPrec (CReal n)) -> CReal n -> ReadPrec (CReal n)
forall a b. (a -> b) -> a -> b
$ Rational -> CReal n
forall a. Fractional a => Rational -> a
fromRational (Rational -> CReal n) -> Rational -> CReal n
forall a b. (a -> b) -> a -> b
$ Rational -> Rational
forall a. Num a => a -> a
negate (Rational -> Rational) -> Rational -> Rational
forall a b. (a -> b) -> a -> b
$ Number -> Rational
L.numberToRational Number
n
          Lexeme
_ -> ReadPrec (CReal n)
forall a. ReadPrec a
pfail
      Lexeme
_ -> ReadPrec (CReal n)
forall a. ReadPrec a
pfail
  {-# INLINE readPrec #-}

  readListPrec :: ReadPrec [CReal n]
readListPrec = ReadPrec [CReal n]
forall a. Read a => ReadPrec [a]
readListPrecDefault
  {-# INLINE readListPrec #-}

  readsPrec :: Int -> ReadS (CReal n)
readsPrec = ReadPrec (CReal n) -> Int -> ReadS (CReal n)
forall a. ReadPrec a -> Int -> ReadS a
readPrec_to_S ReadPrec (CReal n)
forall a. Read a => ReadPrec a
readPrec
  {-# INLINE readsPrec #-}

  readList :: ReadS [CReal n]
readList = ReadPrec [CReal n] -> Int -> ReadS [CReal n]
forall a. ReadPrec a -> Int -> ReadS a
readPrec_to_S ReadPrec [CReal n]
forall a. Read a => ReadPrec [a]
readListPrec Int
0
  {-# INLINE readList #-}

-- | @signum (x :: CReal p)@ returns the sign of @x@ at precision @p@. It's
-- important to remember that this /may not/ represent the actual sign of @x@ if
-- the distance between @x@ and zero is less than 2^-@p@.
--
-- This is a little bit of a fudge, but it's probably better than failing to
-- terminate when trying to find the sign of zero. The class still respects the
-- abs-signum law though.
--
-- >>> signum (0.1 :: CReal 2)
-- 0.0
--
-- >>> signum (0.1 :: CReal 3)
-- 1.0
instance Num (CReal n) where
  {-# INLINE fromInteger #-}
  fromInteger :: Integer -> CReal n
fromInteger Integer
i = let
    !vc :: Cache
vc = Int -> Integer -> Cache
Current Int
0 Integer
i
    in IO (CReal n) -> CReal n
forall a. IO a -> a
unsafePerformIO (IO (CReal n) -> CReal n) -> IO (CReal n) -> CReal n
forall a b. (a -> b) -> a -> b
$ do
      MVar Cache
mvc <- Cache -> IO (MVar Cache)
forall a. a -> IO (MVar a)
newMVar Cache
vc
      CReal n -> IO (CReal n)
forall (m :: * -> *) a. Monad m => a -> m a
return (CReal n -> IO (CReal n)) -> CReal n -> IO (CReal n)
forall a b. (a -> b) -> a -> b
$ MVar Cache -> (Int -> Integer) -> CReal n
forall (n :: Nat). MVar Cache -> (Int -> Integer) -> CReal n
CR MVar Cache
mvc (Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
B.shiftL Integer
i)

  -- @negate@ and @abs@ try to give initial guesses, but don't wait if the
  -- @\'MVar\'@ is being used elsewhere.
  {-# INLINE negate #-}
  negate :: CReal n -> CReal n
negate (CR MVar Cache
mvc Int -> Integer
fn) = IO (CReal n) -> CReal n
forall a. IO a -> a
unsafePerformIO (IO (CReal n) -> CReal n) -> IO (CReal n) -> CReal n
forall a b. (a -> b) -> a -> b
$ do
    Maybe Cache
vcc <- MVar Cache -> IO (Maybe Cache)
forall a. MVar a -> IO (Maybe a)
tryReadMVar MVar Cache
mvc
    let
      !vcn :: Cache
vcn = case Maybe Cache
vcc of
        Maybe Cache
Nothing -> Cache
Never
        Just Cache
Never -> Cache
Never
        Just (Current Int
p Integer
v) -> Int -> Integer -> Cache
Current Int
p (Integer -> Integer
forall a. Num a => a -> a
negate Integer
v)
    MVar Cache
mvn <- Cache -> IO (MVar Cache)
forall a. a -> IO (MVar a)
newMVar Cache
vcn
    CReal n -> IO (CReal n)
forall (m :: * -> *) a. Monad m => a -> m a
return (CReal n -> IO (CReal n)) -> CReal n -> IO (CReal n)
forall a b. (a -> b) -> a -> b
$ MVar Cache -> (Int -> Integer) -> CReal n
forall (n :: Nat). MVar Cache -> (Int -> Integer) -> CReal n
CR MVar Cache
mvn (Integer -> Integer
forall a. Num a => a -> a
negate (Integer -> Integer) -> (Int -> Integer) -> Int -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Integer
fn)

  {-# INLINE abs #-}
  abs :: CReal n -> CReal n
abs (CR MVar Cache
mvc Int -> Integer
fn) = IO (CReal n) -> CReal n
forall a. IO a -> a
unsafePerformIO (IO (CReal n) -> CReal n) -> IO (CReal n) -> CReal n
forall a b. (a -> b) -> a -> b
$ do
    Maybe Cache
vcc <- MVar Cache -> IO (Maybe Cache)
forall a. MVar a -> IO (Maybe a)
tryReadMVar MVar Cache
mvc
    let
      !vcn :: Cache
vcn = case Maybe Cache
vcc of
        Maybe Cache
Nothing -> Cache
Never
        Just Cache
Never -> Cache
Never
        Just (Current Int
p Integer
v) -> Int -> Integer -> Cache
Current Int
p (Integer -> Integer
forall a. Num a => a -> a
abs Integer
v)
    MVar Cache
mvn <- Cache -> IO (MVar Cache)
forall a. a -> IO (MVar a)
newMVar Cache
vcn
    CReal n -> IO (CReal n)
forall (m :: * -> *) a. Monad m => a -> m a
return (CReal n -> IO (CReal n)) -> CReal n -> IO (CReal n)
forall a b. (a -> b) -> a -> b
$ MVar Cache -> (Int -> Integer) -> CReal n
forall (n :: Nat). MVar Cache -> (Int -> Integer) -> CReal n
CR MVar Cache
mvn (Integer -> Integer
forall a. Num a => a -> a
abs (Integer -> Integer) -> (Int -> Integer) -> Int -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Integer
fn)

  {-# INLINE (+) #-}
  CReal n
x1 + :: CReal n -> CReal n -> CReal n
+ CReal n
x2 = (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\Int
p -> let n1 :: Integer
n1 = CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x1 (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2)
                                 n2 :: Integer
n2 = CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x2 (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2)
                                 in (Integer
n1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
n2) Integer -> Int -> Integer
/^ Int
2)

  {-# INLINE (-) #-}
  CReal n
x1 - :: CReal n -> CReal n -> CReal n
- CReal n
x2 = (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\Int
p -> let n1 :: Integer
n1 = CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x1 (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2)
                                 n2 :: Integer
n2 = CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x2 (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2)
                                 in (Integer
n1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
n2) Integer -> Int -> Integer
/^ Int
2)

  {-# INLINE (*) #-}
  CReal n
x1 * :: CReal n -> CReal n -> CReal n
* CReal n
x2 = let
    s1 :: Int
s1 = Integer -> Int
log2 (Integer -> Integer
forall a. Num a => a -> a
abs (CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x1 Int
0) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
2) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
3
    s2 :: Int
s2 = Integer -> Int
log2 (Integer -> Integer
forall a. Num a => a -> a
abs (CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x2 Int
0) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
2) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
3
    in (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\Int
p -> let n1 :: Integer
n1 = CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x1 (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
s2)
                            n2 :: Integer
n2 = CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x2 (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
s1)
                        in (Integer
n1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
n2) Integer -> Int -> Integer
/^ (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
s1 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
s2))

  signum :: CReal n -> CReal n
signum CReal n
x = (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\Int
p -> Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
B.shiftL (Integer -> Integer
forall a. Num a => a -> a
signum (CReal n
x CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
`atPrecision` Int
p)) Int
p)

-- | Taking the reciprocal of zero will not terminate
instance Fractional (CReal n) where
  {-# INLINE fromRational #-}
  -- Use @roundD@ instead of @/.@ because we know @d > 0@ for a valid Rational.
  fromRational :: Rational -> CReal n
fromRational (Integer
n :% Integer
d) = (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\Int
p -> Integer -> Integer -> Integer
roundD (Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
B.shiftL Integer
n Int
p) Integer
d)

  {-# INLINE recip #-}
  -- TODO: Make recip 0 throw an error (if, for example, it would take more
  -- than 4GB of memory to represent the result)
  recip :: CReal n -> CReal n
recip CReal n
x = let
    s :: Int
s = (Int -> Bool) -> Int
findFirstMonotonic ((Integer
3 Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<=) (Integer -> Bool) -> (Int -> Integer) -> Int -> Bool
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Integer
forall a. Num a => a -> a
abs (Integer -> Integer) -> (Int -> Integer) -> Int -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x)
    in (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\Int
p -> let n :: Integer
n = CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
s Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2)
                        in Int -> Integer
forall a. Bits a => Int -> a
bit (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
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
s Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2) Integer -> Integer -> Integer
/. Integer
n)

instance Floating (CReal n) where
  -- TODO: Could we use something faster such as Ramanujan's formula
  pi :: CReal n
pi = CReal n
forall (n :: Nat). CReal n
piBy4 CReal n -> Int -> CReal n
forall (n :: Nat). CReal n -> Int -> CReal n
`shiftL` Int
2

  exp :: CReal n -> CReal n
exp CReal n
x = let o :: CReal n
o = CReal n -> Int -> CReal n
forall (n :: Nat). CReal n -> Int -> CReal n
shiftL (CReal n
x CReal n -> CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n -> CReal n
*. CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
recipBounded (CReal n -> Int -> CReal n
forall (n :: Nat). CReal n -> Int -> CReal n
shiftL CReal n
forall (n :: Nat). CReal n
ln2 Int
1)) Int
1
              l :: Integer
l = CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
o Int
0
              y :: CReal n
y = CReal n
x CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
- Integer -> CReal n
forall a. Num a => Integer -> a
fromInteger Integer
l CReal n -> CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n -> CReal n
*. CReal n
forall (n :: Nat). CReal n
ln2
          in if Integer
l Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0
               then CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
expBounded CReal n
x
               else CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
expBounded CReal n
y CReal n -> Int -> CReal n
forall (n :: Nat). CReal n -> Int -> CReal n
`shiftL` Integer -> Int
forall a. Num a => Integer -> a
fromInteger Integer
l

  -- | Range reduction on the principle that ln (a * b) = ln a + ln b
  log :: CReal n -> CReal n
log CReal n
x = let l :: Int
l = Integer -> Int
log2 (CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x Int
2) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
2
          in if   -- x <= 0.75
                | Int
l Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0  -> - CReal n -> CReal n
forall a. Floating a => a -> a
log (CReal n -> CReal n
forall a. Fractional a => a -> a
recip CReal n
x)
                  -- 0.75 <= x <= 2
                | Int
l Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0 -> CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
logBounded CReal n
x
                  -- x >= 2
                | Bool
otherwise  -> let a :: CReal n
a = CReal n
x CReal n -> Int -> CReal n
forall (n :: Nat). CReal n -> Int -> CReal n
`shiftR` Int
l
                                in CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
logBounded CReal n
a CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
+ Int -> CReal n
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
l CReal n -> CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n -> CReal n
*. CReal n
forall (n :: Nat). CReal n
ln2

  sqrt :: CReal n -> CReal n
sqrt CReal n
x = (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\Int
p -> let n :: Integer
n = CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
p)
                            in Integer -> Integer
isqrt Integer
n)

  -- | This will diverge when the base is not positive
  CReal n
x ** :: CReal n -> CReal n -> CReal n
** CReal n
y = CReal n -> CReal n
forall a. Floating a => a -> a
exp (CReal n -> CReal n
forall a. Floating a => a -> a
log CReal n
x CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
* CReal n
y)

  logBase :: CReal n -> CReal n -> CReal n
logBase CReal n
x CReal n
y = CReal n -> CReal n
forall a. Floating a => a -> a
log CReal n
y CReal n -> CReal n -> CReal n
forall a. Fractional a => a -> a -> a
/ CReal n -> CReal n
forall a. Floating a => a -> a
log CReal n
x

  sin :: CReal n -> CReal n
sin CReal n
x = CReal n -> CReal n
forall a. Floating a => a -> a
cos (CReal n
x CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
- CReal n
forall (n :: Nat). CReal n
piBy2)

  cos :: CReal n -> CReal n
cos CReal n
x = let o :: CReal n
o = CReal n -> Int -> CReal n
forall (n :: Nat). CReal n -> Int -> CReal n
shiftL (CReal n
x CReal n -> CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n -> CReal n
*. CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
recipBounded CReal n
forall a. Floating a => a
pi) Int
2
              s :: Integer
s = CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
o Int
1 Integer -> Int -> Integer
/^ Int
1
              octant :: Int
octant = Integer -> Int
forall a. Num a => Integer -> a
fromInteger (Integer -> Int) -> Integer -> Int
forall a b. (a -> b) -> a -> b
$ Integer
s Integer -> Integer -> Integer
forall a. Bits a => a -> a -> a
.&. Integer
7
              offset :: CReal n
offset = CReal n
x CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
- (Integer -> CReal n
forall a b. (Integral a, Num b) => a -> b
fromIntegral Integer
s CReal n -> CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n -> CReal n
*. CReal n
forall (n :: Nat). CReal n
piBy4)
              fs :: [CReal n -> CReal n]
fs = [          CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
cosBounded
                   , CReal n -> CReal n
forall a. Num a => a -> a
negate (CReal n -> CReal n) -> (CReal n -> CReal n) -> CReal n -> CReal n
forall b c a. (b -> c) -> (a -> b) -> a -> c
. CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
sinBounded (CReal n -> CReal n) -> (CReal n -> CReal n) -> CReal n -> CReal n
forall b c a. (b -> c) -> (a -> b) -> a -> c
. CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
subtract CReal n
forall (n :: Nat). CReal n
piBy4
                   , CReal n -> CReal n
forall a. Num a => a -> a
negate (CReal n -> CReal n) -> (CReal n -> CReal n) -> CReal n -> CReal n
forall b c a. (b -> c) -> (a -> b) -> a -> c
. CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
sinBounded
                   , CReal n -> CReal n
forall a. Num a => a -> a
negate (CReal n -> CReal n) -> (CReal n -> CReal n) -> CReal n -> CReal n
forall b c a. (b -> c) -> (a -> b) -> a -> c
. CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
cosBounded (CReal n -> CReal n) -> (CReal n -> CReal n) -> CReal n -> CReal n
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (CReal n
forall (n :: Nat). CReal n
piBy4CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
-)
                   , CReal n -> CReal n
forall a. Num a => a -> a
negate (CReal n -> CReal n) -> (CReal n -> CReal n) -> CReal n -> CReal n
forall b c a. (b -> c) -> (a -> b) -> a -> c
. CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
cosBounded
                   ,          CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
sinBounded (CReal n -> CReal n) -> (CReal n -> CReal n) -> CReal n -> CReal n
forall b c a. (b -> c) -> (a -> b) -> a -> c
. CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
subtract CReal n
forall (n :: Nat). CReal n
piBy4
                   ,          CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
sinBounded
                   ,          CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
cosBounded (CReal n -> CReal n) -> (CReal n -> CReal n) -> CReal n -> CReal n
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (CReal n
forall (n :: Nat). CReal n
piBy4CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
-)]
          in ([CReal n -> CReal n]
forall (n :: Nat). [CReal n -> CReal n]
fs [CReal n -> CReal n] -> Int -> CReal n -> CReal n
forall a. [a] -> Int -> a
!! Int
octant) CReal n
offset

  tan :: CReal n -> CReal n
tan CReal n
x = CReal n -> CReal n
forall a. Floating a => a -> a
sin CReal n
x CReal n -> CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n -> CReal n
.* CReal n -> CReal n
forall a. Fractional a => a -> a
recip (CReal n -> CReal n
forall a. Floating a => a -> a
cos CReal n
x)

  asin :: CReal n -> CReal n
asin CReal n
x = CReal n -> CReal n
forall a. Floating a => a -> a
atan (CReal n
x CReal n -> CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n -> CReal n
.*. CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
recipBounded (CReal n
1 CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
+ CReal n -> CReal n
forall a. Floating a => a -> a
sqrt (CReal n
1 CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
- CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
squareBounded CReal n
x))) CReal n -> Int -> CReal n
forall (n :: Nat). CReal n -> Int -> CReal n
`shiftL` Int
1

  acos :: CReal n -> CReal n
acos CReal n
x = CReal n
forall (n :: Nat). CReal n
piBy2 CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
- CReal n -> CReal n
forall a. Floating a => a -> a
asin CReal n
x

  atan :: CReal n -> CReal n
atan CReal n
x = let -- q is 4 times x to within 1/4
               q :: Integer
q = CReal n
x CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
`atPrecision` Int
2
           in if   -- x <= -1
                 | Integer
q Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<  -Integer
4 -> CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
atanBounded (CReal n -> CReal n
forall a. Num a => a -> a
negate (CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
recipBounded CReal n
x)) CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
- CReal n
forall (n :: Nat). CReal n
piBy2
                   -- -1.25 <= x <= -0.75
                 | Integer
q Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== -Integer
4 -> -(CReal n
forall (n :: Nat). CReal n
piBy4 CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
+ CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
atanBounded ((CReal n
x CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
+ CReal n
1) CReal n -> CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n -> CReal n
.*. CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
recipBounded (CReal n
x CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
- CReal n
1)))
                   -- 0.75 <= x <= 1.25
                 | Integer
q Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
==  Integer
4 -> CReal n
forall (n :: Nat). CReal n
piBy4 CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
+ CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
atanBounded ((CReal n
x CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
- CReal n
1) CReal n -> CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n -> CReal n
.*. CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
recipBounded (CReal n
x CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
+ CReal n
1))
                   -- x >= 1
                 | Integer
q Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>   Integer
4 -> CReal n
forall (n :: Nat). CReal n
piBy2 CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
- CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
atanBounded (CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
recipBounded CReal n
x)
                   -- -0.75 <= x <= 0.75
                 | Bool
otherwise -> CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
atanBounded CReal n
x

  -- TODO: benchmark replacing these with their series expansion
  sinh :: CReal n -> CReal n
sinh CReal n
x = let (CReal n
expX, CReal n
expNegX) = CReal n -> (CReal n, CReal n)
forall (n :: Nat). CReal n -> (CReal n, CReal n)
expPosNeg CReal n
x
           in (CReal n
expX CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
- CReal n
expNegX) CReal n -> Int -> CReal n
forall (n :: Nat). CReal n -> Int -> CReal n
`shiftR` Int
1
  cosh :: CReal n -> CReal n
cosh CReal n
x = let (CReal n
expX, CReal n
expNegX) = CReal n -> (CReal n, CReal n)
forall (n :: Nat). CReal n -> (CReal n, CReal n)
expPosNeg CReal n
x
           in (CReal n
expX CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
+ CReal n
expNegX) CReal n -> Int -> CReal n
forall (n :: Nat). CReal n -> Int -> CReal n
`shiftR` Int
1
  tanh :: CReal n -> CReal n
tanh CReal n
x = let e2x :: CReal n
e2x = CReal n -> CReal n
forall a. Floating a => a -> a
exp (CReal n
x CReal n -> Int -> CReal n
forall (n :: Nat). CReal n -> Int -> CReal n
`shiftL` Int
1)
           in (CReal n
e2x CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
- CReal n
1) CReal n -> CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n -> CReal n
*. CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
recipBounded (CReal n
e2x CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
+ CReal n
1)

  asinh :: CReal n -> CReal n
asinh CReal n
x = CReal n -> CReal n
forall a. Floating a => a -> a
log (CReal n
x CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
+ CReal n -> CReal n
forall a. Floating a => a -> a
sqrt (CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
square CReal n
x CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
+ CReal n
1))
  acosh :: CReal n -> CReal n
acosh CReal n
x = CReal n -> CReal n
forall a. Floating a => a -> a
log (CReal n
x CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
+ CReal n -> CReal n
forall a. Floating a => a -> a
sqrt (CReal n
x CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
+ CReal n
1) CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
* CReal n -> CReal n
forall a. Floating a => a -> a
sqrt (CReal n
x CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
- CReal n
1))
  atanh :: CReal n -> CReal n
atanh CReal n
x = (CReal n -> CReal n
forall a. Floating a => a -> a
log (CReal n
1 CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
+ CReal n
x) CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
- CReal n -> CReal n
forall a. Floating a => a -> a
log (CReal n
1 CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
- CReal n
x)) CReal n -> Int -> CReal n
forall (n :: Nat). CReal n -> Int -> CReal n
`shiftR` Int
1

-- | 'toRational' returns the CReal n evaluated at a precision of 2^-n
instance KnownNat n => Real (CReal n) where
  toRational :: CReal n -> Rational
toRational CReal n
x = let p :: Int
p = CReal n -> Int
forall (n :: Nat). KnownNat n => CReal n -> Int
crealPrecision CReal n
x
                 in CReal n
x CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
`atPrecision` Int
p Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% Int -> Integer
forall a. Bits a => Int -> a
bit Int
p

instance KnownNat n => RealFrac (CReal n) where
  properFraction :: CReal n -> (b, CReal n)
properFraction CReal n
x = let p :: Int
p = CReal n -> Int
forall (n :: Nat). KnownNat n => CReal n -> Int
crealPrecision CReal n
x
                         v :: Integer
v = CReal n
x CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
`atPrecision` Int
p
                         r :: Integer
r = Integer
v Integer -> Integer -> Integer
forall a. Bits a => a -> a -> a
.&. (Int -> Integer
forall a. Bits a => Int -> a
bit Int
p Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
1)
                         c :: Integer
c = Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftR (Integer
v Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
r) Int
p
                         n :: Integer
n = if Integer
c Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 Bool -> Bool -> Bool
&& Integer
r Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= Integer
0 then Integer
c Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1 else Integer
c
                         f :: CReal n
f = CReal n -> Integer -> CReal n
forall (n :: Nat). CReal n -> Integer -> CReal n
plusInteger CReal n
x (Integer -> Integer
forall a. Num a => a -> a
negate Integer
n)
                     in (Integer -> b
forall a. Num a => Integer -> a
fromInteger Integer
n, CReal n
f)

  truncate :: CReal n -> b
truncate CReal n
x = let p :: Int
p = CReal n -> Int
forall (n :: Nat). KnownNat n => CReal n -> Int
crealPrecision CReal n
x
                   v :: Integer
v = CReal n
x CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
`atPrecision` Int
p
                   r :: Integer
r = Integer
v Integer -> Integer -> Integer
forall a. Bits a => a -> a -> a
.&. (Int -> Integer
forall a. Bits a => Int -> a
bit Int
p Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
1)
                   c :: Integer
c = Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftR (Integer
v Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
r) Int
p
                   n :: Integer
n = if Integer
c Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 Bool -> Bool -> Bool
&& Integer
r Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= Integer
0 then Integer
c Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1 else Integer
c
                   in Integer -> b
forall a. Num a => Integer -> a
fromInteger Integer
n

  round :: CReal n -> b
round CReal n
x = let p :: Int
p = CReal n -> Int
forall (n :: Nat). KnownNat n => CReal n -> Int
crealPrecision CReal n
x
                n :: Integer
n = (CReal n
x CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
`atPrecision` Int
p) Integer -> Int -> Integer
/^ Int
p
                in Integer -> b
forall a. Num a => Integer -> a
fromInteger Integer
n

  ceiling :: CReal n -> b
ceiling CReal n
x = let p :: Int
p = CReal n -> Int
forall (n :: Nat). KnownNat n => CReal n -> Int
crealPrecision CReal n
x
                  v :: Integer
v = CReal n
x CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
`atPrecision` Int
p
                  r :: Integer
r = Integer
v Integer -> Integer -> Integer
forall a. Bits a => a -> a -> a
.&. (Int -> Integer
forall a. Bits a => Int -> a
bit Int
p Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
1)
                  n :: Integer
n = Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftR (Integer
v Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
r) Int
p
                  in if Integer
r Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= Integer
0 then Integer -> b
forall a. Num a => Integer -> a
fromInteger (Integer -> b) -> Integer -> b
forall a b. (a -> b) -> a -> b
$ Integer
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1 else Integer -> b
forall a. Num a => Integer -> a
fromInteger Integer
n

  floor :: CReal n -> b
floor CReal n
x = let p :: Int
p = CReal n -> Int
forall (n :: Nat). KnownNat n => CReal n -> Int
crealPrecision CReal n
x
                v :: Integer
v = CReal n
x CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
`atPrecision` Int
p
                r :: Integer
r = Integer
v Integer -> Integer -> Integer
forall a. Bits a => a -> a -> a
.&. (Int -> Integer
forall a. Bits a => Int -> a
bit Int
p Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
1)
                n :: Integer
n = Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftR (Integer
v Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
r) Int
p
                in Integer -> b
forall a. Num a => Integer -> a
fromInteger Integer
n

-- | Several of the functions in this class ('floatDigits', 'floatRange',
-- 'exponent', 'significand') only make sense for floats represented by a
-- mantissa and exponent. These are bound to error.
--
-- @atan2 y x `atPrecision` p@ performs the comparison to determine the
-- quadrant at precision p. This can cause atan2 to be slightly slower than atan
instance KnownNat n => RealFloat (CReal n) where
  floatRadix :: CReal n -> Integer
floatRadix CReal n
_ = Integer
2
  floatDigits :: CReal n -> Int
floatDigits CReal n
_ = String -> Int
forall a. HasCallStack => String -> a
error String
"Data.CReal.Internal floatDigits"
  floatRange :: CReal n -> (Int, Int)
floatRange CReal n
_ = String -> (Int, Int)
forall a. HasCallStack => String -> a
error String
"Data.CReal.Internal floatRange"
  decodeFloat :: CReal n -> (Integer, Int)
decodeFloat CReal n
x = let p :: Int
p = CReal n -> Int
forall (n :: Nat). KnownNat n => CReal n -> Int
crealPrecision CReal n
x
                  in (CReal n
x CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
`atPrecision` Int
p, -Int
p)
  encodeFloat :: Integer -> Int -> CReal n
encodeFloat Integer
m Int
n = if Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
0
    then Rational -> CReal n
forall a. Fractional a => Rational -> a
fromRational (Integer
m Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% Int -> Integer
forall a. Bits a => Int -> a
bit (Int -> Int
forall a. Num a => a -> a
negate Int
n))
    else Rational -> CReal n
forall a. Fractional a => Rational -> a
fromRational (Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
m Int
n Integer -> Integer -> Rational
forall a. a -> a -> Ratio a
:% Integer
1)
  exponent :: CReal n -> Int
exponent = String -> CReal n -> Int
forall a. HasCallStack => String -> a
error String
"Data.CReal.Internal exponent"
  significand :: CReal n -> CReal n
significand = String -> CReal n -> CReal n
forall a. HasCallStack => String -> a
error String
"Data.CReal.Internal significand"
  scaleFloat :: Int -> CReal n -> CReal n
scaleFloat = (CReal n -> Int -> CReal n) -> Int -> CReal n -> CReal n
forall a b c. (a -> b -> c) -> b -> a -> c
flip CReal n -> Int -> CReal n
forall (n :: Nat). CReal n -> Int -> CReal n
shiftL
  isNaN :: CReal n -> Bool
isNaN CReal n
_ = Bool
False
  isInfinite :: CReal n -> Bool
isInfinite CReal n
_ = Bool
False
  isDenormalized :: CReal n -> Bool
isDenormalized CReal n
_ = Bool
False
  isNegativeZero :: CReal n -> Bool
isNegativeZero CReal n
_ = Bool
False
  isIEEE :: CReal n -> Bool
isIEEE CReal n
_ = Bool
False
  atan2 :: CReal n -> CReal n -> CReal n
atan2 CReal n
y CReal n
x = (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\Int
p ->
    let y' :: Integer
y' = CReal n
y CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
`atPrecision` Int
p
        x' :: Integer
x' = CReal n
x CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
`atPrecision` Int
p
        θ :: CReal n
θ = if | Integer
x' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
0            ->  CReal n -> CReal n
forall a. Floating a => a -> a
atan (CReal n
yCReal n -> CReal n -> CReal n
forall a. Fractional a => a -> a -> a
/CReal n
x)
               | Integer
x' Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 Bool -> Bool -> Bool
&& Integer
y' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
0 ->  CReal n
forall (n :: Nat). CReal n
piBy2
               | Integer
x' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<  Integer
0 Bool -> Bool -> Bool
&& Integer
y' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
0 ->  CReal n
forall a. Floating a => a
pi CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
+ CReal n -> CReal n
forall a. Floating a => a -> a
atan (CReal n
yCReal n -> CReal n -> CReal n
forall a. Fractional a => a -> a -> a
/CReal n
x)
               | Integer
x' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
0 Bool -> Bool -> Bool
&& Integer
y' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 -> -CReal n -> CReal n -> CReal n
forall a. RealFloat a => a -> a -> a
atan2 (-CReal n
y) CReal n
x
               | Integer
y' Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 Bool -> Bool -> Bool
&& Integer
x' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 ->  CReal n
forall a. Floating a => a
pi    -- must be after the previous test on zero y
               | Integer
x'Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
==Integer
0 Bool -> Bool -> Bool
&& Integer
y'Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
==Integer
0    ->  CReal n
0     -- must be after the other double zero tests
               | Bool
otherwise         ->  String -> CReal n
forall a. HasCallStack => String -> a
error String
"Data.CReal.Internal atan2"
    in CReal n
θ CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
`atPrecision` Int
p)

-- | Values of type @CReal p@ are compared for equality at precision @p@. This
-- may cause values which differ by less than 2^-p to compare as equal.
--
-- >>> 0 == (0.1 :: CReal 1)
-- True
instance KnownNat n => Eq (CReal n) where
  -- TODO, should this try smaller values first?
  CR MVar Cache
mvx Int -> Integer
_ == :: CReal n -> CReal n -> Bool
== CR MVar Cache
mvy Int -> Integer
_ | MVar Cache
mvx MVar Cache -> MVar Cache -> Bool
forall a. Eq a => a -> a -> Bool
== MVar Cache
mvy = Bool
True
  CReal n
x == CReal n
y = let p :: Int
p = CReal n -> Int
forall (n :: Nat). KnownNat n => CReal n -> Int
crealPrecision CReal n
x
           in ((CReal n
x CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
- CReal n
y) CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
`atPrecision` Int
p) Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0


-- | Like equality, values of type @CReal p@ are compared at precision @p@.
instance KnownNat n => Ord (CReal n) where
  compare :: CReal n -> CReal n -> Ordering
compare (CR MVar Cache
mvx Int -> Integer
_) (CR MVar Cache
mvy Int -> Integer
_) | MVar Cache
mvx MVar Cache -> MVar Cache -> Bool
forall a. Eq a => a -> a -> Bool
== MVar Cache
mvy = Ordering
EQ
  compare CReal n
x CReal n
y =
    let p :: Int
p = CReal n -> Int
forall (n :: Nat). KnownNat n => CReal n -> Int
crealPrecision CReal n
x
    in Integer -> Integer -> Ordering
forall a. Ord a => a -> a -> Ordering
compare ((CReal n
x CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
- CReal n
y) CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
`atPrecision` Int
p) Integer
0
  max :: CReal n -> CReal n -> CReal n
max CReal n
x CReal n
y = (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\Int
p -> Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
max (CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x Int
p) (CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
y Int
p))
  min :: CReal n -> CReal n -> CReal n
min CReal n
x CReal n
y = (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\Int
p -> Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
min (CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x Int
p) (CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
y Int
p))

-- | The 'Random' instance for @\'CReal\' p@ will return random number with at
-- least @p@ digits of precision, every digit after that is zero.
instance KnownNat n => Random (CReal n) where
  randomR :: (CReal n, CReal n) -> g -> (CReal n, g)
randomR (CReal n
lo, CReal n
hi) g
g = let d :: CReal n
d = CReal n
hi CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
- CReal n
lo
                           l :: Int
l = Int
1 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Integer -> Int
log2 (CReal n -> CReal n
forall a. Num a => a -> a
abs CReal n
d CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
`atPrecision` Int
0)
                           p :: Int
p = Int
l Int -> Int -> Int
forall a. Num a => a -> a -> a
+ CReal n -> Int
forall (n :: Nat). KnownNat n => CReal n -> Int
crealPrecision CReal n
lo
                           (Integer
n, g
g') = (Integer, Integer) -> g -> (Integer, g)
forall a g. (Random a, RandomGen g) => (a, a) -> g -> (a, g)
randomR (Integer
0, Integer
2Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^Int
p) g
g
                           r :: CReal n
r = Rational -> CReal n
forall a. Fractional a => Rational -> a
fromRational (Integer
n Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% Integer
2Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^Int
p)
                       in (CReal n
r CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
* CReal n
d CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
+ CReal n
lo, g
g')
  random :: g -> (CReal n, g)
random g
g = let p :: Int
p = Int
1 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ CReal n -> Int
forall (n :: Nat). KnownNat n => CReal n -> Int
crealPrecision (CReal n
forall a. HasCallStack => a
undefined :: CReal n)
                 (Integer
n, g
g') = (Integer, Integer) -> g -> (Integer, g)
forall a g. (Random a, RandomGen g) => (a, a) -> g -> (a, g)
randomR (Integer
0, Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
max Integer
0 (Integer
2Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^Int
p Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
2)) g
g
                 r :: CReal n
r = Rational -> CReal n
forall a. Fractional a => Rational -> a
fromRational (Integer
n Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% Integer
2Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^Int
p)
             in (CReal n
r, g
g')

--------------------------------------------------------------------------------
-- Some utility functions
--------------------------------------------------------------------------------

--
-- Constants
--

piBy4 :: CReal n
piBy4 :: CReal n
piBy4 = CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
atanBounded (CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
recipBounded CReal n
5) CReal n -> Int -> CReal n
forall (n :: Nat). CReal n -> Int -> CReal n
`shiftL` Int
2 CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
- CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
atanBounded (CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
recipBounded CReal n
239) -- Machin Formula

piBy2 :: CReal n
piBy2 :: CReal n
piBy2 = CReal n
forall (n :: Nat). CReal n
piBy4 CReal n -> Int -> CReal n
forall (n :: Nat). CReal n -> Int -> CReal n
`shiftL` Int
1

ln2 :: CReal n
ln2 :: CReal n
ln2 = CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
logBounded CReal n
2

--
-- Bounded multiplication
--

infixl 7 `mulBounded`, `mulBoundedL`, .*, *., .*.

-- | Alias for @'mulBoundedL'@
(.*) :: CReal n -> CReal n -> CReal n
.* :: CReal n -> CReal n -> CReal n
(.*) = CReal n -> CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n -> CReal n
mulBoundedL

-- | Alias for @flip 'mulBoundedL'@
(*.) :: CReal n -> CReal n -> CReal n
*. :: CReal n -> CReal n -> CReal n
(*.) = (CReal n -> CReal n -> CReal n) -> CReal n -> CReal n -> CReal n
forall a b c. (a -> b -> c) -> b -> a -> c
flip CReal n -> CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n -> CReal n
mulBoundedL

-- | Alias for @'mulBounded'@
(.*.) :: CReal n -> CReal n -> CReal n
.*. :: CReal n -> CReal n -> CReal n
(.*.) = CReal n -> CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n -> CReal n
mulBounded

-- | A more efficient multiply with the restriction that the first argument
-- must be in the closed range [-1..1]
mulBoundedL :: CReal n -> CReal n -> CReal n
mulBoundedL :: CReal n -> CReal n -> CReal n
mulBoundedL CReal n
x1 CReal n
x2 = let
  s1 :: Int
s1 = Int
4
  s2 :: Int
s2 = Integer -> Int
log2 (Integer -> Integer
forall a. Num a => a -> a
abs (CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x2 Int
0) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
2) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
3
  in (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\Int
p -> let n1 :: Integer
n1 = CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x1 (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
s2)
                          n2 :: Integer
n2 = CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x2 (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
s1)
                      in (Integer
n1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
n2) Integer -> Int -> Integer
/^ (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
s1 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
s2))

-- | A more efficient multiply with the restriction that both values must be
-- in the closed range [-1..1]
mulBounded :: CReal n -> CReal n -> CReal n
mulBounded :: CReal n -> CReal n -> CReal n
mulBounded CReal n
x1 CReal n
x2 = let
  s1 :: Int
s1 = Int
4
  s2 :: Int
s2 = Int
4
  in (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\Int
p -> let n1 :: Integer
n1 = CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x1 (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
s2)
                          n2 :: Integer
n2 = CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x2 (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
s1)
                      in (Integer
n1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
n2) Integer -> Int -> Integer
/^ (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
s1 Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
s2))

-- | A more efficient 'recip' with the restriction that the input must have
-- absolute value greater than or equal to 1
recipBounded :: CReal n -> CReal n
recipBounded :: CReal n -> CReal n
recipBounded CReal n
x = (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\Int
p -> let s :: Int
s = Int
2
                                      n :: Integer
n = CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
s Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2)
                                  in Int -> Integer
forall a. Bits a => Int -> a
bit (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
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
s Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2) Integer -> Integer -> Integer
/. Integer
n)

-- | Return the square of the input, more efficient than @('*')@
{-# INLINABLE square #-}
square :: CReal n -> CReal n
square :: CReal n -> CReal n
square CReal n
x = let
  s :: Int
s = Integer -> Int
log2 (Integer -> Integer
forall a. Num a => a -> a
abs (CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x Int
0) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
2) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
3
  in (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\Int
p -> let n :: Integer
n = CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
s)
                      in (Integer
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
n) Integer -> Int -> Integer
/^ (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
s))

-- | A more efficient 'square' with the restrictuion that the value must be in
-- the closed range [-1..1]
{-# INLINABLE squareBounded #-}
squareBounded :: CReal n -> CReal n
squareBounded :: CReal n -> CReal n
squareBounded x :: CReal n
x@(CR MVar Cache
mvc Int -> Integer
_) = IO (CReal n) -> CReal n
forall a. IO a -> a
unsafePerformIO (IO (CReal n) -> CReal n) -> IO (CReal n) -> CReal n
forall a b. (a -> b) -> a -> b
$ do
  Maybe Cache
vcc <- MVar Cache -> IO (Maybe Cache)
forall a. MVar a -> IO (Maybe a)
tryReadMVar MVar Cache
mvc
  let
    !s :: Int
s = Int
4
    !vcn :: Cache
vcn = case Maybe Cache
vcc of
      Maybe Cache
Nothing -> Cache
Never
      Just Cache
Never -> Cache
Never
      Just (Current Int
j Integer
n) -> case Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
s of
        Int
p | Int
p Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 -> Cache
Never
        Int
p -> Int -> Integer -> Cache
Current Int
p ((Integer
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
n) Integer -> Int -> Integer
/^ (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
s))
    fn' :: Int -> Integer
fn' !Int
p = let n :: Integer
n = CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
s)
             in (Integer
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
n) Integer -> Int -> Integer
/^ (Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
s)
  MVar Cache
mvn <- Cache -> IO (MVar Cache)
forall a. a -> IO (MVar a)
newMVar Cache
vcn
  CReal n -> IO (CReal n)
forall (m :: * -> *) a. Monad m => a -> m a
return (CReal n -> IO (CReal n)) -> CReal n -> IO (CReal n)
forall a b. (a -> b) -> a -> b
$ MVar Cache -> (Int -> Integer) -> CReal n
forall (n :: Nat). MVar Cache -> (Int -> Integer) -> CReal n
CR MVar Cache
mvn Int -> Integer
fn'

--
-- Bounded exponential functions and expPosNeg
--

-- | A more efficient 'exp' with the restriction that the input must be in the
-- closed range [-1..1]
expBounded :: CReal n -> CReal n
expBounded :: CReal n -> CReal n
expBounded CReal n
x = let q :: [Rational]
q = (Integer
1Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
%) (Integer -> Rational) -> [Integer] -> [Rational]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> (Integer -> Integer -> Integer)
-> Integer -> [Integer] -> [Integer]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl' Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
(*) Integer
1 [Integer
1..]
               in [Rational] -> (Int -> Int) -> CReal n -> CReal n
forall (n :: Nat). [Rational] -> (Int -> Int) -> CReal n -> CReal n
powerSeries [Rational]
q (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
5) CReal n
x

-- | A more efficient 'log' with the restriction that the input must be in the
-- closed range [2/3..2]
logBounded :: CReal n -> CReal n
logBounded :: CReal n -> CReal n
logBounded CReal n
x = let q :: [Rational]
q = [Integer
1 Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% Integer
n | Integer
n <- [Integer
1..]]
                   y :: CReal n
y = (CReal n
x CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
- CReal n
1) CReal n -> CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n -> CReal n
.* CReal n -> CReal n
forall a. Fractional a => a -> a
recip CReal n
x
               in CReal n
y CReal n -> CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n -> CReal n
.* [Rational] -> (Int -> Int) -> CReal n -> CReal n
forall (n :: Nat). [Rational] -> (Int -> Int) -> CReal n -> CReal n
powerSeries [Rational]
q Int -> Int
forall a. a -> a
id CReal n
y

-- | @expPosNeg x@ returns @(exp x, exp (-x))#
expPosNeg :: CReal n -> (CReal n, CReal n)
expPosNeg :: CReal n -> (CReal n, CReal n)
expPosNeg CReal n
x = let o :: CReal n
o = CReal n
x CReal n -> CReal n -> CReal n
forall a. Fractional a => a -> a -> a
/ CReal n
forall (n :: Nat). CReal n
ln2
                  l :: Integer
l = CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
o Int
0
                  y :: CReal n
y = CReal n
x CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
- Integer -> CReal n
forall a. Num a => Integer -> a
fromInteger Integer
l CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
* CReal n
forall (n :: Nat). CReal n
ln2
              in if Integer
l Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0
                   then (CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
expBounded CReal n
x, CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
expBounded (-CReal n
x))
                   else (CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
expBounded CReal n
y CReal n -> Int -> CReal n
forall (n :: Nat). CReal n -> Int -> CReal n
`shiftL` Integer -> Int
forall a. Num a => Integer -> a
fromInteger Integer
l,
                         CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
expBounded (CReal n -> CReal n
forall a. Num a => a -> a
negate CReal n
y) CReal n -> Int -> CReal n
forall (n :: Nat). CReal n -> Int -> CReal n
`shiftR` Integer -> Int
forall a. Num a => Integer -> a
fromInteger Integer
l)

--
-- Bounded trigonometric functions
--

-- | A more efficient 'sin' with the restriction that the input must be in the
-- closed range [-1..1]
sinBounded :: CReal n -> CReal n
sinBounded :: CReal n -> CReal n
sinBounded CReal n
x = let q :: [Rational]
q = [Rational] -> [Rational]
forall a. Num a => [a] -> [a]
alternateSign ((Rational -> Rational -> Rational)
-> Rational -> [Rational] -> [Rational]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl' Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
(*) Rational
1 [ Integer
1 Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% (Integer
nInteger -> 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
2,Integer
4..]])
               in CReal n
x CReal n -> CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n -> CReal n
.* [Rational] -> (Int -> Int) -> CReal n -> CReal n
forall (n :: Nat). [Rational] -> (Int -> Int) -> CReal n -> CReal n
powerSeries [Rational]
q (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
1) (CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
squareBounded CReal n
x)

-- | A more efficient 'cos' with the restriction that the input must be in the
-- closed range [-1..1]
cosBounded :: CReal n -> CReal n
cosBounded :: CReal n -> CReal n
cosBounded CReal n
x = let q :: [Rational]
q = [Rational] -> [Rational]
forall a. Num a => [a] -> [a]
alternateSign ((Rational -> Rational -> Rational)
-> Rational -> [Rational] -> [Rational]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl' Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
(*) Rational
1 [Integer
1 Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% (Integer
nInteger -> 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
3..]])
               in [Rational] -> (Int -> Int) -> CReal n -> CReal n
forall (n :: Nat). [Rational] -> (Int -> Int) -> CReal n -> CReal n
powerSeries [Rational]
q (Int -> Int -> Int
forall a. Ord a => a -> a -> a
max Int
1) (CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
squareBounded CReal n
x)

-- | A more efficient 'atan' with the restriction that the input must be in the
-- closed range [-1..1]
atanBounded :: CReal n -> CReal n
atanBounded :: CReal n -> CReal n
atanBounded CReal n
x = let q :: [Rational]
q = (Rational -> Rational -> Rational)
-> Rational -> [Rational] -> [Rational]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl' Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
(*) Rational
1 [Integer
n Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% (Integer
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1) | Integer
n <- [Integer
2,Integer
4..]]
                    s :: CReal n
s = CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
squareBounded CReal n
x
                    rd :: CReal n
rd = CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
recipBounded (CReal n -> Integer -> CReal n
forall (n :: Nat). CReal n -> Integer -> CReal n
plusInteger CReal n
s Integer
1)
                in (CReal n
x CReal n -> CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n -> CReal n
.*. CReal n
rd) CReal n -> CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n -> CReal n
.* [Rational] -> (Int -> Int) -> CReal n -> CReal n
forall (n :: Nat). [Rational] -> (Int -> Int) -> CReal n -> CReal n
powerSeries [Rational]
q (Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) (CReal n
s CReal n -> CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n -> CReal n
.*. CReal n
rd)

--
-- Integer addition
--

infixl 6 `plusInteger`

-- | @x \`plusInteger\` n@ is equal to @x + fromInteger n@, but more efficient
{-# INLINE plusInteger #-}
plusInteger :: CReal n -> Integer -> CReal n
plusInteger :: CReal n -> Integer -> CReal n
plusInteger CReal n
x Integer
0 = CReal n
x
plusInteger (CR MVar Cache
mvc Int -> Integer
fn) Integer
n = IO (CReal n) -> CReal n
forall a. IO a -> a
unsafePerformIO (IO (CReal n) -> CReal n) -> IO (CReal n) -> CReal n
forall a b. (a -> b) -> a -> b
$ do
  Maybe Cache
vcc <- MVar Cache -> IO (Maybe Cache)
forall a. MVar a -> IO (Maybe a)
tryReadMVar MVar Cache
mvc
  let
    !vcn :: Cache
vcn = case Maybe Cache
vcc of
      Maybe Cache
Nothing -> Cache
Never
      Just Cache
Never -> Cache
Never
      Just (Current Int
j Integer
v) -> Int -> Integer -> Cache
Current Int
j (Integer
v Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
n Int
j)
    fn' :: Int -> Integer
fn' !Int
p = Int -> Integer
fn Int
p Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
B.shiftL Integer
n Int
p
  MVar Cache
mvc' <- Cache -> IO (MVar Cache)
forall a. a -> IO (MVar a)
newMVar Cache
vcn
  CReal n -> IO (CReal n)
forall (m :: * -> *) a. Monad m => a -> m a
return (CReal n -> IO (CReal n)) -> CReal n -> IO (CReal n)
forall a b. (a -> b) -> a -> b
$ MVar Cache -> (Int -> Integer) -> CReal n
forall (n :: Nat). MVar Cache -> (Int -> Integer) -> CReal n
CR MVar Cache
mvc' Int -> Integer
fn'

--
-- Multiplication with powers of two
--

infixl 8 `shiftL`, `shiftR`

-- | @x \`shiftR\` n@ is equal to @x@ divided by 2^@n@
--
-- @n@ can be negative or zero
--
-- This can be faster than doing the division
shiftR :: CReal n -> Int -> CReal n
shiftR :: CReal n -> Int -> CReal n
shiftR CReal n
x Int
n = (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\Int
p -> let p' :: Int
p' = Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
n
                              in if Int
p' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
0
                                 then CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x Int
p'
                                 else CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x Int
0 Integer -> Int -> Integer
/^ Int -> Int
forall a. Num a => a -> a
negate Int
p')

-- | @x \`shiftL\` n@ is equal to @x@ multiplied by 2^@n@
--
-- @n@ can be negative or zero
--
-- This can be faster than doing the multiplication
shiftL :: CReal n -> Int -> CReal n
shiftL :: CReal n -> Int -> CReal n
shiftL CReal n
x = CReal n -> Int -> CReal n
forall (n :: Nat). CReal n -> Int -> CReal n
shiftR CReal n
x (Int -> CReal n) -> (Int -> Int) -> Int -> CReal n
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Int
forall a. Num a => a -> a
negate


--
-- Showing CReals
--

-- | Return a string representing a decimal number within 2^-p of the value
-- represented by the given @CReal p@.
showAtPrecision :: Int -> CReal n -> String
showAtPrecision :: Int -> CReal n -> String
showAtPrecision Int
p CReal n
x = let places :: Int
places = Int -> Int
decimalDigitsAtPrecision Int
p
                          r :: Rational
r      = CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x Int
p Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% Int -> Integer
forall a. Bits a => Int -> a
bit Int
p
                      in Int -> Rational -> String
rationalToDecimal Int
places Rational
r

-- | How many decimal digits are required to represent a number to within 2^-p
decimalDigitsAtPrecision :: Int -> Int
decimalDigitsAtPrecision :: Int -> Int
decimalDigitsAtPrecision Int
0 = Int
0
decimalDigitsAtPrecision Int
p = Integer -> Int
log10 (Int -> Integer
forall a. Bits a => Int -> a
bit Int
p) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1

-- | @rationalToDecimal p x@ returns a string representing @x@ at @p@ decimal
-- places.
rationalToDecimal :: Int -> Rational -> String
rationalToDecimal :: Int -> Rational -> String
rationalToDecimal Int
places (Integer
n :% Integer
d) = String
p String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
is String -> ShowS
forall a. [a] -> [a] -> [a]
++ if Int
places Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
0 then String
"." String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
fs else String
""
  where p :: String
p = case Integer -> Integer
forall a. Num a => a -> a
signum Integer
n of
              -1 -> String
"-"
              Integer
_  -> String
""
        ds :: String
ds = Integer -> String
forall a. Show a => a -> String
show (Integer -> Integer -> Integer
roundD (Integer -> Integer
forall a. Num a => a -> a
abs Integer
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
10Integer -> Int -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^Int
places) Integer
d)
        l :: Int
l = String -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length String
ds
        (String
is, String
fs) = if Int
l Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
places then (String
"0", Int -> Char -> String
forall a. Int -> a -> [a]
replicate (Int
places Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
l) Char
'0' String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
ds) else Int -> String -> (String, String)
forall a. Int -> [a] -> ([a], [a])
splitAt (Int
l Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
places) String
ds


--
-- Integer operations
--

divZeroErr :: a
divZeroErr :: a
divZeroErr = String -> a
forall a. HasCallStack => String -> a
error String
"Division by zero"
{-# NOINLINE divZeroErr #-}

roundD :: Integer -> Integer -> Integer
roundD :: Integer -> Integer -> Integer
roundD Integer
n Integer
d = case Integer -> Integer -> (Integer, Integer)
forall a. Integral a => a -> a -> (a, a)
divMod Integer
n Integer
d of
  (Integer
q, Integer
r) -> case Integer -> Integer -> Ordering
forall a. Ord a => a -> a -> Ordering
compare (Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
r Int
1) Integer
d of
    Ordering
LT -> Integer
q
    Ordering
EQ -> if Integer -> Int -> Bool
forall a. Bits a => a -> Int -> Bool
testBit Integer
q Int
0 then Integer
q Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1 else Integer
q
    Ordering
GT -> Integer
q Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1
{-# INLINE roundD #-}

infixl 7 /.
-- | Division rounding to the nearest integer and rounding half integers to the
-- nearest even integer.
(/.) :: Integer -> Integer -> Integer
(!Integer
n) /. :: Integer -> Integer -> Integer
/. (!Integer
d) = case Integer -> Integer -> Ordering
forall a. Ord a => a -> a -> Ordering
compare Integer
d Integer
0 of
  Ordering
LT -> Integer -> Integer -> Integer
roundD (Integer -> Integer
forall a. Num a => a -> a
negate Integer
n) (Integer -> Integer
forall a. Num a => a -> a
negate Integer
d)
  Ordering
EQ -> Integer
forall a. a
divZeroErr
  Ordering
GT -> Integer -> Integer -> Integer
roundD Integer
n Integer
d
{-# INLINABLE (/.) #-}

infixl 7 /^
-- | @n /^ p@ is equivalent to @n \'/.\' (2^p)@, but faster, and it works for
-- negative values of p.
(/^) :: Integer -> Int -> Integer
(!Integer
n) /^ :: Integer -> Int -> Integer
/^ (!Int
p) = case Int -> Int -> Ordering
forall a. Ord a => a -> a -> Ordering
compare Int
p Int
0 of
  Ordering
LT -> Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
n (Int -> Int
forall a. Num a => a -> a
negate Int
p)
  Ordering
EQ -> Integer
n
  Ordering
GT -> let
    !bp :: Integer
bp = Int -> Integer
forall a. Bits a => Int -> a
bit Int
p
    !r :: Integer
r = Integer
n Integer -> Integer -> Integer
forall a. Bits a => a -> a -> a
.&. (Integer
bp Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
1)
    !q :: Integer
q = Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftR (Integer
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
r) Int
p
    in case Integer -> Integer -> Ordering
forall a. Ord a => a -> a -> Ordering
compare (Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
r Int
1) Integer
bp of
      Ordering
LT -> Integer
q
      Ordering
EQ -> if Integer -> Int -> Bool
forall a. Bits a => a -> Int -> Bool
testBit Integer
q Int
0 then Integer
q Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1 else Integer
q
      Ordering
GT -> Integer
q Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1

-- | @log2 x@ returns the base 2 logarithm of @x@ rounded towards zero.
--
-- The input must be positive
{-# INLINE log2 #-}
log2 :: Integer -> Int
log2 :: Integer -> Int
log2 Integer
x = Int# -> Int
I# (Integer -> Int#
integerLog2# Integer
x)

-- | @log10 x@ returns the base 10 logarithm of @x@ rounded towards zero.
--
-- The input must be positive
{-# INLINE log10 #-}
log10 :: Integer -> Int
log10 :: Integer -> Int
log10 Integer
x = Int# -> Int
I# (Integer -> Integer -> Int#
integerLogBase# Integer
10 Integer
x)

-- | @isqrt x@ returns the square root of @x@ rounded towards zero.
--
-- The input must not be negative
{-# INLINABLE isqrt #-}
isqrt :: Integer -> Integer
isqrt :: Integer -> Integer
isqrt Integer
x | Integer
x Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0     = String -> Integer
forall a. HasCallStack => String -> a
error String
"Sqrt applied to negative Integer"
        | Integer
x Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0    = Integer
0
        | Bool
otherwise = (Integer -> Bool) -> (Integer -> Integer) -> Integer -> Integer
forall a. (a -> Bool) -> (a -> a) -> a -> a
until Integer -> Bool
satisfied Integer -> Integer
improve Integer
initialGuess
  where improve :: Integer -> Integer
improve Integer
r    = Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftR (Integer
r Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
x Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`div` Integer
r) Int
1
        satisfied :: Integer -> Bool
satisfied Integer
r  = let r2 :: Integer
r2 = Integer
r Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
r in Integer
r2 Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
x Bool -> Bool -> Bool
&& Integer
r2 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer -> Int -> Integer
forall a. Bits a => a -> Int -> a
unsafeShiftL Integer
r Int
1 Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>= Integer
x
        initialGuess :: Integer
initialGuess = Int -> Integer
forall a. Bits a => Int -> a
bit (Int -> Int -> Int
forall a. Bits a => a -> Int -> a
unsafeShiftR (Integer -> Int
log2 Integer
x) Int
1)

--
-- Searching
--

-- | Given a monotonic function
{-# INLINABLE findFirstMonotonic #-}
findFirstMonotonic :: (Int -> Bool) -> Int
findFirstMonotonic :: (Int -> Bool) -> Int
findFirstMonotonic Int -> Bool
p = Int -> Int -> Int
findBounds Int
0 Int
1
  where findBounds :: Int -> Int -> Int
findBounds !Int
l !Int
u = if Int -> Bool
p Int
u then Int -> Int -> Int
binarySearch Int
l Int
u
                                  else Int -> Int -> Int
findBounds Int
u (Int
u Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
2)
        binarySearch :: Int -> Int -> Int
binarySearch !Int
l !Int
u = let !m :: Int
m = Int
l Int -> Int -> Int
forall a. Num a => a -> a -> a
+ (Int
u Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
l) Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2
                             in if | Int
lInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1 Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
u  -> Int
l
                                   | Int -> Bool
p Int
m       -> Int -> Int -> Int
binarySearch Int
l Int
m
                                   | Bool
otherwise -> Int -> Int -> Int
binarySearch Int
m Int
u


--
-- Power series
--

-- | Apply 'negate' to every other element, starting with the second
--
-- >>> alternateSign [1..5]
-- [1,-2,3,-4,5]
{-# INLINABLE alternateSign #-}
alternateSign :: Num a => [a] -> [a]
alternateSign :: [a] -> [a]
alternateSign [a]
ls = (a -> (Bool -> [a]) -> Bool -> [a])
-> (Bool -> [a]) -> [a] -> Bool -> [a]
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr
  (\a
a Bool -> [a]
r Bool
b -> if Bool
b then a -> a
forall a. Num a => a -> a
negate a
a a -> [a] -> [a]
forall a. a -> [a] -> [a]
: Bool -> [a]
r Bool
False else a
a a -> [a] -> [a]
forall a. a -> [a] -> [a]
: Bool -> [a]
r Bool
True)
  ([a] -> Bool -> [a]
forall a b. a -> b -> a
const [])
  [a]
ls
  Bool
False

-- | @powerSeries q f x `atPrecision` p@ will evaluate the power series with
-- coefficients @q@ up to the coefficient at index @f p@ at value @x@
--
-- @f@ should be a function such that the CReal invariant is maintained. This
-- means that if the power series @y = a[0] + a[1] + a[2] + ...@ is evaluated
-- at precision @p@ then the sum of every @a[n]@ for @n > f p@ must be less than
-- 2^-p.
--
-- This is used by all the bounded transcendental functions.
--
-- >>> let (!) x = product [2..x]
-- >>> powerSeries [1 % (n!) | n <- [0..]] (max 5) 1 :: CReal 218
-- 2.718281828459045235360287471352662497757247093699959574966967627724
powerSeries :: [Rational] -> (Int -> Int) -> CReal n -> CReal n
powerSeries :: [Rational] -> (Int -> Int) -> CReal n -> CReal n
powerSeries [Rational]
q Int -> Int
termsAtPrecision CReal n
x = (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize
     (\Int
p -> let t :: Int
t = Int -> Int
termsAtPrecision Int
p
                d :: Int
d = Integer -> Int
log2 (Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
t) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
2
                p' :: Int
p' = Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
d
                p'' :: Int
p'' = Int
p' Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
d
                m :: Integer
m = CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x Int
p''
                xs :: [Rational]
xs = (Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% Integer
1) (Integer -> Rational) -> [Integer] -> [Rational]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> (Integer -> Integer) -> Integer -> [Integer]
forall a. (a -> a) -> a -> [a]
iterate (\Integer
e -> Integer
m Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
e Integer -> Int -> Integer
/^ Int
p'') (Int -> Integer
forall a. Bits a => Int -> a
bit Int
p')
                r :: Integer
r = [Integer] -> Integer
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ([Integer] -> Integer)
-> ([Rational] -> [Integer]) -> [Rational] -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> [Integer] -> [Integer]
forall a. Int -> [a] -> [a]
take (Int
t Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) ([Integer] -> [Integer])
-> ([Rational] -> [Integer]) -> [Rational] -> [Integer]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Rational -> Integer) -> [Rational] -> [Integer]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap (Rational -> Integer
forall a b. (RealFrac a, Integral b) => a -> b
round (Rational -> Integer)
-> (Rational -> Rational) -> Rational -> Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
* Integer -> Rational
forall a. Num a => Integer -> a
fromInteger (Int -> Integer
forall a. Bits a => Int -> a
bit Int
d))) ([Rational] -> Integer) -> [Rational] -> Integer
forall a b. (a -> b) -> a -> b
$ (Rational -> Rational -> Rational)
-> [Rational] -> [Rational] -> [Rational]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith Rational -> Rational -> Rational
forall a. Num a => a -> a -> a
(*) [Rational]
q [Rational]
xs
            in Integer
r Integer -> Int -> Integer
/^ (Int
2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d))