diff --git a/Math/NumberTheory/DirichletCharacters.hs b/Math/NumberTheory/DirichletCharacters.hs --- a/Math/NumberTheory/DirichletCharacters.hs +++ b/Math/NumberTheory/DirichletCharacters.hs @@ -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 diff --git a/Math/NumberTheory/Moduli/JacobiSymbol.hs b/Math/NumberTheory/Moduli/JacobiSymbol.hs --- a/Math/NumberTheory/Moduli/JacobiSymbol.hs +++ b/Math/NumberTheory/Moduli/JacobiSymbol.hs @@ -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 #-} + -- | 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 diff --git a/Math/NumberTheory/Moduli/Singleton.hs b/Math/NumberTheory/Moduli/Singleton.hs --- a/Math/NumberTheory/Moduli/Singleton.hs +++ b/Math/NumberTheory/Moduli/Singleton.hs @@ -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 diff --git a/Math/NumberTheory/Moduli/SomeMod.hs b/Math/NumberTheory/Moduli/SomeMod.hs --- a/Math/NumberTheory/Moduli/SomeMod.hs +++ b/Math/NumberTheory/Moduli/SomeMod.hs @@ -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'. diff --git a/Math/NumberTheory/Primes/Counting/Impl.hs b/Math/NumberTheory/Primes/Counting/Impl.hs --- a/Math/NumberTheory/Primes/Counting/Impl.hs +++ b/Math/NumberTheory/Primes/Counting/Impl.hs @@ -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) diff --git a/Math/NumberTheory/Primes/Factorisation/Montgomery.hs b/Math/NumberTheory/Primes/Factorisation/Montgomery.hs --- a/Math/NumberTheory/Primes/Factorisation/Montgomery.hs +++ b/Math/NumberTheory/Primes/Factorisation/Montgomery.hs @@ -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. diff --git a/Math/NumberTheory/Primes/Sieve/Eratosthenes.hs b/Math/NumberTheory/Primes/Sieve/Eratosthenes.hs --- a/Math/NumberTheory/Primes/Sieve/Eratosthenes.hs +++ b/Math/NumberTheory/Primes/Sieve/Eratosthenes.hs @@ -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 diff --git a/Math/NumberTheory/Primes/Sieve/Indexing.hs b/Math/NumberTheory/Primes/Sieve/Indexing.hs --- a/Math/NumberTheory/Primes/Sieve/Indexing.hs +++ b/Math/NumberTheory/Primes/Sieve/Indexing.hs @@ -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] diff --git a/Math/NumberTheory/Recurrences/Bilinear.hs b/Math/NumberTheory/Recurrences/Bilinear.hs --- a/Math/NumberTheory/Recurrences/Bilinear.hs +++ b/Math/NumberTheory/Recurrences/Bilinear.hs @@ -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 #-} diff --git a/Math/NumberTheory/Utils/DirichletSeries.hs b/Math/NumberTheory/Utils/DirichletSeries.hs --- a/Math/NumberTheory/Utils/DirichletSeries.hs +++ b/Math/NumberTheory/Utils/DirichletSeries.hs @@ -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. diff --git a/arithmoi.cabal b/arithmoi.cabal --- a/arithmoi.cabal +++ b/arithmoi.cabal @@ -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, diff --git a/benchmark/Math/NumberTheory/SequenceBench.hs b/benchmark/Math/NumberTheory/SequenceBench.hs --- a/benchmark/Math/NumberTheory/SequenceBench.hs +++ b/benchmark/Math/NumberTheory/SequenceBench.hs @@ -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 diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -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 diff --git a/test-suite/Math/NumberTheory/DirichletCharactersTests.hs b/test-suite/Math/NumberTheory/DirichletCharactersTests.hs --- a/test-suite/Math/NumberTheory/DirichletCharactersTests.hs +++ b/test-suite/Math/NumberTheory/DirichletCharactersTests.hs @@ -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] diff --git a/test-suite/Math/NumberTheory/Moduli/ClassTests.hs b/test-suite/Math/NumberTheory/Moduli/ClassTests.hs --- a/test-suite/Math/NumberTheory/Moduli/ClassTests.hs +++ b/test-suite/Math/NumberTheory/Moduli/ClassTests.hs @@ -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