arithmoi 0.13.1.0 → 0.13.2.0
raw patch · 15 files changed
+210/−289 lines, 15 filesdep +bitvecdep −arraydep ~deepseqdep ~moddep ~randomPVP ok
version bump matches the API change (PVP)
Dependencies added: bitvec
Dependencies removed: array
Dependency ranges changed: deepseq, mod, random, tasty-bench, transformers
API changes (from Hackage documentation)
Files
- Math/NumberTheory/DirichletCharacters.hs +2/−0
- Math/NumberTheory/Moduli/JacobiSymbol.hs +10/−10
- Math/NumberTheory/Moduli/Singleton.hs +1/−0
- Math/NumberTheory/Moduli/SomeMod.hs +5/−6
- Math/NumberTheory/Primes/Counting/Impl.hs +81/−159
- Math/NumberTheory/Primes/Factorisation/Montgomery.hs +4/−3
- Math/NumberTheory/Primes/Sieve/Eratosthenes.hs +72/−84
- Math/NumberTheory/Primes/Sieve/Indexing.hs +4/−4
- Math/NumberTheory/Recurrences/Bilinear.hs +1/−2
- Math/NumberTheory/Utils/DirichletSeries.hs +13/−5
- arithmoi.cabal +5/−9
- benchmark/Math/NumberTheory/SequenceBench.hs +4/−4
- changelog.md +6/−0
- test-suite/Math/NumberTheory/DirichletCharactersTests.hs +1/−2
- test-suite/Math/NumberTheory/Moduli/ClassTests.hs +1/−1
Math/NumberTheory/DirichletCharacters.hs view
@@ -15,6 +15,8 @@ {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE ViewPatterns #-} +{-# OPTIONS_GHC -Wno-unrecognised-warning-flags -Wno-pattern-namespace-specifier #-}+ module Math.NumberTheory.DirichletCharacters ( -- * An absorbing semigroup
Math/NumberTheory/Moduli/JacobiSymbol.hs view
@@ -41,11 +41,6 @@ Zero -> Zero One -> MinusOne -{-# SPECIALISE symbolToNum :: JacobiSymbol -> Integer,- JacobiSymbol -> Int,- JacobiSymbol -> Word,- JacobiSymbol -> Natural- #-} -- | Convenience function to convert out of a Jacobi symbol symbolToNum :: Num a => JacobiSymbol -> a symbolToNum = \case@@ -53,6 +48,11 @@ One -> 1 MinusOne -> -1 +{-# SPECIALISE symbolToNum :: JacobiSymbol -> Integer #-}+{-# SPECIALISE symbolToNum :: JacobiSymbol -> Int #-}+{-# SPECIALISE symbolToNum :: JacobiSymbol -> Word #-}+{-# SPECIALISE symbolToNum :: JacobiSymbol -> Natural #-}+ -- | <https://en.wikipedia.org/wiki/Jacobi_symbol Jacobi symbol> of two arguments. -- The lower argument (\"denominator\") must be odd and positive, -- this condition is checked.@@ -64,17 +64,17 @@ -- Zero -- >>> jacobi 1001 9907 -- MinusOne-{-# SPECIALISE jacobi :: Integer -> Integer -> JacobiSymbol,- Natural -> Natural -> JacobiSymbol,- Int -> Int -> JacobiSymbol,- Word -> Word -> JacobiSymbol- #-} jacobi :: (Integral a, Bits a) => a -> a -> JacobiSymbol jacobi _ 1 = One jacobi a b | b < 0 = error "Math.NumberTheory.Moduli.jacobi: negative denominator" | evenI b = error "Math.NumberTheory.Moduli.jacobi: even denominator" | otherwise = jacobi' a b -- b odd, > 1++{-# SPECIALISE jacobi :: Integer -> Integer -> JacobiSymbol #-}+{-# SPECIALISE jacobi :: Natural -> Natural -> JacobiSymbol #-}+{-# SPECIALISE jacobi :: Int -> Int -> JacobiSymbol #-}+{-# SPECIALISE jacobi :: Word -> Word -> JacobiSymbol #-} jacobi' :: (Integral a, Bits a) => a -> a -> JacobiSymbol jacobi' 0 _ = Zero
Math/NumberTheory/Moduli/Singleton.hs view
@@ -19,6 +19,7 @@ {-# LANGUAGE TupleSections #-} {-# LANGUAGE TypeOperators #-} {-# LANGUAGE ViewPatterns #-}+{-# OPTIONS_GHC -Wno-unrecognised-warning-flags -Wno-pattern-namespace-specifier #-} module Math.NumberTheory.Moduli.Singleton ( -- * SFactors singleton
Math/NumberTheory/Moduli/SomeMod.hs view
@@ -20,7 +20,7 @@ ) where import Data.Euclidean (GcdDomain(..), Euclidean(..), Field)-import Data.Mod+import Data.Mod (Mod, unMod, invertMod, (^%)) import Data.Proxy import Data.Semiring (Semiring(..), Ring(..)) import Data.Type.Equality@@ -180,11 +180,10 @@ InfMod r -> Just (InfMod (recip r)) {-# INLINABLE [1] invertSomeMod #-} -{-# SPECIALISE [1] powSomeMod ::- SomeMod -> Integer -> SomeMod,- SomeMod -> Natural -> SomeMod,- SomeMod -> Int -> SomeMod,- SomeMod -> Word -> SomeMod #-}+{-# SPECIALISE [1] powSomeMod :: SomeMod -> Integer -> SomeMod #-}+{-# SPECIALISE [1] powSomeMod :: SomeMod -> Natural -> SomeMod #-}+{-# SPECIALISE [1] powSomeMod :: SomeMod -> Int -> SomeMod #-}+{-# SPECIALISE [1] powSomeMod :: SomeMod -> Word -> SomeMod #-} -- | Drop-in replacement for 'Prelude.^', with much better performance. -- When -O is enabled, there is a rewrite rule, which specialises 'Prelude.^' to 'powSomeMod'.
Math/NumberTheory/Primes/Counting/Impl.hs view
@@ -28,11 +28,11 @@ import Math.NumberTheory.Utils.FromIntegral import Control.Monad.ST-import Data.Array.Base-import Data.Array.ST+import Data.Bit import Data.Bits import Data.Int-import Unsafe.Coerce+import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Unboxed.Mutable as MU -- | Maximal allowed argument of 'primeCount'. Currently 8e18. primeCountMaxArg :: Integer@@ -55,7 +55,7 @@ | n < 1000 = intToInteger . length . takeWhile (<= n) . map unPrime . primeList . primeSieve $ max 242 n | n < 30000 = runST $ do ba <- sieveTo n- (s,e) <- getBounds ba+ let (s, e) = (0, MU.length ba - 1) ct <- countFromTo s e ba return (intToInteger $ ct+3) | otherwise =@@ -137,7 +137,7 @@ strt = a + 1 + (a .&. 1) psieves@(PS vO ba:_) = psieveFrom (toInteger strt) rep | o0 < 0 = 0- | otherwise = sum [1 | i <- [0 .. r2], ba `unsafeAt` i]+ | otherwise = sum [1 | i <- [0 .. r2], unBit (ba `U.unsafeIndex` i)] where o0 = toInteger strt - vO - 9 -- (strt - 2) - v0 - 7 r0 = fromInteger o0 `rem` 30@@ -157,32 +157,32 @@ picr = approxPrimeCount cr diff = pisr - picr size = int64ToInt (diff + diff `quot` 50) + 30- store <- unsafeNewArray_ (0,size-1) :: ST s (STUArray s Int Int64)- let feed :: Int64 -> Int -> Int -> UArray Int Bool -> [PrimeSieve] -> ST s Integer+ store <- MU.unsafeNew size :: ST s (MU.MVector s Int64)+ let feed :: Int64 -> Int -> Int -> U.Vector Bit -> [PrimeSieve] -> ST s Integer feed voff !wi !ri uar sves | ri == sieveBits = case sves of (PS vO ba : more) -> feed (fromInteger vO) wi 0 ba more _ -> error "prime stream ended prematurely" | pval > sr = do- stu <- unsafeThaw uar+ stu <- U.unsafeThaw uar eat 0 0 voff (wi-1) ri stu sves- | uar `unsafeAt` ri = do- unsafeWrite store wi (ub `quot` pval)+ | unBit (uar `U.unsafeIndex` ri) = do+ MU.unsafeWrite store wi (ub `quot` pval) feed voff (wi+1) (ri+1) uar sves | otherwise = feed voff wi (ri+1) uar sves where pval = voff + toPrim ri- eat :: Integer -> Integer -> Int64 -> Int -> Int -> STUArray s Int Bool -> [PrimeSieve] -> ST s Integer+ eat :: Integer -> Integer -> Int64 -> Int -> Int -> MU.MVector s Bit -> [PrimeSieve] -> ST s Integer eat !acc !btw voff !wi !si stu sves | si == sieveBits = case sves of [] -> error "Premature end of prime stream" (PS vO ba : more) -> do- nstu <- unsafeThaw ba+ nstu <- U.unsafeThaw ba eat acc btw (fromInteger vO) wi 0 nstu more | wi < 0 = return acc | otherwise = do- qb <- unsafeRead store wi+ qb <- MU.unsafeRead store wi let dist = qb - voff - 7 if dist < intToInt64 sieveRange then do@@ -196,7 +196,7 @@ (b,j) = idxPr (fds+7) !li = (b `shiftL` 3) .|. j ctLoop !lac 0 (PS vO ba : more) = do- nstu <- unsafeThaw ba+ nstu <- U.unsafeThaw ba new <- countFromTo 0 li nstu let nbtw = btw + lac + 1 + intToInteger new eat (acc+nbtw) nbtw (integerToInt64 vO) (wi-1) (li+1) nstu more@@ -216,29 +216,29 @@ calcST :: forall s. Int64 -> Int64 -> ST s Integer calcST lim plim = do !parr <- sieveTo (int64ToInteger plim)- (plo,phi) <- getBounds parr+ let (plo, phi) = (0, MU.length parr - 1) !pct <- countFromTo plo phi parr- !ar1 <- unsafeNewArray_ (0,end-1)- unsafeWrite ar1 0 lim- unsafeWrite ar1 1 1- !ar2 <- unsafeNewArray_ (0,end-1)- let go :: Int -> Int -> STUArray s Int Int64 -> STUArray s Int Int64 -> ST s Integer+ !ar1 <- MU.unsafeNew end+ MU.unsafeWrite ar1 0 lim+ MU.unsafeWrite ar1 1 1+ !ar2 <- MU.unsafeNew end+ let go :: Int -> Int -> MU.MVector s Int64 -> MU.MVector s Int64 -> ST s Integer go cap pix old new | pix == 2 = coll cap old | otherwise = do- isp <- unsafeRead parr pix+ Bit isp <- MU.unsafeRead parr pix if isp then do let !n = fromInteger (toPrim pix) !ncap <- treat cap n old new go ncap (pix-1) new old else go cap (pix-1) old new- coll :: Int -> STUArray s Int Int64 -> ST s Integer+ coll :: Int -> MU.MVector s Int64 -> ST s Integer coll stop ar = let cgo !acc i | i < stop = do- !k <- unsafeRead ar i- !v <- unsafeRead ar (i+1)+ !k <- MU.unsafeRead ar i+ !v <- MU.unsafeRead ar (i+1) cgo (acc + int64ToInteger v*cp6 k) (i+2) | otherwise = return (acc+intToInteger pct+2) in cgo 0 0@@ -249,37 +249,37 @@ !size = int64ToInt $ integerSquareRoot lim `quot` 4 !end = 2*size -treat :: Int -> Int64 -> STUArray s Int Int64 -> STUArray s Int Int64 -> ST s Int+treat :: Int -> Int64 -> MU.MVector s Int64 -> MU.MVector s Int64 -> ST s Int treat end n old new = do qi0 <- locate n 0 (end `quot` 2 - 1) old let collect stop !acc ix | ix < end = do- !k <- unsafeRead old ix+ !k <- MU.unsafeRead old ix if k < stop then do- v <- unsafeRead old (ix+1)+ v <- MU.unsafeRead old (ix+1) collect stop (acc-v) (ix+2) else return (acc,ix) | otherwise = return (acc,ix) goTreat !wi !ci qi | qi < end = do- !key <- unsafeRead old qi- !val <- unsafeRead old (qi+1)+ !key <- MU.unsafeRead old qi+ !val <- MU.unsafeRead old (qi+1) let !q0 = key `quot` n !r0 = int64ToInt (q0 `rem` 30030)- !nkey = q0 - int8ToInt64 (cpDfAr `unsafeAt` r0)- nk0 = q0 + int8ToInt64 (cpGpAr `unsafeAt` (r0+1) + 1)+ !nkey = q0 - int8ToInt64 (cpDfAr `U.unsafeIndex` r0)+ nk0 = q0 + int8ToInt64 (cpGpAr `U.unsafeIndex` (r0+1) + 1) !nlim = n*nk0 (wi1,ci1) <- copyTo end nkey old ci new wi- ckey <- unsafeRead old ci1+ ckey <- MU.unsafeRead old ci1 (!acc, !ci2) <- if ckey == nkey then do- !ov <- unsafeRead old (ci1+1)+ !ov <- MU.unsafeRead old (ci1+1) return (ov-val,ci1+2) else return (-val,ci1) (!tot, !nqi) <- collect nlim acc (qi+2)- unsafeWrite new wi1 nkey- unsafeWrite new (wi1+1) tot+ MU.unsafeWrite new wi1 nkey+ MU.unsafeWrite new (wi1+1) tot goTreat (wi1+2) ci2 nqi | otherwise = copyRem end old ci new wi goTreat 0 0 qi0@@ -288,12 +288,12 @@ -- Auxiliaries -- -------------------------------------------------------------------------------- -locate :: Int64 -> Int -> Int -> STUArray s Int Int64 -> ST s Int+locate :: Int64 -> Int -> Int -> MU.MVector s Int64 -> ST s Int locate p low high arr = do let go lo hi | lo < hi = do let !md = (lo+hi) `quot` 2- v <- unsafeRead arr (2*md)+ v <- MU.unsafeRead arr (2*md) case compare p v of LT -> go lo md EQ -> return (2*md)@@ -302,58 +302,55 @@ go low high {-# INLINE copyTo #-}-copyTo :: Int -> Int64 -> STUArray s Int Int64 -> Int- -> STUArray s Int Int64 -> Int -> ST s (Int,Int)+copyTo :: Int -> Int64 -> MU.MVector s Int64 -> Int+ -> MU.MVector s Int64 -> Int -> ST s (Int,Int) copyTo end lim old oi new ni = do let go ri wi | ri < end = do- ok <- unsafeRead old ri+ ok <- MU.unsafeRead old ri if ok < lim then do- !ov <- unsafeRead old (ri+1)- unsafeWrite new wi ok- unsafeWrite new (wi+1) ov+ !ov <- MU.unsafeRead old (ri+1)+ MU.unsafeWrite new wi ok+ MU.unsafeWrite new (wi+1) ov go (ri+2) (wi+2) else return (wi,ri) | otherwise = return (wi,ri) go oi ni {-# INLINE copyRem #-}-copyRem :: Int -> STUArray s Int Int64 -> Int -> STUArray s Int Int64 -> Int -> ST s Int+copyRem :: Int -> MU.MVector s Int64 -> Int -> MU.MVector s Int64 -> Int -> ST s Int copyRem end old oi new ni = do- let go ri wi- | ri < end = do- unsafeRead old ri >>= unsafeWrite new wi- go (ri+1) (wi+1)- | otherwise = return wi- go oi ni+ let len = end - oi+ MU.copy (MU.slice ni len new) (MU.slice oi len old)+ pure $ ni + len {-# INLINE cp6 #-} cp6 :: Int64 -> Integer cp6 k = case k `quotRem` 30030 of (q,r) -> 5760*int64ToInteger q +- int16ToInteger (cpCtAr `unsafeAt` int64ToInt r)+ int16ToInteger (cpCtAr `U.unsafeIndex` int64ToInt r) cop :: Int64 -> Int64-cop m = m - int8ToInt64 (cpDfAr `unsafeAt` int64ToInt (m `rem` 30030))+cop m = m - int8ToInt64 (cpDfAr `U.unsafeIndex` int64ToInt (m `rem` 30030)) -------------------------------------------------------------------------------- -- Ugly helper arrays -- -------------------------------------------------------------------------------- -cpCtAr :: UArray Int Int16-cpCtAr = runSTUArray $ do- ar <- newArray (0,30029) 1+cpCtAr :: U.Vector Int16+cpCtAr = runST $ do+ ar <- MU.replicate 30030 1 let zilch s i- | i < 30030 = unsafeWrite ar i 0 >> zilch s (i+s)+ | i < 30030 = MU.unsafeWrite ar i 0 >> zilch s (i+s) | otherwise = return () accumulate ct i | i < 30030 = do- v <- unsafeRead ar i+ v <- MU.unsafeRead ar i let !ct' = ct+v- unsafeWrite ar i ct'+ MU.unsafeWrite ar i ct' accumulate ct' (i+1) | otherwise = return ar zilch 2 0@@ -362,20 +359,21 @@ zilch 14 7 zilch 22 11 zilch 26 13- accumulate 1 2+ vec <- accumulate 1 2+ U.unsafeFreeze vec -cpDfAr :: UArray Int Int8-cpDfAr = runSTUArray $ do- ar <- newArray (0,30029) 0+cpDfAr :: U.Vector Int8+cpDfAr = runST $ do+ ar <- MU.replicate 30030 0 let note s i- | i < 30029 = unsafeWrite ar i 1 >> note s (i+s)+ | i < 30029 = MU.unsafeWrite ar i 1 >> note s (i+s) | otherwise = return () accumulate d i | i < 30029 = do- v <- unsafeRead ar i+ v <- MU.unsafeRead ar i if v == 0 then accumulate 2 (i+2)- else do unsafeWrite ar i d+ else do MU.unsafeWrite ar i d accumulate (d+1) (i+1) | otherwise = return ar note 2 0@@ -384,22 +382,23 @@ note 14 7 note 22 11 note 26 13- accumulate 2 3+ vec <- accumulate 2 3+ U.unsafeFreeze vec -cpGpAr :: UArray Int Int8-cpGpAr = runSTUArray $ do- ar <- newArray (0,30030) 0- unsafeWrite ar 30030 1+cpGpAr :: U.Vector Int8+cpGpAr = runST $ do+ ar <- MU.replicate 30031 0+ MU.unsafeWrite ar 30030 1 let note s i- | i < 30029 = unsafeWrite ar i 1 >> note s (i+s)+ | i < 30029 = MU.unsafeWrite ar i 1 >> note s (i+s) | otherwise = return () accumulate d i | i < 1 = return ar | otherwise = do- v <- unsafeRead ar i+ v <- MU.unsafeRead ar i if v == 0 then accumulate 2 (i-2)- else do unsafeWrite ar i d+ else do MU.unsafeWrite ar i d accumulate (d+1) (i-1) note 2 0 note 6 3@@ -407,103 +406,26 @@ note 14 7 note 22 11 note 26 13- accumulate 2 30027+ vec <- accumulate 2 30027+ U.unsafeFreeze vec ------------------------------------------------------------------------------- -- Prime counting -rMASK :: Int-rMASK = finiteBitSize (0 :: Word) - 1--wSHFT :: (Bits a, Num a) => a-wSHFT = if finiteBitSize (0 :: Word) == 64 then 6 else 5--tOPB :: Int-tOPB = finiteBitSize (0 :: Word) `shiftR` 1--tOPM :: (Bits a, Num a) => a-tOPM = (1 `shiftL` tOPB) - 1- -- find the n-th set bit in a list of PrimeSieves, -- aka find the (n+3)-rd prime countToNth :: Int -> [PrimeSieve] -> Integer countToNth !_ [] = error "countToNth: Prime stream ended prematurely"-countToNth !n (PS v0 bs : more) = go n 0- where- wa :: UArray Int Word- wa = coerceArray bs-- go !k i- | i > snd (bounds wa)- = countToNth k more- | otherwise- = let w = unsafeAt wa i- bc = popCount w- in if bc < k- then go (k-bc) (i+1)- else let j = bc - k- px = top w j bc- in v0 + toPrim (px + (i `shiftL` wSHFT))+countToNth !n (PS v0 bs : more) = case nthBitIndex (Bit True) n bs of+ Just i -> v0 + toPrim i+ Nothing -> countToNth (n - countBits bs) more -- count all set bits in a chunk, do it wordwise for speed. countAll :: PrimeSieve -> Int-countAll (PS _ bs) = go 0 0- where- wa :: UArray Int Word- wa = coerceArray bs-- go !ct i- | i > snd (bounds wa)- = ct- | otherwise- = go (ct + popCount (unsafeAt wa i)) (i+1)---- This is dangerous: we rely on the fact that reading a whole word--- at the end of ByteArray is permissible, even if its length--- is not aligned with word boundaries.-coerceArray :: UArray Int Bool -> UArray Int Word-coerceArray (UArray 0 u n ba)- | u + 1 == n, n' > 0 = UArray 0 (n' - 1) n' ba- where- n' = (n + rMASK) `shiftR` wSHFT-coerceArray (UArray l u n _) =- error $ "Cannot coerce array of Bool to array of Word: " ++ show (l, u, n)---- Find the j-th highest of bc set bits in the Word w.-top :: Word -> Int -> Int -> Int-top w j bc = go 0 tOPB tOPM bn w- where- !bn = bc-j- go !_ _ !_ !_ 0 = error "Too few bits set"- go bs 0 _ _ wd = if wd .&. 1 == 0 then error "Too few bits, shift 0" else bs- go bs a msk ix wd =- case popCount (wd .&. msk) of- lc | lc < ix -> go (bs+a) a msk (ix-lc) (wd `unsafeShiftR` a)- | otherwise ->- let !na = a `shiftR` 1- in go bs na (msk `unsafeShiftR` na) ix wd+countAll (PS _ bs) = countBits bs -- count set bits between two indices (inclusive) -- start and end must both be valid indices and start <= end-countFromTo :: Int -> Int -> STUArray s Int Bool -> ST s Int-countFromTo start end ba = do- wa <- (castSTUArray :: STUArray s Int Bool -> ST s (STUArray s Int Word)) ba- let !sb = start `shiftR` wSHFT- !si = start .&. rMASK- !eb = end `shiftR` wSHFT- !ei = end .&. rMASK- count !acc i- | i == eb = do- w <- unsafeRead wa i- return (acc + popCount (w `shiftL` (rMASK - ei)))- | otherwise = do- w <- unsafeRead wa i- count (acc + popCount w) (i+1)- if sb < eb- then do- w <- unsafeRead wa sb- count (popCount (w `shiftR` si)) (sb+1)- else do- w <- unsafeRead wa sb- let !w1 = w `shiftR` si- return (popCount (w1 `shiftL` (rMASK - ei + si)))+countFromTo :: Int -> Int -> MU.MVector s Bit -> ST s Int+countFromTo start end =+ fmap countBits . U.unsafeFreeze . MU.slice start (end - start + 1)
Math/NumberTheory/Primes/Factorisation/Montgomery.hs view
@@ -36,7 +36,7 @@ import Prelude hiding (Foldable(..)) import Control.Arrow import Control.Monad.Trans.State.Lazy-import Data.Array.Base (bounds, unsafeAt)+import Data.Bit import Data.Bits import Data.Foldable import Data.IntMap (IntMap)@@ -45,6 +45,7 @@ import Data.Mod import Data.Proxy import Data.Traversable+import qualified Data.Vector.Unboxed as U import GHC.Exts import GHC.Num.BigNat import GHC.Natural@@ -339,8 +340,8 @@ -- generate list of primes from arrays list :: [PrimeSieve] -> [Word]-list sieves = concat [[off + toPrim i | i <- [0 .. li], unsafeAt bs i]- | PS vO bs <- sieves, let { (_,li) = bounds bs; off = fromInteger vO; }]+list sieves = concat [[off + toPrim i | i <- [0 .. li], unBit (U.unsafeIndex bs i)]+ | PS vO bs <- sieves, let { li = U.length bs - 1; off = fromInteger vO; }] -- | @'smallFactors' n@ finds all prime divisors of @n > 1@ up to 2^16 by trial division and returns the -- list of these together with their multiplicities, and a possible remaining factor which may be composite.
Math/NumberTheory/Primes/Sieve/Eratosthenes.hs view
@@ -26,11 +26,12 @@ import Control.Monad (when) import Control.Monad.ST-import Data.Array.Base-import Data.Array.ST+import Data.Bit import Data.Bits import Data.Coerce import Data.Proxy+import qualified Data.Vector.Unboxed as U+import qualified Data.Vector.Unboxed.Mutable as MU import Data.Word import Math.NumberTheory.Primes.Sieve.Indexing@@ -74,11 +75,8 @@ sieveRange :: Int sieveRange = 30 * sieveBytes -wSHFT :: (Bits a, Num a) => a-wSHFT = if finiteBitSize (0 :: Word) == 64 then 6 else 5- -- | Compact store of primality flags.-data PrimeSieve = PS !Integer {-# UNPACK #-} !(UArray Int Bool)+data PrimeSieve = PS !Integer {-# UNPACK #-} !(U.Vector Bit) -- | Sieve primes up to (and including) a bound (or 7, if bound is smaller). -- For small enough bounds, this is more efficient than@@ -90,7 +88,7 @@ -- is often within memory limits, so don't give bounds larger than -- @8*10^9@ there. primeSieve :: Integer -> PrimeSieve-primeSieve bound = PS 0 (runSTUArray $ sieveTo bound)+primeSieve bound = PS 0 (runST $ sieveTo bound >>= U.unsafeFreeze) -- | Generate a list of primes for consumption from a -- 'PrimeSieve'.@@ -106,9 +104,9 @@ primeListInternal :: Num a => PrimeSieve -> [a] primeListInternal (PS v0 bs) = map ((+ fromInteger v0) . toPrim)- $ filter (unsafeAt bs) [lo..hi]+ $ filter (unBit . U.unsafeIndex bs) [lo..hi] where- (lo, hi) = bounds bs+ (lo, hi) = (0, U.length bs - 1) -- | Returns true if integer is beyond representation range of type a. doesNotFit :: forall a. Integral a => Proxy a -> Integer -> Bool@@ -157,13 +155,13 @@ where plim = 4801 -- prime #647, 644 of them to use sqlim = plim*plim- cache = runSTUArray $ do+ cache = runST $ do sieve <- sieveTo (4801 :: Integer)- new <- unsafeNewArray_ (0,1287) :: ST s (STUArray s Int Word64)+ new <- MU.unsafeNew 1288 :: ST s (MU.MVector s Word64) let fill j indx | 1279 < indx = return new -- index of 4801 = 159*30 + 31 ~> 159*8+7 | otherwise = do- p <- unsafeRead sieve indx+ Bit p <- MU.unsafeRead sieve indx if p then do let !i = indx .&. jMASK@@ -172,20 +170,21 @@ !strt = intToWord64 (strt1 .&. iXMASK) !skip = intToWord64 (strt1 `shiftR` iXBITS) !ixes = intToWord64 indx `shiftL` iXJBITS + strt `shiftL` jBITS + intToWord64 i- unsafeWrite new j skip- unsafeWrite new (j+1) ixes+ MU.unsafeWrite new j skip+ MU.unsafeWrite new (j+1) ixes fill (j+2) (indx+1) else fill j (indx+1)- fill 0 0+ vec <- fill 0 0+ U.unsafeFreeze vec -makeSieves :: Integer -> Integer -> Integer -> Integer -> UArray Int Word64 -> [PrimeSieve]+makeSieves :: Integer -> Integer -> Integer -> Integer -> U.Vector Word64 -> [PrimeSieve] makeSieves plim sqlim bitOff valOff cache | valOff' < sqlim = let (nc, bs) = runST $ do- cch <- unsafeThaw cache :: ST s (STUArray s Int Word64)+ cch <- U.unsafeThaw cache :: ST s (MU.MVector s Word64) bs0 <- slice cch- fcch <- unsafeFreeze cch- fbs0 <- unsafeFreeze bs0+ fcch <- U.unsafeFreeze cch+ fbs0 <- U.unsafeFreeze bs0 return (fcch, fbs0) in PS valOff bs : makeSieves plim sqlim bitOff' valOff' nc | otherwise =@@ -194,26 +193,26 @@ (nc,bs) = runST $ do cch <- growCache bitOff plim cache bs0 <- slice cch- fcch <- unsafeFreeze cch- fbs0 <- unsafeFreeze bs0+ fcch <- U.unsafeFreeze cch+ fbs0 <- U.unsafeFreeze bs0 return (fcch, fbs0) in PS valOff bs : makeSieves plim' sqlim' bitOff' valOff' nc where valOff' = valOff + intToInteger sieveRange bitOff' = bitOff + intToInteger sieveBits -slice :: STUArray s Int Word64 -> ST s (STUArray s Int Bool)+slice :: MU.MVector s Word64 -> ST s (MU.MVector s Bit) slice cache = do- hi <- snd `fmap` getBounds cache- sieve <- newArray (0,lastIndex) True+ let hi = MU.length cache - 1+ sieve <- MU.replicate (lastIndex + 1) (Bit True) let treat pr | hi < pr = return sieve | otherwise = do- w <- unsafeRead cache pr+ w <- MU.unsafeRead cache pr if w /= 0- then unsafeWrite cache pr (w-1)+ then MU.unsafeWrite cache pr (w-1) else do- ixes <- unsafeRead cache (pr+1)+ ixes <- MU.unsafeRead cache (pr+1) let !stj = word64ToInt ixes .&. iXJMASK -- position of multiple and index of cofactor !ixw = word64ToInt (ixes `shiftR` iXJBITS) -- prime data, up to 41 bits !i = ixw .&. jMASK@@ -224,19 +223,19 @@ (n, u) <- tick k o j s let !skip = intToWord64 (n `shiftR` iXBITS) !strt = intToWord64 (n .&. iXMASK)- unsafeWrite cache pr skip- unsafeWrite cache (pr+1) ((ixes .&. complement iXJMASK) .|. strt `shiftL` jBITS .|. intToWord64 u)+ MU.unsafeWrite cache pr skip+ MU.unsafeWrite cache (pr+1) ((ixes .&. complement iXJMASK) .|. strt `shiftL` jBITS .|. intToWord64 u) treat (pr+2) tick stp off j ix | lastIndex < ix = return (ix - sieveBits, j) | otherwise = do- p <- unsafeRead sieve ix- when p (unsafeWrite sieve ix False)+ Bit p <- MU.unsafeRead sieve ix+ when p (MU.unsafeWrite sieve ix (Bit False)) tick stp off ((j+1) .&. jMASK) (ix + stp*delta j + tau (off+j)) treat 0 -- | Sieve up to bound in one go.-sieveTo :: Integer -> ST s (STUArray s Int Bool)+sieveTo :: Integer -> ST s (MU.MVector s Bit) sieveTo bound = arr where (bytes,lidx) = idxPr bound@@ -247,18 +246,18 @@ (kr,r) = idxPr mxsve !svbd = 8*kr+r arr = do- ar <- newArray (0,mxidx) True+ ar <- MU.replicate (mxidx + 1) (Bit True) let start k i = 8*(k*(30*k+2*rho i) + byte i) + idx i tick stp off j ix | mxidx < ix = return () | otherwise = do- p <- unsafeRead ar ix- when p (unsafeWrite ar ix False)+ Bit p <- MU.unsafeRead ar ix+ when p (MU.unsafeWrite ar ix (Bit False)) tick stp off ((j+1) .&. jMASK) (ix + stp*delta j + tau (off+j)) sift ix | svbd < ix = return ar | otherwise = do- p <- unsafeRead ar ix+ Bit p <- MU.unsafeRead ar ix when p (do let i = ix .&. jMASK k = ix `shiftR` jBITS !off = i `shiftL` jBITS@@ -267,26 +266,26 @@ sift (ix+1) sift 0 -growCache :: Integer -> Integer -> UArray Int Word64 -> ST s (STUArray s Int Word64)+growCache :: Integer -> Integer -> U.Vector Word64 -> ST s (MU.MVector s Word64) growCache offset plim old = do- let (_,num) = bounds old+ let num = U.length old - 1 (bt,ix) = idxPr plim !start = 8*bt+ix+1 !nlim = plim+4800 sieve <- sieveTo nlim -- Implement SieveFromTo for this, it's pretty wasteful when nlim isn't- (_,hi) <- getBounds sieve -- very small anymore+ let hi = MU.length sieve - 1 more <- countFromToWd start hi sieve- new <- unsafeNewArray_ (0,num+2*more) :: ST s (STUArray s Int Word64)+ new <- MU.unsafeNew (1 + num + 2 * more) :: ST s (MU.MVector s Word64) let copy i | num < i = return () | otherwise = do- unsafeWrite new i (old `unsafeAt` i)+ MU.unsafeWrite new i (old `U.unsafeIndex` i) copy (i+1) copy 0 let fill j indx | hi < indx = return new | otherwise = do- p <- unsafeRead sieve indx+ Bit p <- MU.unsafeRead sieve indx if p then do let !i = indx .&. jMASK@@ -299,27 +298,16 @@ !strt = integerToWord64 strt1 .&. iXMASK !skip = integerToWord64 (strt1 `shiftR` iXBITS) !ixes = intToWord64 indx `shiftL` iXJBITS .|. strt `shiftL` jBITS .|. intToWord64 i- unsafeWrite new j skip- unsafeWrite new (j+1) ixes+ MU.unsafeWrite new j skip+ MU.unsafeWrite new (j+1) ixes fill (j+2) (indx+1) else fill j (indx+1) fill (num+1) start --- Danger: relies on start and end being the first resp. last--- index in a Word--- Do not use except in growCache and psieveFrom {-# INLINE countFromToWd #-}-countFromToWd :: Int -> Int -> STUArray s Int Bool -> ST s Int-countFromToWd start end ba = do- wa <- (castSTUArray :: STUArray s Int Bool -> ST s (STUArray s Int Word)) ba- let !sb = start `shiftR` wSHFT- !eb = end `shiftR` wSHFT- count !acc i- | eb < i = return acc- | otherwise = do- w <- unsafeRead wa i- count (acc + popCount w) (i+1)- count 0 sb+countFromToWd :: Int -> Int -> MU.MVector s Bit -> ST s Int+countFromToWd start end =+ fmap countBits . U.unsafeFreeze . MU.slice start (end - start + 1) -- | @'psieveFrom' n@ creates the list of 'PrimeSieve's starting roughly -- at @n@. Due to the organisation of the sieve, the list may contain@@ -338,15 +326,15 @@ plim0 = integerSquareRoot end1 plim = plim0 + 4801 - (plim0 `rem` 4800) sqlim = plim*plim- cache = runSTUArray $ do+ cache = runST $ do sieve <- sieveTo plim- (lo,hi) <- getBounds sieve+ let (lo,hi) = (0, MU.length sieve - 1) pct <- countFromToWd lo hi sieve- new <- unsafeNewArray_ (0,2*pct-1) :: ST s (STUArray s Int Word64)+ new <- MU.unsafeNew (2 * pct) :: ST s (MU.MVector s Word64) let fill j indx | hi < indx = return new | otherwise = do- isPr <- unsafeRead sieve indx+ Bit isPr <- MU.unsafeRead sieve indx if isPr then do let !i = indx .&. jMASK@@ -376,26 +364,26 @@ !strt = integerToWord64 strt2 .&. iXMASK !skip = integerToWord64 (strt2 `shiftR` iXBITS) !ixes = intToWord64 indx `shiftL` iXJBITS .|. strt `shiftL` jBITS .|. intToWord64 r2- unsafeWrite new j skip- unsafeWrite new (j+1) ixes+ MU.unsafeWrite new j skip+ MU.unsafeWrite new (j+1) ixes fill (j+2) (indx+1) else fill j (indx+1)- fill 0 0-+ vec <- fill 0 0+ U.unsafeFreeze vec {-# INLINE delta #-} delta :: Int -> Int-delta = unsafeAt deltas+delta = U.unsafeIndex deltas -deltas :: UArray Int Int-deltas = listArray (0,7) [4,2,4,2,4,6,2,6]+deltas :: U.Vector Int+deltas = U.fromList [4,2,4,2,4,6,2,6] {-# INLINE tau #-} tau :: Int -> Int-tau = unsafeAt taus+tau = U.unsafeIndex taus -taus :: UArray Int Int-taus = listArray (0,63)+taus :: U.Vector Int+taus = U.fromList [ 7, 4, 7, 4, 7, 12, 3, 12 , 12, 6, 11, 6, 12, 18, 5, 18 , 14, 7, 13, 7, 14, 21, 7, 21@@ -408,28 +396,28 @@ {-# INLINE byte #-} byte :: Int -> Int-byte = unsafeAt startByte+byte = U.unsafeIndex startByte -startByte :: UArray Int Int-startByte = listArray (0,7) [1,3,5,9,11,17,27,31]+startByte :: U.Vector Int+startByte = U.fromList [1,3,5,9,11,17,27,31] {-# INLINE idx #-} idx :: Int -> Int-idx = unsafeAt startIdx+idx = U.unsafeIndex startIdx -startIdx :: UArray Int Int-startIdx = listArray (0,7) [4,7,4,4,7,4,7,7]+startIdx :: U.Vector Int+startIdx = U.fromList [4,7,4,4,7,4,7,7] {-# INLINE mu #-} mu :: Int -> Int-mu = unsafeAt mArr+mu = U.unsafeIndex mArr {-# INLINE nu #-} nu :: Int -> Int-nu = unsafeAt nArr+nu = U.unsafeIndex nArr -mArr :: UArray Int Int-mArr = listArray (0,63)+mArr :: U.Vector Int+mArr = U.fromList [ 1, 2, 2, 3, 4, 5, 6, 7 , 2, 3, 4, 6, 6, 8, 10, 11 , 2, 4, 5, 7, 8, 9, 12, 13@@ -440,8 +428,8 @@ , 7, 11, 13, 17, 19, 23, 29, 31 ] -nArr :: UArray Int Int-nArr = listArray (0,63)+nArr :: U.Vector Int+nArr = U.fromList [ 4, 3, 7, 6, 2, 1, 5, 0 , 3, 7, 5, 0, 6, 2, 4, 1 , 7, 5, 4, 1, 0, 6, 3, 2
Math/NumberTheory/Primes/Sieve/Indexing.hs view
@@ -12,8 +12,8 @@ , rho ) where -import Data.Array.Base import Data.Bits+import qualified Data.Vector.Unboxed as U {-# INLINE idxPr #-} idxPr :: Integral a => a -> (Int,Int)@@ -37,7 +37,7 @@ {-# INLINE rho #-} rho :: Int -> Int-rho = unsafeAt residues+rho = U.unsafeIndex residues -residues :: UArray Int Int-residues = listArray (0,7) [7,11,13,17,19,23,29,31]+residues :: U.Vector Int+residues = U.fromList [7,11,13,17,19,23,29,31]
Math/NumberTheory/Recurrences/Bilinear.hs view
@@ -343,5 +343,4 @@ helperForBEEP g xs = Inf.map (f . g) stirling2 where f = sum . zipWith4 (\sgn fact x stir -> sgn * fact * x * stir) (cycle [1, -1]) (Inf.toList factorial) (Inf.toList xs)-{-# SPECIALIZE helperForBEEP :: (forall b. [b] -> [b]) -> Infinite (Ratio Int) -> Infinite (Ratio Int) #-}-{-# SPECIALIZE helperForBEEP :: (forall b. [b] -> [b]) -> Infinite (Rational) -> Infinite (Rational) #-}+{-# INLINABLE helperForBEEP #-}
Math/NumberTheory/Utils/DirichletSeries.hs view
@@ -22,14 +22,22 @@ , timesAndCrop ) where -import Prelude hiding (filter, last, rem, quot, snd, lookup)-import Data.Coerce-import Data.Euclidean+import Prelude (Integer, Num, (*))+import Data.Bool (Bool)+import Data.Coerce (coerce)+import Data.Euclidean (Euclidean, quot, divide)+import Data.Function (($), (.))+import Data.Int (Int)+import Data.List (takeWhile) import Data.Map (Map) import qualified Data.Map.Strict as M-import Data.Maybe+import Data.Maybe (isJust)+import Data.Ord (Ord, (<=)) import Data.Semiring (Semiring(..))-import Numeric.Natural+import Data.Tuple (fst)+import Data.Word (Word)+import Numeric.Natural (Natural)+import Text.Show (Show) -- Sparse Dirichlet series are represented by an ascending list of pairs. -- For instance, [(a, b), (c, d)] stands for 1 + b/a^s + d/c^s.
arithmoi.cabal view
@@ -1,5 +1,5 @@ name: arithmoi-version: 0.13.1.0+version: 0.13.2.0 cabal-version: 2.0 build-type: Simple license: MIT@@ -17,7 +17,7 @@ powers (integer roots and tests, modular exponentiation). category: Math, Algorithms, Number Theory author: Andrew Lelechenko, Daniel Fischer-tested-with: GHC ==9.0.2 GHC ==9.2.8 GHC ==9.4.8 GHC ==9.6.7 GHC ==9.8.4 GHC ==9.10.1 GHC ==9.12.2+tested-with: GHC ==9.0.2 GHC ==9.2.8 GHC ==9.4.8 GHC ==9.6.7 GHC ==9.8.4 GHC ==9.10.2 GHC ==9.12.2 GHC ==9.14.1 extra-doc-files: changelog.md @@ -28,7 +28,7 @@ library build-depends: base >=4.15 && <5,- array >=0.5 && <0.6,+ bitvec >=1.1 && <1.2, containers >=0.5.11 && <0.9, chimera >=0.3 && <0.5, constraints <0.15,@@ -38,7 +38,7 @@ infinite-list <0.2, integer-logarithms >=1.0 && <1.1, integer-roots >=1.0 && <1.1,- mod <0.3,+ mod >=0.2.1 && <0.3, random >=1.0 && <1.4, transformers >=0.4 && <0.7, semirings >=0.5.2 && <0.8,@@ -113,9 +113,8 @@ infinite-list, integer-roots >=1.0, mod,- QuickCheck >=2.10 && <2.16,+ QuickCheck >=2.10 && <2.17, quickcheck-classes >=0.6.3 && <0.7,- random, semirings >=0.2, smallcheck >=1.2 && <1.3, tasty >=0.10 && <1.6,@@ -123,7 +122,6 @@ tasty-quickcheck >=0.9 && <0.12, tasty-rerun >=1.1.17 && <1.2, tasty-smallcheck >=0.8 && <0.9,- transformers >=0.5, vector other-modules: Math.NumberTheory.ArithmeticFunctionsTests@@ -176,10 +174,8 @@ build-depends: base, arithmoi,- array, constraints, containers,- deepseq, infinite-list, integer-logarithms, mod,
benchmark/Math/NumberTheory/SequenceBench.hs view
@@ -6,8 +6,8 @@ import Test.Tasty.Bench -import Data.Array.Unboxed import Data.Bits+import qualified Data.Vector.Unboxed as U import Math.NumberTheory.Primes (Prime(..), nextPrime, precPrime) import Math.NumberTheory.Primes.Testing@@ -43,10 +43,10 @@ -- Utils copypasted from internal modules rho :: Int -> Int-rho i = residues ! i+rho i = residues U.! i -residues :: UArray Int Int-residues = listArray (0,7) [7,11,13,17,19,23,29,31]+residues :: U.Vector Int+residues = U.fromList [7,11,13,17,19,23,29,31] toIdx :: Integral a => a -> Int toIdx n = 8*fromIntegral q+r2
changelog.md view
@@ -1,5 +1,11 @@ # Changelog +## 0.13.2.0++### Changed++* Migrate implementation of prime sieves from `Array Bool` to `Vector Bit`.+ ## 0.13.1.0 ### Fixed
test-suite/Math/NumberTheory/DirichletCharactersTests.hs view
@@ -141,9 +141,8 @@ -- | Primitive checker is correct (in both directions) primitiveCheck :: forall n. KnownNat n => DirichletCharacter n -> Bool primitiveCheck chi = isJust (isPrimitive chi) == isPrimitive'- where isPrimitive' = all testModulus possibleModuli+ where isPrimitive' = all (\m -> m == n || testModulus m) (divisorsList n) n = fromIntegral (natVal @n Proxy) :: Int- possibleModuli = init (divisorsList n) table = evalAll chi testModulus d = not $ null [a | a <- [1..n-1], gcd a n == 1, a `mod` d == 1 `mod` d, table V.! a /= mempty]
test-suite/Math/NumberTheory/Moduli/ClassTests.hs view
@@ -35,7 +35,7 @@ invertModProperty (AnySign k) (Positive m) = case invertMod k m of Nothing -> k `rem` m == 0 || gcd k m > 1 Just InfMod{} -> False- Just (SomeMod inv) -> gcd k m == 1 && k * getVal inv `mod` m == 1+ Just (SomeMod inv) -> gcd k m == 1 && k * getVal inv `mod` m == 1 `mod` m -- | Check that 'powerMod' is multiplicative by first argument. powerModProperty2 :: (Integral a) => NonNegative a -> AnySign Integer -> AnySign Integer -> Positive Integer -> Bool