exact-real 0.12.3 → 0.12.4

raw patch · 4 files changed

+323/−122 lines, 4 filesdep −memoizePVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies removed: memoize

API changes (from Hackage documentation)

Files

• .gitignore +2/−0
• changelog.md +5/−0
• exact-real.cabal +2/−3
• src/Data/CReal/Internal.hs +314/−119

@@ -14,3 +14,5 @@ *.aux *.hp .stack-work/+result+.hdevtools

@@ -2,5 +2,10 @@  ## WIP +## [0.12.4] - 2020-06-07+  - Big speedup (orbits testsuite about 9 times faster)++Big thanks for @Zemyla for the new memoization scheme+ ## [0.12.3] - 2020-05-29   - More relaxed version bounds

@@ -4,10 +4,10 @@ -- -- see: https://github.com/sol/hpack ----- hash: 4207194af32e5bff84d5b39c43a8685b42b56f671421be69c390c71e933be2b5+-- hash: 3b8926fece094589ffa737495e0936d123b8c2d9a4e8341be340b08916967c77  name:           exact-real-version:        0.12.3+version:        0.12.4 synopsis:       Exact real arithmetic description:    A type to represent exact real numbers using fast binary Cauchy sequences. category:       Math@@ -49,7 +49,6 @@   build-depends:       base >=4.8 && <5     , integer-gmp-    , memoize >=0.7     , random >=1.0   default-language: Haskell2010 

@@ -1,9 +1,11 @@ {-# 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@@ -14,12 +16,16 @@   (     -- * 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   , (.*.)@@ -30,6 +36,7 @@   , shiftL   , shiftR   , square+  , squareBounded      -- ** Exponential   , expBounded@@ -48,6 +55,7 @@      -- ** Integer operations   , (/.)+  , (/^)   , log2   , log10   , isqrt@@ -59,15 +67,20 @@   ) where  import Data.List (scanl')-import Data.Ratio (numerator,denominator,(%))+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 Numeric (readSigned, readFloat)-import Data.Function.Memoize (memoize)-import System.Random (Random(..))+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" #-}+{-# ANN module ("HLint: ignore Reduce duplication" :: String) #-}  -- $setup -- >>> :set -XDataKinds@@ -75,17 +88,28 @@  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 (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.-newtype CReal (n :: Nat) = CR (Int -> Integer)+-- 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 = CR . memoize+crMemoize fn = unsafePerformIO $ do+  mvc <- newMVar Never+  return $ CR mvc fn  -- | crealPrecision x returns the type level parameter representing x's default -- precision.@@ -101,7 +125,16 @@ -- >>> 10 `atPrecision` 10 -- 10240 atPrecision :: CReal n -> Int -> Integer-(CR x) `atPrecision` p = x p+(CR mvc f) `atPrecision` (!p) = unsafePerformIO $ modifyMVar mvc $ \vc -> do+  vc' <- evaluate vc+  case vc' of+    Current j v | j >= p -> do+      pure (vc', v /^ (j - p))+    _ -> do+      v <- evaluate $ f p+      let !vcn = Current p v+      pure (vcn, 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.@@ -116,9 +149,28 @@  -- | The instance of Read will read an optionally signed number expressed in -- decimal scientific notation-instance KnownNat n => Read (CReal n) where-  readsPrec _ = readSigned readFloat+instance Read (CReal n) where+  readPrec = parens $ do+    lit <- lexP+    case lit of+      Number n -> return $ fromRational $ L.numberToRational n+      Symbol "-" -> prec 6 $ do+        lit' <- lexP+        case lit' of+          Number n -> return $ fromRational $ negate $ L.numberToRational n+          _ -> pfail+      _ -> pfail+  {-# INLINE readPrec #-} +  readListPrec = readListPrecDefault+  {-# INLINE readListPrec #-}++  readsPrec = readPrec_to_S readPrec+  {-# INLINE readsPrec #-}++  readList = readPrec_to_S 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@.@@ -134,53 +186,83 @@ -- 1.0 instance Num (CReal n) where   {-# INLINE fromInteger #-}-  fromInteger i = crMemoize (\p -> i * 2 ^ p)+  fromInteger i = let+    !vc = Current 0 i+    in unsafePerformIO $ do+      mvc <- newMVar vc+      return $ CR mvc (B.shiftL i) +  -- @negate@ and @abs@ try to give initial guesses, but don't wait if the+  -- @\'MVar\'@ is being used elsewhere.   {-# INLINE negate #-}-  negate (CR x) = crMemoize (negate . x)+  negate (CR mvc fn) = unsafePerformIO $ do+    vcc <- tryReadMVar mvc+    let+      !vcn = case vcc of+        Nothing -> Never+        Just Never -> Never+        Just (Current p v) -> Current p (negate v)+    mvn <- newMVar vcn+    return $ CR mvn (negate . fn)    {-# INLINE abs #-}-  abs (CR x) = crMemoize (abs . x)+  abs (CR mvc fn) = unsafePerformIO $ do+    vcc <- tryReadMVar mvc+    let+      !vcn = case vcc of+        Nothing -> Never+        Just Never -> Never+        Just (Current p v) -> Current p (abs v)+    mvn <- newMVar vcn+    return $ CR mvn (abs . fn)    {-# INLINE (+) #-}-  CR x1 + CR x2 = crMemoize (\p -> let n1 = x1 (p + 2)-                                       n2 = x2 (p + 2)-                                   in (n1 + n2) /. 4)+  x1 + x2 = crMemoize (\p -> let n1 = atPrecision x1 (p + 2)+                                 n2 = atPrecision x2 (p + 2)+                                 in (n1 + n2) /^ 2) +  {-# INLINE (-) #-}+  x1 - x2 = crMemoize (\p -> let n1 = atPrecision x1 (p + 2)+                                 n2 = atPrecision x2 (p + 2)+                                 in (n1 - n2) /^ 2)+   {-# INLINE (*) #-}-  CR x1 * CR x2 = crMemoize (\p -> let s1 = log2 (abs (x1 0) + 2) + 3-                                       s2 = log2 (abs (x2 0) + 2) + 3-                                       n1 = x1 (p + s2)-                                       n2 = x2 (p + s1)-                                   in (n1 * n2) /. 2^(p + s1 + s2)  )+  x1 * x2 = let+    s1 = log2 (abs (atPrecision x1 0) + 2) + 3+    s2 = log2 (abs (atPrecision x2 0) + 2) + 3+    in crMemoize (\p -> let n1 = atPrecision x1 (p + s2)+                            n2 = atPrecision x2 (p + s1)+                        in (n1 * n2) /^ (p + s1 + s2)) -  signum x = crMemoize (\p -> signum (x `atPrecision` p) * 2^p)+  signum x = crMemoize (\p -> B.shiftL (signum (x `atPrecision` p)) p)  -- | Taking the reciprocal of zero will not terminate instance Fractional (CReal n) where-  fromRational n = fromInteger (numerator n) *. recipBounded (fromInteger (denominator n))+  {-# INLINE fromRational #-}+  -- Use @roundD@ instead of @/.@ because we know @d > 0@ for a valid Rational.+  fromRational (n :% d) = crMemoize (\p -> roundD (B.shiftL n p) 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 (CR x) = crMemoize (\p -> let s = findFirstMonotonic ((3 <=) . abs . x)-                                      n = x (p + 2 * s + 2)-                                  in 2^(2 * p + 2 * s + 2) /. n)+  recip x = let+    s = findFirstMonotonic ((3 <=) . abs . atPrecision x)+    in crMemoize (\p -> let n = atPrecision x (p + 2 * s + 2)+                        in bit (2 * p + 2 * s + 2) /. n)  instance Floating (CReal n) where   -- TODO: Could we use something faster such as Ramanujan's formula-  pi = 4 * piBy4+  pi = piBy4 `shiftL` 2 -  exp x = let CR o = x / ln2-              l = o 0-              y = x - fromInteger l * ln2+  exp x = let o = shiftL (x *. recipBounded (shiftL ln2 1)) 1+              l = atPrecision o 0+              y = x - fromInteger l *. ln2           in if l == 0                then expBounded x                else expBounded y `shiftL` fromInteger l    -- | Range reduction on the principle that ln (a * b) = ln a + ln b-  log x = let CR o = x-              l = log2 (o 2) - 2+  log x = let l = log2 (atPrecision x 2) - 2           in if   -- x <= 0.75                 | l < 0  -> - log (recip x)                   -- 0.75 <= x <= 2@@ -189,20 +271,20 @@                 | l > 0  -> let a = x `shiftR` l                             in logBounded a + fromIntegral l *. ln2 -  sqrt (CR x) = crMemoize (\p -> let n = x (2 * p)-                                 in isqrt n)+  sqrt x = crMemoize (\p -> let n = atPrecision x (2 * p)+                            in isqrt n)    -- | This will diverge when the base is not positive   x ** y = exp (log x * y)    logBase x y = log y / log x -  sin x = cos (x - pi / 2)+  sin x = cos (x - piBy2) -  cos x = let CR o = x / piBy4-              s = o 1 /. 2-              octant = fromInteger $ s `mod` 8-              offset = x - (fromIntegral s * piBy4)+  cos x = let o = shiftL (x *. recipBounded pi) 2+              s = atPrecision o 1 /^ 1+              octant = fromInteger $ s .&. 7+              offset = x - (fromIntegral s *. piBy4)               fs = [          cosBounded                    , negate . sinBounded . subtract piBy4                    , negate . sinBounded@@ -215,46 +297,74 @@    tan x = sin x .* recip (cos x) -  asin x = 2 * atan (x .*. recipBounded (1 + sqrt (1 - x.*.x)))+  asin x = (atan (x .*. recipBounded (1 + sqrt (1 - squareBounded x)))) `shiftL` 1 -  acos x = pi/2 - asin x+  acos x = piBy2 - asin x    atan x = let -- q is 4 times x to within 1/4                q = x `atPrecision` 2            in if   -- x <= -1-                 | q <  -4 -> atanBounded (negate (recipBounded x)) - pi / 2+                 | q <  -4 -> atanBounded (negate (recipBounded x)) - piBy2                    -- -1.25 <= x <= -0.75-                 | q == -4 -> -pi / 4 - atanBounded ((x + 1) .*. recipBounded (x - 1))+                 | q == -4 -> -(piBy4 + atanBounded ((x + 1) .*. recipBounded (x - 1)))                    -- 0.75 <= x <= 1.25-                 | q ==  4 -> pi / 4 + atanBounded ((x - 1) .*. recipBounded (x + 1))+                 | q ==  4 -> piBy4 + atanBounded ((x - 1) .*. recipBounded (x + 1))                    -- x >= 1-                 | q >   4 -> pi / 2 - atanBounded (recipBounded x)+                 | q >   4 -> piBy2 - atanBounded (recipBounded x)                    -- -0.75 <= x <= 0.75                  | otherwise -> atanBounded x    -- TODO: benchmark replacing these with their series expansion   sinh x = let (expX, expNegX) = expPosNeg x-           in (expX - expNegX) / 2+           in (expX - expNegX) `shiftR` 1   cosh x = let (expX, expNegX) = expPosNeg x-           in (expX + expNegX) / 2-  tanh x = let e2x = exp (2 * x)+           in (expX + expNegX) `shiftR` 1+  tanh x = let e2x = exp (x `shiftL` 1)            in (e2x - 1) *. recipBounded (e2x + 1) -  asinh x = log (x + sqrt (x * x + 1))+  asinh x = log (x + sqrt (square x + 1))   acosh x = log (x + sqrt (x + 1) * sqrt (x - 1))-  atanh x = (log (1 + x) - log (1 - x)) / 2+  atanh x = (log (1 + x) - log (1 - x)) `shiftR` 1  -- | 'toRational' returns the CReal n evaluated at a precision of 2^-n instance KnownNat n => Real (CReal n) where   toRational x = let p = crealPrecision x-                 in x `atPrecision` p % 2^p+                 in x `atPrecision` p % bit p  instance KnownNat n => RealFrac (CReal n) where   properFraction x = let p = crealPrecision x-                         n = (x `atPrecision` p) `quot` 2^p-                         f =  x - fromInteger n+                         v = x `atPrecision` p+                         r = v .&. (bit p - 1)+                         c = unsafeShiftR (v - r) p+                         n = if c < 0 && r /= 0 then c + 1 else c+                         f = plusInteger x (negate n)                      in (fromInteger n, f) +  truncate x = let p = crealPrecision x+                   v = x `atPrecision` p+                   r = v .&. (bit p - 1)+                   c = unsafeShiftR (v - r) p+                   n = if c < 0 && r /= 0 then c + 1 else c+                   in fromInteger n++  round x = let p = crealPrecision x+                n = (x `atPrecision` p) /^ p+                in fromInteger n++  ceiling x = let p = crealPrecision x+                  v = x `atPrecision` p+                  r = v .&. (bit p - 1)+                  n = unsafeShiftR (v - r) p+                  in case r /= 0 of+                    True -> fromInteger $ n + 1+                    _    -> fromInteger n++  floor x = let p = crealPrecision x+                v = x `atPrecision` p+                r = v .&. (bit p - 1)+                n = unsafeShiftR (v - r) p+                in fromInteger 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.@@ -267,7 +377,9 @@   floatRange _ = error "Data.CReal.Internal floatRange"   decodeFloat x = let p = crealPrecision x                   in (x `atPrecision` p, -p)-  encodeFloat m n = fromRational (m % 2^(-n))+  encodeFloat m n = if n <= 0+    then fromRational (m % bit (negate n))+    else fromRational (unsafeShiftL m n :% 1)   exponent = error "Data.CReal.Internal exponent"   significand = error "Data.CReal.Internal significand"   scaleFloat = flip shiftL@@ -280,7 +392,7 @@     let y' = y `atPrecision` p         x' = x `atPrecision` p         θ = if | x' > 0            ->  atan (y/x)-               | x' == 0 && y' > 0 ->  pi/2+               | x' == 0 && y' > 0 ->  piBy2                | x' <  0 && y' > 0 ->  pi + atan (y/x)                | x' <= 0 && y' < 0 -> -atan2 (-y) x                | y' == 0 && x' < 0 ->  pi    -- must be after the previous test on zero y@@ -295,15 +407,17 @@ -- True instance KnownNat n => Eq (CReal n) where   -- TODO, should this try smaller values first?-  x == y = let p = crealPrecision x-           in (x - y) `atPrecision` p == 0+  CR mvx _ == CR mvy _ | mvx == mvy = True+  x == y = let p = crealPrecision x + 2+           in (atPrecision x p - atPrecision y p) /^ 2 == 0  -- | Like equality values of type @CReal p@ are compared at precision @p@. instance KnownNat n => Ord (CReal n) where-  compare x y = let p = crealPrecision x-                in compare ((x - y) `atPrecision` p) 0-  max (CR x) (CR y) = crMemoize (\p -> max (x p) (y p))-  min (CR x) (CR y) = crMemoize (\p -> min (x p) (y p))+  compare (CR mvx _) (CR mvy _) | mvx == mvy = EQ+  compare x y = let p = crealPrecision x + 2+                in compare ((atPrecision x p - atPrecision y p) /^ 2) 0+  max x y = crMemoize (\p -> max (atPrecision x p) (atPrecision y p))+  min x y = crMemoize (\p -> min (atPrecision x p) (atPrecision y p))  -- | The 'Random' instance for @\'CReal\' p@ will return random number with at -- least @p@ digits of precision, every digit after that is zero.@@ -328,8 +442,11 @@ --  piBy4 :: CReal n-piBy4 = 4 * atanBounded (recipBounded 5) - atanBounded (recipBounded 239) -- Machin Formula+piBy4 = (atanBounded (recipBounded 5) `shiftL` 2) - atanBounded (recipBounded 239) -- Machin Formula +piBy2 :: CReal n+piBy2 = piBy4 `shiftL` 1+ ln2 :: CReal n ln2 = logBounded 2 @@ -347,41 +464,64 @@ (*.) :: CReal n -> CReal n -> CReal n (*.) = flip mulBoundedL --- | Alias for @'mulBoundedL'@+-- | Alias for @'mulBounded'@ (.*.) :: 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 (CR x1) (CR x2) = crMemoize (\p -> let s1 = 4-                                                   s2 = log2 (abs (x2 0) + 2) + 3-                                                   n1 = x1 (p + s2)-                                                   n2 = x2 (p + s1)-                                               in (n1 * n2) /. 2^(p + s1 + s2))+mulBoundedL x1 x2 = let+  s1 = 4+  s2 = log2 (abs (atPrecision x2 0) + 2) + 3+  in crMemoize (\p -> let n1 = atPrecision x1 (p + s2)+                          n2 = atPrecision x2 (p + s1)+                      in (n1 * n2) /^ (p + s1 + 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 (CR x1) (CR x2) = crMemoize (\p -> let s1 = 4-                                                  s2 = 4-                                                  n1 = x1 (p + s2)-                                                  n2 = x2 (p + s1)-                                              in (n1 * n2) /. 2^(p + s1 + s2))+mulBounded x1 x2 = let+  s1 = 4+  s2 = 4+  in crMemoize (\p -> let n1 = atPrecision x1 (p + s2)+                          n2 = atPrecision x2 (p + s1)+                      in (n1 * n2) /^ (p + s1 + 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 (CR x) = crMemoize (\p -> let s = 2-                                           n = x (p + 2 * s + 2)-                                       in 2^(2 * p + 2 * s + 2) /. n)+recipBounded x = crMemoize (\p -> let s = 2+                                      n = atPrecision x (p + 2 * s + 2)+                                  in bit (2 * p + 2 * s + 2) /. n)  -- | Return the square of the input, more efficient than @('*')@+{-# INLINABLE square #-} square :: CReal n -> CReal n-square (CR x) = crMemoize (\p -> let s = log2 (abs (x 0) + 2) + 3-                                     n = x (p + s)-                                 in (n * n) /. 2^(p + 2 * s))+square x = let+  s = log2 (abs (atPrecision x 0) + 2) + 3+  in crMemoize (\p -> let n = atPrecision x (p + s)+                      in (n * n) /^ (p + 2 * 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 x@(CR mvc _) = unsafePerformIO $ do+  vcc <- tryReadMVar mvc+  let+    !s = 4+    !vcn = case vcc of+      Nothing -> Never+      Just Never -> Never+      Just (Current j n) -> case j - s of+        p | p < 0 -> Never+        p -> Current p ((n * n) /^ (p + 2 * s))+    fn' !p = let n = atPrecision x (p + s)+             in (n * n) /^ (p + 2 * s)+  mvn <- newMVar vcn+  return $ CR mvn fn'+ -- -- Bounded exponential functions and expPosNeg --@@ -401,8 +541,8 @@  -- | @expPosNeg x@ returns @(exp x, exp (-x))# expPosNeg :: CReal n -> (CReal n, CReal n)-expPosNeg x = let CR o = x / ln2-                  l = o 0+expPosNeg x = let o = x / ln2+                  l = atPrecision o 0                   y = x - fromInteger l * ln2               in if l == 0                    then (expBounded x, expBounded (-x))@@ -417,23 +557,44 @@ -- closed range [-1..1] sinBounded :: CReal n -> CReal n sinBounded x = let q = alternateSign (scanl' (*) 1 [ 1 % (n*(n+1)) | n <- [2,4..]])-               in x * powerSeries q (max 1) (x .*. x)+               in x .* powerSeries q (max 1) (squareBounded 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 x = let q = alternateSign (scanl' (*) 1 [1 % (n*(n+1)) | n <- [1,3..]])-               in powerSeries q (max 1) (x .*. x)+               in powerSeries q (max 1) (squareBounded 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 x = let q = scanl' (*) 1 [n % (n + 1) | n <- [2,4..]]-                    d = 1 + x .*. x-                    rd = recipBounded d-                in (x .*. rd) .* powerSeries q (+1) (x .*. x .*. rd)+                    s = squareBounded x+                    rd = recipBounded (plusInteger s 1)+                in (x .*. rd) .* powerSeries q (+1) (s .*. 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 x 0 = x+plusInteger (CR mvc fn) n = unsafePerformIO $ do+  vcc <- tryReadMVar mvc+  let+    !vcn = case vcc of+      Nothing -> Never+      Just Never -> Never+      Just (Current j v) -> Current j (v + unsafeShiftL n j)+    fn' !p = fn p + B.shiftL n p+  mvc' <- newMVar vcn+  return $ CR mvc' fn'++-- -- Multiplication with powers of two -- @@ -445,10 +606,10 @@ -- -- This can be faster than doing the division shiftR :: CReal n -> Int -> CReal n-shiftR (CR x) n = crMemoize (\p -> let p' = p - n-                                   in if p' >= 0-                                        then x p'-                                        else x 0 /. 2^(-p'))+shiftR x n = crMemoize (\p -> let p' = p - n+                              in if p' >= 0+                                 then atPrecision x p'+                                 else atPrecision x 0 /^ (negate p'))  -- | @x \`shiftL\` n@ is equal to @x@ multiplied by 2^@n@ --@@ -466,77 +627,110 @@ -- | 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 p (CR x) = let places = decimalDigitsAtPrecision p-                               r = x p % 2^p-                           in rationalToDecimal places r+showAtPrecision p x = let places = decimalDigitsAtPrecision p+                          r      = atPrecision x p % bit p+                      in rationalToDecimal places r  -- | How many decimal digits are required to represent a number to within 2^-p decimalDigitsAtPrecision :: Int -> Int decimalDigitsAtPrecision 0 = 0-decimalDigitsAtPrecision p = log10 (2^p) + 1+decimalDigitsAtPrecision p = log10 (bit p) + 1  -- | @rationalToDecimal p x@ returns a string representing @x@ at @p@ decimal -- places. rationalToDecimal :: Int -> Rational -> String-rationalToDecimal places r = p ++ is ++ if places > 0 then "." ++ fs else ""-  where r' = abs r-        p = case signum r of+rationalToDecimal places (n :% d) = p ++ is ++ if places > 0 then "." ++ fs else ""+  where p = case signum n of               -1 -> "-"               _  -> ""-        ds = show ((numerator r' * 10^places) /. denominator r')+        ds = show (roundD (abs n * 10^places) d)         l = length ds         (is, fs) = if | l <= places -> ("0", replicate (places - l) '0' ++ ds)-                      | otherwise -> splitAt (length ds - places) ds+                      | otherwise -> splitAt (l - places) ds   -- -- Integer operations -- +divZeroErr :: a+divZeroErr = error "Division by zero"+{-# NOINLINE divZeroErr #-}++roundD :: Integer -> Integer -> Integer+roundD n d = case divMod n d of+  (q, r) -> case compare (unsafeShiftL r 1) d of+    LT -> q+    EQ -> if testBit q 0 then q + 1 else q+    GT -> q + 1+{-# INLINE roundD #-}+ infixl 7 /. -- | Division rounding to the nearest integer and rounding half integers to the -- nearest even integer. (/.) :: Integer -> Integer -> Integer-n /. d = round (n % d)+(!n) /. (!d) = case compare d 0 of+  LT -> roundD (negate n) (negate d)+  EQ -> divZeroErr+  GT -> roundD n 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+(!n) /^ (!p) = case compare p 0 of+  LT -> unsafeShiftL n (negate p)+  EQ -> n+  GT -> let+    !bp = bit p+    !r = n .&. (bp - 1)+    !q = unsafeShiftR (n - r) p+    in case compare (unsafeShiftL r 1) bp of+      LT -> q+      EQ -> if testBit q 0 then q + 1 else q+      GT -> q + 1+ -- | @log2 x@ returns the base 2 logarithm of @x@ rounded towards zero. -- -- The input must be positive+{-# INLINE log2 #-} log2 :: Integer -> Int log2 x = I# (integerLog2# x)  -- | @log10 x@ returns the base 10 logarithm of @x@ rounded towards zero. -- -- The input must be positive+{-# INLINE log10 #-} log10 :: Integer -> Int log10 x = I# (integerLogBase# 10 x)  -- | @isqrt x@ returns the square root of @x@ rounded towards zero. -- -- The input must not be negative+{-# INLINABLE isqrt #-} isqrt :: Integer -> Integer isqrt x | x < 0     = error "Sqrt applied to negative Integer"         | x == 0    = 0         | otherwise = until satisfied improve initialGuess-  where improve r    = (r + (x `div` r)) `div` 2-        satisfied r  = sq r <= x && sq (r + 1) > x-        initialGuess = 2 ^ (log2 x `div` 2)-        sq r         = r * r+  where improve r    = unsafeShiftR (r + (x `div` r)) 1+        satisfied r  = let r2 = r * r in r2 <= x && r2 + unsafeShiftL r 1 >= x+        initialGuess = bit (unsafeShiftR (log2 x) 1)  -- -- Searching --  -- | Given a monotonic function+{-# INLINABLE findFirstMonotonic #-} findFirstMonotonic :: (Int -> Bool) -> Int-findFirstMonotonic p = binarySearch l' u'-  where (l', u') = findBounds 0 1-        findBounds l u = if p u then (l, u)-                                else findBounds u (u*2)-        binarySearch l u = let m = l + ((u - l) `div` 2)-                           in if | l+1 == u  -> l-                                 | p m       -> binarySearch l m-                                 | otherwise -> binarySearch m u+findFirstMonotonic p = findBounds 0 1+  where findBounds !l !u = if p u then binarySearch l u+                                  else findBounds u (u * 2)+        binarySearch !l !u = let !m = l + ((u - l) `div` 2)+                             in if | l+1 == u  -> l+                                   | p m       -> binarySearch l m+                                   | otherwise -> binarySearch m u   --@@ -547,8 +741,9 @@ -- -- >>> alternateSign [1..5] -- [1,-2,3,-4,5]+{-# INLINABLE alternateSign #-} alternateSign :: Num a => [a] -> [a]-alternateSign = zipWith ($) (cycle [id, negate])+alternateSign = \ls -> foldr (\a r b -> if b then (negate a):r False else a:r True) (const []) ls 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@@@ -564,13 +759,13 @@ -- >>> powerSeries [1 % (n!) | n <- [0..]] (max 5) 1 :: CReal 218 -- 2.718281828459045235360287471352662497757247093699959574966967627724 powerSeries :: [Rational] -> (Int -> Int) -> CReal n -> CReal n-powerSeries q termsAtPrecision (CR x) = crMemoize+powerSeries q termsAtPrecision x = crMemoize      (\p -> let t = termsAtPrecision p                 d = log2 (toInteger t) + 2                 p' = p + d                 p'' = p' + d-                m = x p''-                xs = (%1) <$> iterate (\e -> m * e /. 2^p'') (2^p')-                r = sum . take (t + 1) . fmap (round . (* (2^d))) $ zipWith (*) q xs-            in r /. 4^d)+                m = atPrecision x p''+                xs = (%1) <$> iterate (\e -> m * e /^ p'') (bit p')+                r = sum . take (t + 1) . fmap (round . (* fromInteger (bit d))) $ zipWith (*) q xs+            in r /^ (2 * d))