{-# 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(..), RandomGen(..))
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 fn :: 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 mvc :: MVar Cache
mvc f :: 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
$ \vc :: Cache
vc -> do
  Cache
vc' <- Cache -> IO Cache
forall a. a -> IO a
evaluate Cache
vc
  case Cache
vc' of
    Current j :: Int
j v :: Integer
v | Int
j Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
p -> do
      (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))
    _ -> 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 x :: 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 n :: 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 "-" -> Int -> ReadPrec (CReal n) -> ReadPrec (CReal n)
forall a. Int -> ReadPrec a -> ReadPrec a
prec 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 n :: 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
          _ -> ReadPrec (CReal n)
forall a. ReadPrec a
pfail
      _ -> 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 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 i :: Integer
i = let
    !vc :: Cache
vc = Int -> Integer -> Cache
Current 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 mvc :: MVar Cache
mvc fn :: 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
        Nothing -> Cache
Never
        Just Never -> Cache
Never
        Just (Current p :: Int
p v :: 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 mvc :: MVar Cache
mvc fn :: 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
        Nothing -> Cache
Never
        Just Never -> Cache
Never
        Just (Current p :: Int
p v :: 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 (+) #-}
  x1 :: CReal n
x1 + :: CReal n -> CReal n -> CReal n
+ x2 :: CReal n
x2 = (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\p :: 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
+ 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
+ 2)
                                 in (Integer
n1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
n2) Integer -> Int -> Integer
/^ 2)

  {-# INLINE (-) #-}
  x1 :: CReal n
x1 - :: CReal n -> CReal n -> CReal n
- x2 :: CReal n
x2 = (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\p :: 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
+ 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
+ 2)
                                 in (Integer
n1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
n2) Integer -> Int -> Integer
/^ 2)

  {-# INLINE (*) #-}
  x1 :: CReal n
x1 * :: CReal n -> CReal n -> CReal n
* x2 :: 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 0) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ 2) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ 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 0) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ 2) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ 3
    in (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\p :: 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 x :: CReal n
x = (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\p :: 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 (n :: Integer
n :% d :: Integer
d) = (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\p :: 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 x :: CReal n
x = let
    s :: Int
s = (Int -> Bool) -> Int
findFirstMonotonic ((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 (\p :: 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
+ 2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
s Int -> Int -> Int
forall a. Num a => a -> a -> a
+ 2)
                        in Int -> Integer
forall a. Bits a => Int -> a
bit (2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ 2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
s Int -> Int -> Int
forall a. Num a => a -> a -> a
+ 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` 2

  exp :: CReal n -> CReal n
exp x :: 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 1)) 1
              l :: Integer
l = CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
o 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
== 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 x :: 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 2) Int -> Int -> Int
forall a. Num a => a -> a -> a
- 2
          in if   -- x <= 0.75
                | Int
l Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< 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
== 0 -> CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
logBounded CReal n
x
                  -- x >= 2
                | Int
l Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> 0  -> 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 x :: CReal n
x = (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\p :: Int
p -> let n :: Integer
n = CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x (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
  x :: CReal n
x ** :: CReal n -> CReal n -> CReal n
** y :: 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 x :: CReal n
x y :: 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 x :: 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 x :: 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) 2
              s :: Integer
s = CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
o 1 Integer -> Int -> Integer
/^ 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
.&. 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 x :: 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 x :: 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 (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 (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` 1

  acos :: CReal n -> CReal n
acos x :: 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 x :: 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` 2
           in if   -- x <= -1
                 | Integer
q Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<  -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
== -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
+ 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
- 1)))
                   -- 0.75 <= x <= 1.25
                 | Integer
q Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
==  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
- 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
+ 1))
                   -- x >= 1
                 | Integer
q Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
>   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 x :: CReal n
x = let (expX :: CReal n
expX, expNegX :: 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` 1
  cosh :: CReal n -> CReal n
cosh x :: CReal n
x = let (expX :: CReal n
expX, expNegX :: 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` 1
  tanh :: CReal n -> CReal n
tanh x :: 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` 1)
           in (CReal n
e2x CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
- 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
+ 1)

  asinh :: CReal n -> CReal n
asinh x :: 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
+ 1))
  acosh :: CReal n -> CReal n
acosh x :: 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
+ 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
- 1))
  atanh :: CReal n -> CReal n
atanh x :: CReal n
x = (CReal n -> CReal n
forall a. Floating a => a -> a
log (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 (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` 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 x :: 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 x :: 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
- 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
< 0 Bool -> Bool -> Bool
&& Integer
r Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= 0 then Integer
c Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ 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 x :: 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
- 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
< 0 Bool -> Bool -> Bool
&& Integer
r Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= 0 then Integer
c Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ 1 else Integer
c
                   in Integer -> b
forall a. Num a => Integer -> a
fromInteger Integer
n

  round :: CReal n -> b
round x :: 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 x :: 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
- 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 case Integer
r Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= 0 of
                    True -> 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
+ 1
                    _    -> Integer -> b
forall a. Num a => Integer -> a
fromInteger Integer
n

  floor :: CReal n -> b
floor x :: 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
- 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 _ = 2
  floatDigits :: CReal n -> Int
floatDigits _ = String -> Int
forall a. HasCallStack => String -> a
error "Data.CReal.Internal floatDigits"
  floatRange :: CReal n -> (Int, Int)
floatRange _ = String -> (Int, Int)
forall a. HasCallStack => String -> a
error "Data.CReal.Internal floatRange"
  decodeFloat :: CReal n -> (Integer, Int)
decodeFloat x :: 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 m :: Integer
m n :: Int
n = if Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= 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
:% 1)
  exponent :: CReal n -> Int
exponent = String -> CReal n -> Int
forall a. HasCallStack => String -> a
error "Data.CReal.Internal exponent"
  significand :: CReal n -> CReal n
significand = String -> CReal n -> CReal n
forall a. HasCallStack => String -> a
error "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 _ = Bool
False
  isInfinite :: CReal n -> Bool
isInfinite _ = Bool
False
  isDenormalized :: CReal n -> Bool
isDenormalized _ = Bool
False
  isNegativeZero :: CReal n -> Bool
isNegativeZero _ = Bool
False
  isIEEE :: CReal n -> Bool
isIEEE _ = Bool
False
  atan2 :: CReal n -> CReal n -> CReal n
atan2 y :: CReal n
y x :: CReal n
x = (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\p :: 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
> 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
== 0 Bool -> Bool -> Bool
&& Integer
y' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> 0 ->  CReal n
forall (n :: Nat). CReal n
piBy2
               | Integer
x' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<  0 Bool -> Bool -> Bool
&& Integer
y' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> 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
<= 0 Bool -> Bool -> Bool
&& Integer
y' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< 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
== 0 Bool -> Bool -> Bool
&& Integer
x' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< 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
==0 Bool -> Bool -> Bool
&& Integer
y'Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
==0    ->  0     -- must be after the other double zero tests
               | Bool
otherwise         ->  String -> CReal n
forall a. HasCallStack => String -> a
error "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 mvx :: MVar Cache
mvx _ == :: CReal n -> CReal n -> Bool
== CR mvy :: MVar Cache
mvy _ | MVar Cache
mvx MVar Cache -> MVar Cache -> Bool
forall a. Eq a => a -> a -> Bool
== MVar Cache
mvy = Bool
True
  x :: CReal n
x == y :: CReal n
y = let p :: Int
p = CReal n -> Int
forall (n :: Nat). KnownNat n => CReal n -> Int
crealPrecision CReal n
x Int -> Int -> Int
forall a. Num a => a -> a -> a
+ 2
           in (CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x Int
p Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
y Int
p) Integer -> Int -> Integer
/^ 2 Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 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 mvx :: MVar Cache
mvx _) (CR mvy :: MVar Cache
mvy _) | MVar Cache
mvx MVar Cache -> MVar Cache -> Bool
forall a. Eq a => a -> a -> Bool
== MVar Cache
mvy = Ordering
EQ
  compare x :: CReal n
x y :: CReal n
y = let p :: Int
p = CReal n -> Int
forall (n :: Nat). KnownNat n => CReal n -> Int
crealPrecision CReal n
x Int -> Int -> Int
forall a. Num a => a -> a -> a
+ 2
                in Integer -> Integer -> Ordering
forall a. Ord a => a -> a -> Ordering
compare ((CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
x Int
p Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- CReal n -> Int -> Integer
forall (n :: Nat). CReal n -> Int -> Integer
atPrecision CReal n
y Int
p) Integer -> Int -> Integer
/^ 2) 0
  max :: CReal n -> CReal n -> CReal n
max x :: CReal n
x y :: CReal n
y = (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\p :: 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 x :: CReal n
x y :: CReal n
y = (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\p :: 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 (lo :: CReal n
lo, hi :: CReal n
hi) g :: 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 = 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` 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
                           (n :: Integer
n, g' :: g
g') = (Integer, Integer) -> g -> (Integer, g)
forall a g. (Random a, RandomGen g) => (a, a) -> g -> (a, g)
randomR (0, 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
% 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
g = let p :: Int
p = 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)
                 (n :: Integer
n, g' :: g
g') = (Integer, Integer) -> g -> (Integer, g)
forall a g. (Random a, RandomGen g) => (a, a) -> g -> (a, g)
randomR (0, Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
max 0 (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
- 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
% 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 5) CReal n -> Int -> CReal n
forall (n :: Nat). CReal n -> Int -> CReal n
`shiftL` 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 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` 1

ln2 :: CReal n
ln2 :: CReal n
ln2 = CReal n -> CReal n
forall (n :: Nat). CReal n -> CReal n
logBounded 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 x1 :: CReal n
x1 x2 :: CReal n
x2 = let
  s1 :: Int
s1 = 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 0) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ 2) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ 3
  in (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\p :: 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 x1 :: CReal n
x1 x2 :: CReal n
x2 = let
  s1 :: Int
s1 = 4
  s2 :: Int
s2 = 4
  in (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\p :: 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 x :: CReal n
x = (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\p :: Int
p -> let s :: Int
s = 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
+ 2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
s Int -> Int -> Int
forall a. Num a => a -> a -> a
+ 2)
                                  in Int -> Integer
forall a. Bits a => Int -> a
bit (2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
p Int -> Int -> Int
forall a. Num a => a -> a -> a
+ 2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
s Int -> Int -> Int
forall a. Num a => a -> a -> a
+ 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 x :: 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 0) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ 2) Int -> Int -> Int
forall a. Num a => a -> a -> a
+ 3
  in (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\p :: 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
+ 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 mvc :: MVar Cache
mvc _) = 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 = 4
    !vcn :: Cache
vcn = case Maybe Cache
vcc of
      Nothing -> Cache
Never
      Just Never -> Cache
Never
      Just (Current j :: Int
j n :: Integer
n) -> case Int
j Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
s of
        p :: Int
p | Int
p Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< 0 -> Cache
Never
        p :: 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
+ 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
+ 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 x :: CReal n
x = let q :: [Rational]
q = (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
(*) 1 [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 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 x :: CReal n
x = let q :: [Rational]
q = [1 Integer -> Integer -> Rational
forall a. Integral a => a -> a -> Ratio a
% Integer
n | Integer
n <- [1..]]
                   y :: CReal n
y = (CReal n
x CReal n -> CReal n -> CReal n
forall a. Num a => a -> a -> a
- 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 x :: 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 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
== 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 x :: 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
(*) 1 [ 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
+1)) | Integer
n <- [2,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 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 x :: 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
(*) 1 [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
+1)) | Integer
n <- [1,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 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 x :: 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
(*) 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
+ 1) | Integer
n <- [2,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 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
+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 x :: CReal n
x 0 = CReal n
x
plusInteger (CR mvc :: MVar Cache
mvc fn :: Int -> Integer
fn) n :: 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
      Nothing -> Cache
Never
      Just Never -> Cache
Never
      Just (Current j :: Int
j v :: 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 x :: CReal n
x n :: Int
n = (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize (\p :: 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
>= 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 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 x :: 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 p :: Int
p x :: 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 0 = 0
decimalDigitsAtPrecision p :: 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
+ 1

-- | @rationalToDecimal p x@ returns a string representing @x@ at @p@ decimal
-- places.
rationalToDecimal :: Int -> Rational -> String
rationalToDecimal :: Int -> Rational -> String
rationalToDecimal places :: Int
places (n :: Integer
n :% d :: 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
> 0 then "." String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
fs else ""
  where p :: String
p = case Integer -> Integer
forall a. Num a => a -> a
signum Integer
n of
              -1 -> "-"
              _  -> ""
        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
* 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
        (is :: String
is, fs :: String
fs) = if | Int
l Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
places -> ("0", Int -> Char -> String
forall a. Int -> a -> [a]
replicate (Int
places Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
l) '0' String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
ds)
                      | Bool
otherwise -> 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 "Division by zero"
{-# NOINLINE divZeroErr #-}

roundD :: Integer -> Integer -> Integer
roundD :: Integer -> Integer -> Integer
roundD n :: Integer
n d :: Integer
d = case Integer -> Integer -> (Integer, Integer)
forall a. Integral a => a -> a -> (a, a)
divMod Integer
n Integer
d of
  (q :: Integer
q, r :: 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 1) Integer
d of
    LT -> Integer
q
    EQ -> if Integer -> Int -> Bool
forall a. Bits a => a -> Int -> Bool
testBit Integer
q 0 then Integer
q Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ 1 else Integer
q
    GT -> Integer
q Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ 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 0 of
  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)
  EQ -> Integer
forall a. a
divZeroErr
  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 0 of
  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)
  EQ -> Integer
n
  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
- 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 1) Integer
bp of
      LT -> Integer
q
      EQ -> if Integer -> Int -> Bool
forall a. Bits a => a -> Int -> Bool
testBit Integer
q 0 then Integer
q Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ 1 else Integer
q
      GT -> Integer
q Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ 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 x :: 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 x :: Integer
x = Int# -> Int
I# (Integer -> Integer -> Int#
integerLogBase# 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 x :: Integer
x | Integer
x Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< 0     = String -> Integer
forall a. HasCallStack => String -> a
error "Sqrt applied to negative Integer"
        | Integer
x Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== 0    = 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 r :: 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)) 1
        satisfied :: Integer -> Bool
satisfied r :: 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 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) 1)

--
-- Searching
--

-- | Given a monotonic function
{-# INLINABLE findFirstMonotonic #-}
findFirstMonotonic :: (Int -> Bool) -> Int
findFirstMonotonic :: (Int -> Bool) -> Int
findFirstMonotonic p :: Int -> Bool
p = Int -> Int -> Int
findBounds 0 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
* 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` 2)
                             in if | Int
lInt -> Int -> Int
forall a. Num a => a -> a -> a
+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 = \ls :: [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
a r :: Bool -> [a]
r b :: 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
aa -> [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 q :: [Rational]
q termsAtPrecision :: Int -> Int
termsAtPrecision x :: CReal n
x = (Int -> Integer) -> CReal n
forall (n :: Nat). (Int -> Integer) -> CReal n
crMemoize
     (\p :: 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
+ 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
%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 (\e :: 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
+ 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
/^ (2 Int -> Int -> Int
forall a. Num a => a -> a -> a
* Int
d))