diff --git a/.gitignore b/.gitignore --- a/.gitignore +++ b/.gitignore @@ -14,3 +14,5 @@ *.aux *.hp .stack-work/ +result +.hdevtools diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -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 diff --git a/exact-real.cabal b/exact-real.cabal --- a/exact-real.cabal +++ b/exact-real.cabal @@ -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 diff --git a/src/Data/CReal/Internal.hs b/src/Data/CReal/Internal.hs --- a/src/Data/CReal/Internal.hs +++ b/src/Data/CReal/Internal.hs @@ -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))