diff --git a/Math/NumberTheory/ArithmeticFunctions/Inverse.hs b/Math/NumberTheory/ArithmeticFunctions/Inverse.hs --- a/Math/NumberTheory/ArithmeticFunctions/Inverse.hs +++ b/Math/NumberTheory/ArithmeticFunctions/Inverse.hs @@ -153,7 +153,7 @@ pksLarge :: Map (Prime a) (Set Word) pksLarge = M.unionsWith (<>) - [ maybe mempty (flip M.singleton (S.singleton e)) (isPrime p) + [ maybe mempty (`M.singleton` S.singleton e) (isPrime p) | d <- divs , e <- [1 .. intToWord (quot (integerLogBase (toInteger lim) (toInteger d)) (wordToInt k)) ] , let p = integerRoot (e * k) (d - 1) @@ -206,7 +206,7 @@ -- ^ Semiring element, representing preimages invertFunction point f invF n = DS.lookup n - $ foldl' (\ds b -> uncurry processBatch b ds) (DS.fromDistinctAscList []) batches + $ foldl' (flip (uncurry processBatch)) (DS.fromDistinctAscList []) batches where factors = factorise n batches = strategy f factors $ invF factors @@ -326,6 +326,7 @@ -- -- >>> import qualified Data.Set as S -- >>> import Data.Semigroup +-- >>> :set -XFlexibleContexts -- >>> S.mapMonotonic getProduct (inverseSigma (S.singleton . Product) 120) -- fromList [54,56,87,95] -- @@ -366,8 +367,9 @@ -- -- >>> import qualified Data.Set as S -- >>> import Data.Semigroup +-- >>> :set -XFlexibleContexts -- >>> S.mapMonotonic getProduct (inverseSigmaK 2 (S.singleton . Product) 850) --- fromList [24, 26] +-- fromList [24,26] -- -- Similarly to 'inverseSigma', it is possible to count and sum preimages, or -- get the maximum/minimum preimage. diff --git a/Math/NumberTheory/ArithmeticFunctions/Mertens.hs b/Math/NumberTheory/ArithmeticFunctions/Mertens.hs --- a/Math/NumberTheory/ArithmeticFunctions/Mertens.hs +++ b/Math/NumberTheory/ArithmeticFunctions/Mertens.hs @@ -7,7 +7,6 @@ -- Values of . -- -{-# LANGUAGE CPP #-} {-# LANGUAGE LambdaCase #-} module Math.NumberTheory.ArithmeticFunctions.Mertens @@ -18,6 +17,7 @@ import Math.NumberTheory.Roots import Math.NumberTheory.ArithmeticFunctions.Moebius +import Math.NumberTheory.Utils.FromIntegral -- | Compute individual values of Mertens function in O(n^(2/3)) time and space. -- @@ -30,14 +30,14 @@ mertens 1 = 1 mertens x = sumMultMoebius lookupMus (\n -> sfunc (x `quot` n)) [1 .. x `quot` u] where - u = (integerSquareRoot x + 1) `max` ((integerCubeRoot x) ^ (2 :: Word) `quot` 2) + u = (integerSquareRoot x + 1) `max` (integerCubeRoot x ^ (2 :: Word) `quot` 2) sfunc :: Word -> Int sfunc y = 1 - - sum [ U.unsafeIndex mes (fromIntegral $ y `quot` n) | n <- [y `quot` u + 1 .. kappa] ] - + fromIntegral kappa * U.unsafeIndex mes (fromIntegral nu) - - sumMultMoebius lookupMus (\n -> fromIntegral $ y `quot` n) [1 .. nu] + - sum [ U.unsafeIndex mes (wordToInt $ y `quot` n) | n <- [y `quot` u + 1 .. kappa] ] + + wordToInt kappa * U.unsafeIndex mes (wordToInt nu) + - sumMultMoebius lookupMus (\n -> wordToInt $ y `quot` n) [1 .. nu] where nu = integerSquareRoot y kappa = y `quot` (nu + 1) @@ -51,7 +51,7 @@ mus = sieveBlockMoebius 1 cacheSize lookupMus :: Word -> Moebius - lookupMus i = U.unsafeIndex mus (fromIntegral (i - 1)) + lookupMus i = U.unsafeIndex mus (wordToInt (i - 1)) -- 0-based index mes :: U.Vector Int diff --git a/Math/NumberTheory/ArithmeticFunctions/Moebius.hs b/Math/NumberTheory/ArithmeticFunctions/Moebius.hs --- a/Math/NumberTheory/ArithmeticFunctions/Moebius.hs +++ b/Math/NumberTheory/ArithmeticFunctions/Moebius.hs @@ -19,7 +19,7 @@ , sieveBlockMoebius ) where -import Control.Monad (forM_, liftM) +import Control.Monad (forM_) import Control.Monad.ST (runST) import Data.Bits import Data.Int @@ -38,7 +38,7 @@ import Math.NumberTheory.Roots (integerSquareRoot) import Math.NumberTheory.Primes -import Math.NumberTheory.Utils.FromIntegral (wordToInt) +import Math.NumberTheory.Utils.FromIntegral import Math.NumberTheory.Logarithms @@ -54,11 +54,11 @@ runMoebius m = fromInteger (S# (dataToTag# m -# 1#)) fromMoebius :: Moebius -> Int8 -fromMoebius m = fromIntegral $ I# (dataToTag# m) +fromMoebius m = intToInt8 $ I# (dataToTag# m) {-# INLINE fromMoebius #-} toMoebius :: Int8 -> Moebius -toMoebius i = let !(I# i#) = fromIntegral i in tagToEnum# i# +toMoebius i = let !(I# i#) = int8ToInt i in tagToEnum# i# {-# INLINE toMoebius #-} newtype instance U.MVector s Moebius = MV_Moebius (P.MVector s Int8) @@ -82,16 +82,16 @@ basicLength (MV_Moebius v) = M.basicLength v basicUnsafeSlice i n (MV_Moebius v) = MV_Moebius $ M.basicUnsafeSlice i n v basicOverlaps (MV_Moebius v1) (MV_Moebius v2) = M.basicOverlaps v1 v2 - basicUnsafeNew n = MV_Moebius `liftM` M.basicUnsafeNew n + basicUnsafeNew n = MV_Moebius <$> M.basicUnsafeNew n basicInitialize (MV_Moebius v) = M.basicInitialize v - basicUnsafeReplicate n x = MV_Moebius `liftM` M.basicUnsafeReplicate n (fromMoebius x) - basicUnsafeRead (MV_Moebius v) i = toMoebius `liftM` M.basicUnsafeRead v i + basicUnsafeReplicate n x = MV_Moebius <$> M.basicUnsafeReplicate n (fromMoebius x) + basicUnsafeRead (MV_Moebius v) i = toMoebius <$> M.basicUnsafeRead v i basicUnsafeWrite (MV_Moebius v) i x = M.basicUnsafeWrite v i (fromMoebius x) basicClear (MV_Moebius v) = M.basicClear v basicSet (MV_Moebius v) x = M.basicSet v (fromMoebius x) basicUnsafeCopy (MV_Moebius v1) (MV_Moebius v2) = M.basicUnsafeCopy v1 v2 basicUnsafeMove (MV_Moebius v1) (MV_Moebius v2) = M.basicUnsafeMove v1 v2 - basicUnsafeGrow (MV_Moebius v) n = MV_Moebius `liftM` M.basicUnsafeGrow v n + basicUnsafeGrow (MV_Moebius v) n = MV_Moebius <$> M.basicUnsafeGrow v n instance G.Vector U.Vector Moebius where {-# INLINE basicUnsafeFreeze #-} @@ -100,11 +100,11 @@ {-# INLINE basicUnsafeSlice #-} {-# INLINE basicUnsafeIndexM #-} {-# INLINE elemseq #-} - basicUnsafeFreeze (MV_Moebius v) = V_Moebius `liftM` G.basicUnsafeFreeze v - basicUnsafeThaw (V_Moebius v) = MV_Moebius `liftM` G.basicUnsafeThaw v + basicUnsafeFreeze (MV_Moebius v) = V_Moebius <$> G.basicUnsafeFreeze v + basicUnsafeThaw (V_Moebius v) = MV_Moebius <$> G.basicUnsafeThaw v basicLength (V_Moebius v) = G.basicLength v basicUnsafeSlice i n (V_Moebius v) = V_Moebius $ G.basicUnsafeSlice i n v - basicUnsafeIndexM (V_Moebius v) i = toMoebius `liftM` G.basicUnsafeIndexM v i + basicUnsafeIndexM (V_Moebius v) i = toMoebius <$> G.basicUnsafeIndexM v i basicUnsafeCopy (MV_Moebius mv) (V_Moebius v) = G.basicUnsafeCopy mv v elemseq _ = seq @@ -138,12 +138,12 @@ let offset = negate lowIndex `mod` p offset2 = negate lowIndex `mod` (p * p) l :: Word8 - l = fromIntegral $ intLog2 p .|. 1 - forM_ [offset, offset + p .. len - 1] $ \ix -> do - MU.unsafeModify as (\y -> y + l) ix - forM_ [offset2, offset2 + p * p .. len - 1] $ \ix -> do + l = intToWord8 $ intLog2 p .|. 1 + forM_ [offset, offset + p .. len - 1] $ + MU.unsafeModify as (+ l) + forM_ [offset2, offset2 + p * p .. len - 1] $ \ix -> MU.unsafeWrite as ix 0 - forM_ [0 .. len - 1] $ \ix -> do + forM_ [0 .. len - 1] $ \ix -> MU.unsafeModify as (mapper ix) ix U.unsafeFreeze as @@ -166,7 +166,7 @@ mapper ix val | val .&. 0x80 == 0x00 = 1 - | fromIntegral (val .&. 0x7F) < intLog2 (ix + lowIndex) - 5 + | word8ToInt (val .&. 0x7F) < intLog2 (ix + lowIndex) - 5 - (if ix + lowIndex >= 0x100000 then 2 else 0) - (if ix + lowIndex >= 0x10000000 then 1 else 0) = (val .&. 1) `shiftL` 1 diff --git a/Math/NumberTheory/ArithmeticFunctions/NFreedom.hs b/Math/NumberTheory/ArithmeticFunctions/NFreedom.hs --- a/Math/NumberTheory/ArithmeticFunctions/NFreedom.hs +++ b/Math/NumberTheory/ArithmeticFunctions/NFreedom.hs @@ -25,7 +25,7 @@ import Math.NumberTheory.Roots import Math.NumberTheory.Primes -import Math.NumberTheory.Utils.FromIntegral (wordToInt) +import Math.NumberTheory.Utils.FromIntegral -- | Evaluate the `Math.NumberTheory.ArithmeticFunctions.isNFree` function over a block. -- Value at @0@, if zero falls into block, is undefined. @@ -71,7 +71,7 @@ -- by @nFrees@, see the comments in it. indices :: [a] indices = [offset, offset + pPow .. len - 1] - forM_ indices $ \ix -> do + forM_ indices $ \ix -> MU.write as (fromIntegral ix) False U.freeze as @@ -97,7 +97,7 @@ -- ^ Generated infinite list of @n@-free numbers. nFrees 0 = [1] nFrees 1 = [1] -nFrees n = concatMap (\(lo, len) -> nFreesBlock n lo len) $ zip bounds strides +nFrees n = concatMap (uncurry (nFreesBlock n)) $ zip bounds strides where -- The 56th element of @iterate (2 *) 256@ is @2^64 :: Word == 0@, so to -- avoid overflow only the first 55 elements of this list are used. @@ -105,7 +105,7 @@ -- this value is just repeated. This means after a few dozen iterations, -- the sieve will stop increasing in size. strides :: [Word] - strides = take 55 (iterate (2 *) 256) ++ repeat (fromIntegral (maxBound :: Int)) + strides = take 55 (iterate (2 *) 256) ++ repeat (intToWord (maxBound :: Int)) -- Infinite list of lower bounds at which @sieveBlockNFree@ will be -- applied. This has type @Integral a => a@ instead of @Word@ because diff --git a/Math/NumberTheory/ArithmeticFunctions/SieveBlock.hs b/Math/NumberTheory/ArithmeticFunctions/SieveBlock.hs --- a/Math/NumberTheory/ArithmeticFunctions/SieveBlock.hs +++ b/Math/NumberTheory/ArithmeticFunctions/SieveBlock.hs @@ -11,7 +11,6 @@ {-# LANGUAGE BangPatterns #-} {-# LANGUAGE MagicHash #-} {-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE UnboxedTuples #-} module Math.NumberTheory.ArithmeticFunctions.SieveBlock ( runFunctionOverBlock @@ -109,12 +108,12 @@ -- -- For example, following code lists smallest prime factors: -- --- >>> sieveBlock (SieveBlockConfig maxBound (\p _ -> unPrime p) min) 2 10 +-- >>> sieveBlock (SieveBlockConfig maxBound (\p _ -> unPrime p) min) 2 10 :: Data.Vector.Vector Word -- [2,3,2,5,2,7,2,3,2,11] -- -- And this is how to factorise all numbers in a block: -- --- >>> sieveBlock (SieveBlockConfig [] (\p k -> [(unPrime p, k)]) (++)) 2 10 +-- >>> sieveBlock (SieveBlockConfig [] (\p k -> [(unPrime p, k)]) (++)) 2 10 :: Data.Vector.Vector [(Word, Word)] -- [[(2,1)],[(3,1)],[(2,2)],[(5,1)],[(2,1),(3,1)],[(7,1)],[(2,3)],[(3,2)],[(2,1),(5,1)],[(11,1)]] sieveBlock :: forall v a. @@ -168,8 +167,8 @@ (q, r) = np `quotRem` p if r /= 0 then do - MU.unsafeModify as (\x -> x * p') (I# ix#) - MG.unsafeModify bs (\y -> y `append` f0) (I# ix#) + MU.unsafeModify as (* p') (I# ix#) + MG.unsafeModify bs (`append` f0) (I# ix#) else do let pow = highestPowerDividing p q MU.unsafeModify as (\x -> x * p' ^ (pow + 2)) (I# ix#) diff --git a/Math/NumberTheory/ArithmeticFunctions/Standard.hs b/Math/NumberTheory/ArithmeticFunctions/Standard.hs --- a/Math/NumberTheory/ArithmeticFunctions/Standard.hs +++ b/Math/NumberTheory/ArithmeticFunctions/Standard.hs @@ -144,8 +144,8 @@ -- > sigmaA 0 = tauA sigmaA :: (Integral n, Num a, GcdDomain a) => Word -> ArithmeticFunction n a sigmaA 0 = tauA -sigmaA 1 = multiplicative $ sigmaHelper . fromIntegral . unPrime -sigmaA a = multiplicative $ sigmaHelper . (^ wordToInt a) . fromIntegral . unPrime +sigmaA 1 = multiplicative $ sigmaHelper . fromIntegral' . unPrime +sigmaA a = multiplicative $ sigmaHelper . (^ wordToInt a) . fromIntegral' . unPrime {-# INLINABLE sigmaA #-} sigmaHelper :: (Num a, GcdDomain a) => a -> Word -> a @@ -202,7 +202,7 @@ tp = ramanujanHelper p 1 paPowers = iterate (* (-pa)) 1 binomials = scanl (\acc j -> acc * (k' - 2 * j) * (k' - 2 * j - 1) `quot` (k' - j) `quot` (j + 1)) 1 [0 .. k' `quot` 2 - 1] - k' = fromIntegral k + k' = wordToInteger k tpPowers = reverse $ take (length binomials) $ iterate (* tp^(2::Int)) (if even k then 1 else tp) {-# INLINE ramanujanHelper #-} diff --git a/Math/NumberTheory/Diophantine.hs b/Math/NumberTheory/Diophantine.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Diophantine.hs @@ -0,0 +1,66 @@ +-- Module for Diophantine Equations and related functions + +module Math.NumberTheory.Diophantine + ( cornacchiaPrimitive + , cornacchia + ) +where + +import Math.NumberTheory.Moduli.Sqrt ( sqrtsModFactorisation ) +import Math.NumberTheory.Primes ( factorise + , unPrime + , UniqueFactorisation + ) +import Math.NumberTheory.Roots ( integerSquareRoot ) +import Math.NumberTheory.Utils.FromIntegral + +-- | See `cornacchiaPrimitive`, this is the internal algorithm implementation +-- | as described at https://en.wikipedia.org/wiki/Cornacchia%27s_algorithm +cornacchiaPrimitive' :: Integer -> Integer -> [(Integer, Integer)] +cornacchiaPrimitive' d m = concatMap + (findSolution . head . dropWhile (\r -> r * r >= m) . gcdSeq m) + roots + where + roots = filter (<= m `div` 2) $ sqrtsModFactorisation (m - d) (factorise m) + gcdSeq a b = a : gcdSeq b (mod a b) + -- If s = sqrt((m - r*r) / d) is an integer then (r, s) is a solution + findSolution r = [ (r, s) | rem1 == 0 && s * s == s2 ] + where + (s2, rem1) = divMod (m - r * r) d + s = integerSquareRoot s2 + +-- | Finds all primitive solutions (x,y) to the diophantine equation +-- | x^2 + d*y^2 = m +-- | when 1 <= d < m and gcd(d,m)=1 +-- | Given m is square free these are all the positive integer solutions +cornacchiaPrimitive :: Integer -> Integer -> [(Integer, Integer)] +cornacchiaPrimitive d m + | not (1 <= d && d < m) = error "precondition failed: 1 <= d < m" + | gcd d m /= 1 = error "precondition failed: d and m coprime" + | + -- If d=1 then the algorithm doesn't generate symmetrical pairs + d == 1 = concatMap genPairs solutions + | otherwise = solutions + where + solutions = cornacchiaPrimitive' d m + genPairs (x, y) = if x == y then [(x, y)] else [(x, y), (y, x)] + +-- Find numbers whose square is a factor of the input +squareFactors :: UniqueFactorisation a => a -> [a] +squareFactors = foldl squareProducts [1] . factorise + where + squareProducts acc f = [ a * b | a <- acc, b <- squarePowers f ] + squarePowers (p, a) = map (unPrime p ^) [0 .. wordToInt a `div` 2] + +-- | Finds all positive integer solutions (x,y) to the +-- | diophantine equation: +-- | x^2 + d*y^2 = m +-- | when 1 <= d < m and gcd(d,m)=1 +cornacchia :: Integer -> Integer -> [(Integer, Integer)] +cornacchia d m + | not (1 <= d && d < m) = error "precondition failed: 1 <= d < m" + | gcd d m /= 1 = error "precondition failed: d and m coprime" + | otherwise = concatMap solve $ filter ((> d) . snd) candidates + where + candidates = map (\sf -> (sf, m `div` (sf * sf))) (squareFactors m) + solve (sf, m') = map (\(x, y) -> (x * sf, y * sf)) (cornacchiaPrimitive d m') diff --git a/Math/NumberTheory/DirichletCharacters.hs b/Math/NumberTheory/DirichletCharacters.hs --- a/Math/NumberTheory/DirichletCharacters.hs +++ b/Math/NumberTheory/DirichletCharacters.hs @@ -63,11 +63,13 @@ import Control.Applicative (liftA2) #endif import Data.Bits (Bits(..)) +import Data.Constraint import Data.Foldable (for_) import Data.Functor.Identity (Identity(..)) import Data.Kind import Data.List (mapAccumL, foldl', sort, find, unfoldr) import Data.Maybe (mapMaybe, fromJust, fromMaybe) +import Data.Mod #if MIN_VERSION_base(4,12,0) import Data.Monoid (Ap(..)) #endif @@ -77,20 +79,18 @@ import qualified Data.Vector as V import qualified Data.Vector.Mutable as MV import Data.Vector (Vector, (!)) -import GHC.TypeNats (Nat, SomeNat(..), natVal, someNatVal) +import GHC.TypeNats (KnownNat, Nat, SomeNat(..), natVal, someNatVal) import Numeric.Natural (Natural) import Math.NumberTheory.ArithmeticFunctions (totient) import Math.NumberTheory.Moduli.Chinese -import Math.NumberTheory.Moduli.Class (KnownNat, Mod, getVal) -import Math.NumberTheory.Moduli.Internal (isPrimitiveRoot', discreteLogarithmPP) -import Math.NumberTheory.Moduli.Multiplicative (MultMod(..), isMultElement) -import Math.NumberTheory.Moduli.Singleton (Some(..), cyclicGroupFromFactors) -import Math.NumberTheory.Powers.Modular (powMod) -import Math.NumberTheory.Primes (Prime(..), UniqueFactorisation, factorise, nextPrime) +import Math.NumberTheory.Moduli.Internal (discreteLogarithmPP) +import Math.NumberTheory.Moduli.Multiplicative +import Math.NumberTheory.Moduli.Singleton +import Math.NumberTheory.Primes import Math.NumberTheory.RootsOfUnity -import Math.NumberTheory.Utils.FromIntegral (wordToInt) import Math.NumberTheory.Utils +import Math.NumberTheory.Utils.FromIntegral -- | A Dirichlet character mod \(n\) is a group homomorphism from \((\mathbb{Z}/n\mathbb{Z})^*\) -- to \(\mathbb{C}^*\), represented abstractly by `DirichletCharacter`. In particular, they take @@ -138,29 +138,33 @@ Two == Two = True _ == _ = False --- | For primes, define the canonical primitive root as the smallest such. For prime powers \(p^k\), --- either the smallest primitive root \(g\) mod \(p\) works, or \(g+p\) works. -generator :: (Integral a, UniqueFactorisation a) => Prime a -> Word -> a -generator p k - | k == 1 = modP - | otherwise = if powMod modP (p'-1) (p'*p') == 1 then modP + p' else modP - where p' = unPrime p - modP = case cyclicGroupFromFactors [(p,k)] of - Just (Some cg) -> head $ filter (isPrimitiveRoot' cg) [2..p'-1] - _ -> error "illegal" +-- | For primes, define the canonical primitive root as the smallest such. +generator :: Prime Natural -> Word -> Natural +generator p k = case cyclicGroupFromFactors [(p, k)] of + Nothing -> error "illegal" + Just (Some cg) + | Sub Dict <- proofFromCyclicGroup cg -> + unMod $ multElement $ unPrimitiveRoot $ head $ + mapMaybe (isPrimitiveRoot cg) [2..maxBound] -- | Implement the function \(\lambda\) from page 5 of -- https://www2.eecs.berkeley.edu/Pubs/TechRpts/1984/CSD-84-186.pdf lambda :: Integer -> Int -> Integer -lambda x e = ((powMod x (2*modulus) largeMod - 1) `shiftR` (e+1)) .&. (modulus - 1) - where modulus = bit (e-2) - largeMod = bit (2*e - 1) +lambda x e = ((xPower - 1) `shiftR` (e+1)) .&. (modulus - 1) + where + modulus = 1 `shiftL` (e - 2) + largeMod = 1 `shiftL` (2 * e - 1) + xPower = case someNatVal largeMod of + SomeNat (_ :: Proxy largeMod) -> + toInteger (unMod (fromInteger x ^ (2 * modulus) :: Mod largeMod)) + -- | For elements of the multiplicative group \((\mathbb{Z}/n\mathbb{Z})^*\), a Dirichlet -- character evaluates to a root of unity. eval :: DirichletCharacter n -> MultMod n -> RootOfUnity eval (Generated ds) m = foldMap (evalFactor m') ds - where m' = getVal $ multElement m + where + m' = toInteger $ unMod $ multElement m -- | Evaluate each factor of the Dirichlet character. evalFactor :: Integer -> DirichletFactor -> RootOfUnity @@ -218,8 +222,8 @@ -- | We define `succ` and `pred` with more efficient implementations than -- @`toEnum` . (+1) . `fromEnum`@. instance KnownNat n => Enum (DirichletCharacter n) where - toEnum = indexToChar . fromIntegral - fromEnum = fromIntegral . characterNumber + toEnum = indexToChar . intToNatural + fromEnum = integerToInt . characterNumber succ x = makeChar x (characterNumber x + 1) pred x = makeChar x (characterNumber x - 1) @@ -236,12 +240,12 @@ -- | We have a (non-canonical) enumeration of dirichlet characters. characterNumber :: DirichletCharacter n -> Integer characterNumber (Generated y) = foldl' go 0 y - where go x (OddPrime p k _ a) = x * m + numerator (fromRootOfUnity a * fromIntegral m) - where p' = fromIntegral (unPrime p) + where go x (OddPrime p k _ a) = x * m + numerator (fromRootOfUnity a * (m % 1)) + where p' = naturalToInteger (unPrime p) m = p'^(k-1)*(p'-1) go x (TwoPower k a b) = x' * 2 + numerator (fromRootOfUnity a * 2) where m = bit (k-2) :: Integer - x' = x `shiftL` (k-2) + numerator (fromRootOfUnity b * fromIntegral m) + x' = x `shiftL` (k-2) + numerator (fromRootOfUnity b * (m % 1)) go x Two = x -- | Give the dirichlet character from its number. @@ -268,7 +272,7 @@ -- | Use one character to make many more: better than indicesToChars since it avoids recalculating -- some primitive roots bulkMakeChars :: (Integral a, Functor f) => DirichletCharacter n -> f a -> f (DirichletCharacter n) -bulkMakeChars x = fmap (Generated . unroll t . (`mod` m) . fromIntegral) +bulkMakeChars x = fmap (Generated . unroll t . (`mod` m) . fromIntegral') where (Product m, t) = templateFromCharacter x -- We assign each natural a unique Template, which can be decorated (eg in `unroll`) to @@ -312,9 +316,9 @@ unroll :: [Template] -> Natural -> [DirichletFactor] unroll t m = snd (mapAccumL func m t) where func :: Natural -> Template -> (Natural, DirichletFactor) - func a (OddTemplate p k g n) = (a1, OddPrime p k g (toRootOfUnity $ (toInteger a2) % (toInteger n))) + func a (OddTemplate p k g n) = (a1, OddPrime p k g (toRootOfUnity $ toInteger a2 % toInteger n)) where (a1,a2) = quotRem a n - func a (TwoPTemplate k n) = (b1, TwoPower k (toRootOfUnity $ (toInteger a2) % 2) (toRootOfUnity $ (toInteger b2) % (toInteger n))) + func a (TwoPTemplate k n) = (b1, TwoPower k (toRootOfUnity $ toInteger a2 % 2) (toRootOfUnity $ toInteger b2 % toInteger n)) where (a1,a2) = quotRem a 2 (b1,b2) = quotRem a1 n func a TwoTemplate = (a, Two) @@ -328,11 +332,9 @@ -- reduced to \(a \bmod{d}\). Thus, the multiplicative function on \(\mathbb{Z}/d\mathbb{Z}\) -- induces a multiplicative function on \(\mathbb{Z}/n\mathbb{Z}\). -- --- >>> :set -XTypeApplications +-- >>> :set -XTypeApplications -XDataKinds -- >>> chi = indexToChar 5 :: DirichletCharacter 45 --- >>> chi2 = induced @135 chi --- >>> :t chi2 --- Maybe (DirichletCharacter 135) +-- >>> chi2 = induced @135 chi :: Maybe (DirichletCharacter 135) induced :: forall n d. (KnownNat d, KnownNat n) => DirichletCharacter d -> Maybe (DirichletCharacter n) induced (Generated start) = if n `rem` d == 0 then Just (Generated (combine (snd $ mkTemplate n) start)) @@ -355,7 +357,6 @@ newFactor TwoTemplate = Two newFactor (TwoPTemplate k _) = TwoPower k mempty mempty newFactor (OddTemplate p k g _) = OddPrime p k g mempty - -- rest (p,k) = OddPrime p k (generator p k) mempty -- | The gives a real Dirichlet -- character for odd moduli. @@ -365,7 +366,7 @@ else Nothing where n = natVal (Proxy :: Proxy n) go :: Template -> DirichletFactor - go (OddTemplate p k g _) = OddPrime p k g $ toRootOfUnity ((toInteger k) % 2) + go (OddTemplate p k g _) = OddPrime p k g $ toRootOfUnity (toInteger k % 2) -- jacobi symbol of a primitive root mod p over p is always -1 go _ = error "internal error in jacobiCharacter: please report this as a bug" -- every factor of n should be odd @@ -401,7 +402,7 @@ validChar :: forall n. KnownNat n => DirichletCharacter n -> Bool validChar (Generated xs) = correctDecomposition && all correctPrimitiveRoot xs && all validValued xs where correctDecomposition = sort (factorise n) == map getPP xs - getPP (TwoPower k _ _) = (two, fromIntegral k) + getPP (TwoPower k _ _) = (two, intToWord k) getPP (OddPrime p k _ _) = (p, k) getPP Two = (two,1) correctPrimitiveRoot (OddPrime p k g _) = g == generator p k @@ -502,7 +503,7 @@ -- of a dirichlet character are required, `evalAll` will be better than `evalGeneral`, since -- computations can be shared. evalAll :: forall n. KnownNat n => DirichletCharacter n -> Vector (OrZero RootOfUnity) -evalAll (Generated xs) = V.generate (fromIntegral n) func +evalAll (Generated xs) = V.generate (naturalToInt n) func where n = natVal (Proxy :: Proxy n) vectors = map mkVector xs func :: Int -> OrZero RootOfUnity @@ -511,10 +512,10 @@ go (modulus,v) = v ! (m `mod` modulus) mkVector :: DirichletFactor -> (Int, Vector (OrZero RootOfUnity)) mkVector Two = (2, V.fromList [Zero, mempty]) - mkVector (OddPrime p k (fromIntegral -> g) a) = (modulus, w) + mkVector (OddPrime p k (naturalToInt -> g) a) = (modulus, w) where p' = unPrime p - modulus = fromIntegral (p'^k) :: Int + modulus = naturalToInt (p'^k) :: Int w = V.create $ do v <- MV.replicate modulus Zero -- TODO: we're in the ST monad here anyway, could be better to use STRefs to manage @@ -546,11 +547,11 @@ -- | Attempt to construct a character from its table of values. -- An inverse to `evalAll`, defined only on its image. fromTable :: forall n. KnownNat n => Vector (OrZero RootOfUnity) -> Maybe (DirichletCharacter n) -fromTable v = if length v == fromIntegral n - then Generated <$> traverse makeFactor tmpl >>= check +fromTable v = if length v == naturalToInt n + then traverse makeFactor tmpl >>= check . Generated else Nothing where n = natVal (Proxy :: Proxy n) - n' = fromIntegral n :: Integer + n' = naturalToInteger n :: Integer tmpl = snd (mkTemplate n) check :: DirichletCharacter n -> Maybe (DirichletCharacter n) check chi = if evalAll chi == v then Just chi else Nothing @@ -558,8 +559,8 @@ makeFactor TwoTemplate = Just Two makeFactor (TwoPTemplate k _) = TwoPower k <$> getValue (-1,bit k) <*> getValue (exp4 k, bit k) makeFactor (OddTemplate p k g _) = OddPrime p k g <$> getValue (toInteger g, toInteger (unPrime p)^k) - getValue :: (Integer,Integer) -> Maybe RootOfUnity - getValue (g,m) = getAp (v ! fromInteger (fromJust (chinese (g,m) (1,n' `quot` m)) `mod` n')) + getValue :: (Integer, Integer) -> Maybe RootOfUnity + getValue (g, m) = getAp (v ! fromInteger (fst (fromJust (chinese (g, m) (1, n' `quot` m))) `mod` n')) exp4terms :: [Rational] exp4terms = [4^k % product [1..k] | k <- [0..]] @@ -569,4 +570,8 @@ -- In particular, lambda (exp4 n) n == 1 (for n >= 3) -- I've verified this for 3 <= n <= 2000, so the reasoning in fromTable should be accurate for moduli below 2^2000 exp4 :: Int -> Integer -exp4 n = (`mod` bit n) $ sum $ map (`mod` bit n) $ map (\q -> numerator q * fromMaybe (error "error in exp4") (recipMod (denominator q) (bit n))) $ take n $ exp4terms +exp4 n + = (`mod` bit n) + $ sum + $ map (\q -> (numerator q * fromMaybe (error "error in exp4") (recipMod (denominator q) (bit n))) `mod` bit n) + $ take n exp4terms diff --git a/Math/NumberTheory/Euclidean.hs b/Math/NumberTheory/Euclidean.hs deleted file mode 100644 --- a/Math/NumberTheory/Euclidean.hs +++ /dev/null @@ -1,64 +0,0 @@ --- | --- Module: Math.NumberTheory.Euclidean --- Copyright: (c) 2018 Alexandre Rodrigues Baldé --- Licence: MIT --- Maintainer: Alexandre Rodrigues Baldé --- Description: Deprecated --- --- This module exports a class to represent Euclidean domains. --- - -{-# LANGUAGE BangPatterns #-} -{-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE GeneralizedNewtypeDeriving #-} -{-# LANGUAGE MagicHash #-} -{-# LANGUAGE ScopedTypeVariables #-} - -module Math.NumberTheory.Euclidean {-# DEPRECATED "Use Data.Euclidean instead" #-} - ( GcdDomain(..) - , Euclidean(..) - , WrappedIntegral(..) - , extendedGCD - , isUnit - ) where - -import Prelude hiding (divMod, div, gcd, lcm, mod, quotRem, quot, rem) -import Data.Euclidean -import Data.Maybe -import Data.Semiring (Semiring(..), isZero) - --- | Check whether an element is a unit of the ring. -isUnit :: (Eq a, GcdDomain a) => a -> Bool -isUnit x = not (isZero x) && isJust (one `divide` x) - --- | Calculate the greatest common divisor of two numbers and coefficients --- for the linear combination. --- --- For signed types satisfies: --- --- > case extendedGCD a b of --- > (d, u, v) -> u*a + v*b == d --- > && d == gcd a b --- --- For unsigned and bounded types the property above holds, but since @u@ and @v@ must also be unsigned, --- the result may look weird. E. g., on 64-bit architecture --- --- > extendedGCD (2 :: Word) (3 :: Word) == (1, 2^64-1, 1) --- --- For unsigned and unbounded types (like 'Numeric.Natural.Natural') the result is undefined. --- --- For signed types we also have --- --- > abs u < abs b || abs b <= 1 --- > --- > abs v < abs a || abs a <= 1 --- --- (except if one of @a@ and @b@ is 'minBound' of a signed type). -extendedGCD :: (Eq a, Num a, Euclidean a) => a -> a -> (a, a, a) -extendedGCD a b = (d, x * signum a, y * signum b) - where - (d, x, y) = eGCD 0 1 1 0 (abs a) (abs b) - eGCD !n1 o1 !n2 o2 r s - | s == 0 = (r, o1, o2) - | otherwise = case r `quotRem` s of - (q, t) -> eGCD (o1 - q*n1) n1 (o2 - q*n2) n2 s t diff --git a/Math/NumberTheory/Euclidean/Coprimes.hs b/Math/NumberTheory/Euclidean/Coprimes.hs --- a/Math/NumberTheory/Euclidean/Coprimes.hs +++ b/Math/NumberTheory/Euclidean/Coprimes.hs @@ -52,7 +52,7 @@ g = gcd x y (x', g', xgs) = doPair (x `unsafeDivide` g) xm g (xm + ym) - xgs' = if isUnit g' then xgs else ((g', xm + ym) : xgs) + xgs' = if isUnit g' then xgs else (g', xm + ym) : xgs (y', rests) = mapAccumL go (y `unsafeDivide` g) xgs' go w (t, tm) = (w', if isUnit t' || tm == 0 then acc else (t', tm) : acc) @@ -64,11 +64,11 @@ = isJust (x `divide` x') && isJust (y `divide` y') && coprime x' y' - && all (coprime x') (map fst rest) - && all (coprime y') (map fst rest) - && all (not . isUnit) (map fst rest) + && all (coprime x' . fst) rest + && all (coprime y' . fst) rest + && not (any (isUnit . fst) rest) && and [ coprime s t | (s, _) : ts <- tails rest, (t, _) <- ts ] - && abs ((x ^ xm) * (y ^ ym)) == abs ((x' ^ xm) * (y' ^ ym) * product (map (\(r, k) -> r ^ k) rest)) + && abs ((x ^ xm) * (y ^ ym)) == abs ((x' ^ xm) * (y' ^ ym) * product (map (uncurry (^)) rest)) where (x', y', rest) = doPair x xm y ym diff --git a/Math/NumberTheory/Moduli/Cbrt.hs b/Math/NumberTheory/Moduli/Cbrt.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Moduli/Cbrt.hs @@ -0,0 +1,160 @@ +-- | +-- Module: Math.NumberTheory.Moduli.Cbrt +-- Copyright: (c) 2020 Federico Bongiorno +-- Licence: MIT +-- Maintainer: Federico Bongiorno +-- +-- +-- of two Eisenstein Integers. + +{-# LANGUAGE LambdaCase #-} + +module Math.NumberTheory.Moduli.Cbrt + ( CubicSymbol(..) + , cubicSymbol + , symbolToNum + ) where + +import Math.NumberTheory.Quadratic.EisensteinIntegers +import Math.NumberTheory.Utils.FromIntegral +import qualified Data.Euclidean as A +import Math.NumberTheory.Utils +import Data.Semigroup + +-- | Represents the +-- +-- It is either @0@, @ω@, @ω²@ or @1@. +data CubicSymbol = Zero | Omega | OmegaSquare | One deriving (Eq) + +-- | The set of cubic symbols form a semigroup. Note `stimes` +-- is allowed to take non-positive values. In other words, the set +-- of non-zero cubic symbols is regarded as a group. +-- +-- >>> import Data.Semigroup +-- >>> stimes (-1) Omega +-- ω² +-- >>> stimes 0 Zero +-- 1 +instance Semigroup CubicSymbol where + Zero <> _ = Zero + _ <> Zero = Zero + One <> y = y + x <> One = x + Omega <> Omega = OmegaSquare + Omega <> OmegaSquare = One + OmegaSquare <> Omega = One + OmegaSquare <> OmegaSquare = Omega + stimes k n = case (k `mod` 3, n) of + (0, _) -> One + (1, symbol) -> symbol + (2, Omega) -> OmegaSquare + (2, OmegaSquare) -> Omega + (2, symbol) -> symbol + _ -> error "Math.NumberTheory.Moduli.Cbrt: exponentiation undefined." + +instance Show CubicSymbol where + show = \case + Zero -> "0" + Omega -> "ω" + OmegaSquare -> "ω²" + One -> "1" + +-- | Converts a +-- +-- to an Eisenstein Integer. +symbolToNum :: CubicSymbol -> EisensteinInteger +symbolToNum = \case + Zero -> 0 + Omega -> ω + OmegaSquare -> -1 - ω + One -> 1 + +-- The algorithm `cubicSymbol` is adapted from +-- . +-- It is divided in the following steps. +-- +-- (1) Check whether @beta@ is coprime to 3. +-- (2) Replace @alpha@ by the remainder of @alpha@ mod @beta@ +-- This does not affect the cubic symbol. +-- (3) Replace @alpha@ and @beta@ by their associated primary +-- divisors and keep track of how their cubic residue changes. +-- (4) Check if any of the two numbers is a zero or a unit. In this +-- case, return their cubic residue. +-- (5) Otherwise, invoke cubic reciprocity by swapping @alpha@ and +-- @beta@. Note both numbers have to be primary. +-- Return to Step 2. + +-- | +-- of two Eisenstein Integers. +-- The first argument is the numerator and the second argument +-- is the denominator. The latter must be coprime to @3@. +-- This condition is checked. +-- +-- If the arguments have a common factor, the result +-- is 'Zero', otherwise it is either 'Omega', 'OmegaSquare' or 'One'. +-- +-- >>> cubicSymbol (45 + 23*ω) (11 - 30*ω) +-- 0 +-- >>> cubicSymbol (31 - ω) (1 +10*ω) +-- ω +cubicSymbol :: EisensteinInteger -> EisensteinInteger -> CubicSymbol +cubicSymbol alpha beta = case beta `A.rem` (1 - ω) of + -- This checks whether beta is coprime to 3, i.e. divisible by @1 - ω@ + -- In particular, it returns an error if @beta == 0@ + 0 -> error "Math.NumberTheory.Moduli.Cbrt: denominator is not coprime to 3." + _ -> cubicSymbolHelper alpha beta + +cubicSymbolHelper :: EisensteinInteger -> EisensteinInteger -> CubicSymbol +cubicSymbolHelper alpha beta = cubicReciprocity primaryRemainder primaryBeta <> newSymbol + where + (primaryRemainder, primaryBeta, newSymbol) = extractPrimaryContributions remainder beta + remainder = A.rem alpha beta + +cubicReciprocity :: EisensteinInteger -> EisensteinInteger -> CubicSymbol +-- Note @cubicReciprocity 0 1 = One@. It is better to adopt this convention. +cubicReciprocity _ 1 = One +-- Checks if first argument is zero. Note the second argument is never zero. +cubicReciprocity 0 _ = Zero +-- This checks if the first argument is a unit. Because it's primary, +-- it is enough to pattern match with 1. +cubicReciprocity 1 _ = One +-- Otherwise, cubic reciprocity is called. +cubicReciprocity alpha beta = cubicSymbolHelper beta alpha + +-- | This function takes two Eisenstein intgers @alpha@ and @beta@ and returns +-- three arguments @(gamma, delta, newSymbol)@. @gamma@ and @delta@ are the +-- associated primary numbers of alpha and beta respectively. @newSymbol@ +-- is the cubic symbol measuring the discrepancy between the cubic residue +-- of @alpha@ and @beta@, and the cubic residue of @gamma@ and @delta@. +extractPrimaryContributions :: EisensteinInteger -> EisensteinInteger -> (EisensteinInteger, EisensteinInteger, CubicSymbol) +extractPrimaryContributions alpha beta = (gamma, delta, newSymbol) + where + newSymbol = stimes (j * m) Omega <> stimes (- m - n) i + m :+ n = A.quot (delta - 1) 3 + (i, gamma) = getPrimaryDecomposition alphaThreeFree + (_, delta) = getPrimaryDecomposition beta + j = wordToInteger jIntWord + -- This function outputs data such that + -- @(1 - ω)^jIntWord * alphaThreeFree = alpha@. + (jIntWord, alphaThreeFree) = splitOff (1 - ω) alpha + +-- | This function takes an Eisenstein number @e@ and returns @(symbol, delta)@ +-- where @delta@ is its associated primary integer and @symbol@ is the +-- cubic symbol discrepancy between @e@ and @delta@. @delta@ is defined to be +-- the unique associated Eisenstein Integer to @e@ such that +-- \( \textrm{delta} \equiv 1 (\textrm{mod} 3) \). +-- Note that @delta@ exists if and only if @e@ is coprime to 3. In this +-- case, an error message is displayed. +getPrimaryDecomposition :: EisensteinInteger -> (CubicSymbol, EisensteinInteger) +-- This is the case where a common factor between @alpha@ and @beta@ is detected. +-- In this instance @cubicReciprocity@ will return `Zero`. +-- Strictly speaking, this is not a primary decomposition. +getPrimaryDecomposition 0 = (Zero, 0) +getPrimaryDecomposition e = case e `A.rem` 3 of + 1 -> (One, e) + 1 :+ 1 -> (OmegaSquare, -ω * e) + 0 :+ 1 -> (Omega, (-1 - ω) * e) + (-1) :+ 0 -> (One, -e) + (-1) :+ (-1) -> (OmegaSquare, ω * e) + 0 :+ (-1) -> (Omega, (1 + ω) * e) + _ -> error "Math.NumberTheory.Moduli.Cbrt: primary decomposition failed." diff --git a/Math/NumberTheory/Moduli/Chinese.hs b/Math/NumberTheory/Moduli/Chinese.hs --- a/Math/NumberTheory/Moduli/Chinese.hs +++ b/Math/NumberTheory/Moduli/Chinese.hs @@ -7,34 +7,17 @@ -- Chinese remainder theorem -- -{-# LANGUAGE BangPatterns #-} -{-# LANGUAGE CPP #-} {-# LANGUAGE RankNTypes #-} {-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE TupleSections #-} -{-# LANGUAGE TypeApplications #-} -{-# LANGUAGE TypeOperators #-} -#if __GLASGOW_HASKELL__ > 805 -{-# LANGUAGE NoStarIsType #-} -#endif - module Math.NumberTheory.Moduli.Chinese ( -- * Safe interface chinese - , chineseCoprime , chineseSomeMod - , chineseCoprimeSomeMod - - , -- * Unsafe interface - chineseRemainder - , chineseRemainder2 ) where import Prelude hiding ((^), (+), (-), (*), rem, mod, quot, gcd, lcm) -import qualified Prelude -import Control.Monad (foldM) import Data.Euclidean import Data.Mod import Data.Ratio @@ -42,94 +25,42 @@ import GHC.TypeNats (KnownNat, natVal) import Math.NumberTheory.Moduli.SomeMod -import Math.NumberTheory.Utils (recipMod) --- | 'chineseCoprime' @(n1, m1)@ @(n2, m2)@ returns @n@ such that --- @n \`mod\` m1 == n1@ and @n \`mod\` m2 == n2@. --- Moduli @m1@ and @m2@ must be coprime, otherwise 'Nothing' is returned. --- --- This function is slightly faster than 'chinese', but more restricted. --- --- >>> chineseCoprime (1, 2) (2, 3) --- Just 5 --- >>> chineseCoprime (3, 4) (5, 6) --- Nothing -- moduli must be coprime -chineseCoprime :: (Eq a, Ring a, Euclidean a) => (a, a) -> (a, a) -> Maybe a -chineseCoprime (n1, m1) (n2, m2) - | d == one - = Just $ (v * m2 * n1 + u * m1 * n2) `rem` (m1 * m2) - | otherwise = Nothing - where - (d, u, v) = extendedGCD m1 m2 -{-# DEPRECATED chineseCoprime "Use 'chinese' instead" #-} - --- | 'chinese' @(n1, m1)@ @(n2, m2)@ returns @n@ such that +-- | 'chinese' @(n1, m1)@ @(n2, m2)@ returns @(n, lcm m1 m2)@ such that -- @n \`mod\` m1 == n1@ and @n \`mod\` m2 == n2@, if exists. -- Moduli @m1@ and @m2@ are allowed to have common factors. -- -- >>> chinese (1, 2) (2, 3) --- Just 5 +-- Just (-1, 6) -- >>> chinese (3, 4) (5, 6) --- Just 11 +-- Just (-1, 12) -- >>> chinese (3, 4) (2, 6) -- Nothing -chinese :: forall a. (Eq a, Ring a, Euclidean a) => (a, a) -> (a, a) -> Maybe a +chinese :: forall a. (Eq a, Ring a, Euclidean a) => (a, a) -> (a, a) -> Maybe (a, a) chinese (n1, m1) (n2, m2) | d == one - = Just $ (v * m2 * n1 + u * m1 * n2) `rem` (m1 * m2) + = Just ((v * m2 * n1 + u * m1 * n2) `rem` m, m) | (n1 - n2) `rem` d == zero - = Just $ (v * (m2 `quot` d) * n1 + u * (m1 `quot` d) * n2) `rem` ((m1 `quot` d) * m2) + = Just ((v * (m2 `quot` d) * n1 + u * (m1 `quot` d) * n2) `rem` m, m) | otherwise = Nothing where (d, u, v) = extendedGCD m1 m2 + m = if d == one then m1 * m2 else (m1 `quot` d) * m2 -{-# SPECIALISE chinese :: (Int, Int) -> (Int, Int) -> Maybe Int #-} -{-# SPECIALISE chinese :: (Word, Word) -> (Word, Word) -> Maybe Word #-} -{-# SPECIALISE chinese :: (Integer, Integer) -> (Integer, Integer) -> Maybe Integer #-} +{-# SPECIALISE chinese :: (Int, Int) -> (Int, Int) -> Maybe (Int, Int) #-} +{-# SPECIALISE chinese :: (Word, Word) -> (Word, Word) -> Maybe (Word, Word) #-} +{-# SPECIALISE chinese :: (Integer, Integer) -> (Integer, Integer) -> Maybe (Integer, Integer) #-} isCompatible :: KnownNat m => Mod m -> Rational -> Bool isCompatible n r = case invertMod (fromInteger (denominator r)) of Nothing -> False Just r' -> r' * fromInteger (numerator r) == n -chineseWrap - :: (Integer -> Integer -> Integer) - -> ((Integer, Integer) -> (Integer, Integer) -> Maybe Integer) - -> SomeMod - -> SomeMod - -> Maybe SomeMod -chineseWrap f g (SomeMod n1) (SomeMod n2) - = fmap (`modulo` fromInteger (f m1 m2)) (g (toInteger $ unMod n1, m1) (toInteger $ unMod n2, m2)) - where - m1 = toInteger $ natVal n1 - m2 = toInteger $ natVal n2 -chineseWrap _ _ (SomeMod n) (InfMod r) - | isCompatible n r = Just $ InfMod r - | otherwise = Nothing -chineseWrap _ _ (InfMod r) (SomeMod n) - | isCompatible n r = Just $ InfMod r - | otherwise = Nothing -chineseWrap _ _ (InfMod r1) (InfMod r2) - | r1 == r2 = Just $ InfMod r1 - | otherwise = Nothing - --- | Same as 'chineseCoprime', but operates on residues. --- --- >>> :set -XDataKinds --- >>> import Math.NumberTheory.Moduli.Class --- >>> (1 `modulo` 2) `chineseCoprimeSomeMod` (2 `modulo` 3) --- Just (5 `modulo` 6) --- >>> (3 `modulo` 4) `chineseCoprimeSomeMod` (5 `modulo` 6) --- Nothing -chineseCoprimeSomeMod :: SomeMod -> SomeMod -> Maybe SomeMod -chineseCoprimeSomeMod = chineseWrap (*) chineseCoprime -{-# DEPRECATED chineseCoprimeSomeMod "Use 'chineseSomeMod' instead" #-} - -- | Same as 'chinese', but operates on residues. -- -- >>> :set -XDataKinds --- >>> import Math.NumberTheory.Moduli.Class +-- >>> import Data.Mod -- >>> (1 `modulo` 2) `chineseSomeMod` (2 `modulo` 3) -- Just (5 `modulo` 6) -- >>> (3 `modulo` 4) `chineseSomeMod` (5 `modulo` 6) @@ -137,44 +68,19 @@ -- >>> (3 `modulo` 4) `chineseSomeMod` (2 `modulo` 6) -- Nothing chineseSomeMod :: SomeMod -> SomeMod -> Maybe SomeMod -chineseSomeMod = chineseWrap lcm chinese - -------------------------------------------------------------------------------- --- Unsafe interface - --- | Given a list @[(r_1,m_1), ..., (r_n,m_n)]@ of @(residue,modulus)@ --- pairs, @chineseRemainder@ calculates the solution to the simultaneous --- congruences --- --- > --- > r ≡ r_k (mod m_k) --- > --- --- if all moduli are positive and pairwise coprime. Otherwise --- the result is @Nothing@ regardless of whether --- a solution exists. -chineseRemainder :: [(Integer, Integer)] -> Maybe Integer -chineseRemainder remainders = foldM addRem 0 remainders - where - !modulus = product (map snd remainders) - addRem acc (_,1) = Just acc - addRem acc (r,m) = do - let cf = modulus `quot` m - inv <- recipMod cf m - Just $! (acc + inv*cf*r) `rem` modulus -{-# DEPRECATED chineseRemainder "Use 'chinese' instead" #-} - --- | @chineseRemainder2 (r_1,m_1) (r_2,m_2)@ calculates the solution of --- --- > --- > r ≡ r_k (mod m_k) --- --- if @m_1@ and @m_2@ are coprime. -chineseRemainder2 :: (Integer, Integer) -> (Integer, Integer) -> Integer -chineseRemainder2 (n1, m1) (n2, m2) = ((1 - u * m1) * n1 + (1 - v * m2) * n2) `Prelude.mod` (m1 * m2) - where - (_, u, v) = extendedGCD m1 m2 -{-# DEPRECATED chineseRemainder2 "Use 'chinese' instead" #-} +chineseSomeMod (SomeMod n1) (SomeMod n2) + = (\(n, m) -> n `modulo` fromInteger m) <$> chinese + (toInteger $ unMod n1, toInteger $ natVal n1) + (toInteger $ unMod n2, toInteger $ natVal n2) +chineseSomeMod (SomeMod n) (InfMod r) + | isCompatible n r = Just $ InfMod r + | otherwise = Nothing +chineseSomeMod (InfMod r) (SomeMod n) + | isCompatible n r = Just $ InfMod r + | otherwise = Nothing +chineseSomeMod (InfMod r1) (InfMod r2) + | r1 == r2 = Just $ InfMod r1 + | otherwise = Nothing ------------------------------------------------------------------------------- -- Utils diff --git a/Math/NumberTheory/Moduli/Class.hs b/Math/NumberTheory/Moduli/Class.hs --- a/Math/NumberTheory/Moduli/Class.hs +++ b/Math/NumberTheory/Moduli/Class.hs @@ -7,17 +7,11 @@ -- Safe modular arithmetic with modulo on type level. -- -{-# LANGUAGE BangPatterns #-} {-# LANGUAGE DataKinds #-} {-# LANGUAGE GADTs #-} -{-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE KindSignatures #-} -{-# LANGUAGE LambdaCase #-} -{-# LANGUAGE MagicHash #-} {-# LANGUAGE RankNTypes #-} {-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE StandaloneDeriving #-} -{-# LANGUAGE UnboxedTuples #-} module Math.NumberTheory.Moduli.Class ( -- * Known modulo diff --git a/Math/NumberTheory/Moduli/DiscreteLogarithm.hs b/Math/NumberTheory/Moduli/DiscreteLogarithm.hs deleted file mode 100644 --- a/Math/NumberTheory/Moduli/DiscreteLogarithm.hs +++ /dev/null @@ -1,14 +0,0 @@ --- | --- Module: Math.NumberTheory.Moduli.DiscreteLogarithm --- Copyright: (c) 2018 Bhavik Mehta --- License: MIT --- Maintainer: Andrew Lelechenko --- Description: Deprecated --- - -module Math.NumberTheory.Moduli.DiscreteLogarithm - {-# DEPRECATED "Use Math.NumberTheory.Moduli.Multiplicative instead" #-} - ( discreteLogarithm - ) where - -import Math.NumberTheory.Moduli.Multiplicative diff --git a/Math/NumberTheory/Moduli/Equations.hs b/Math/NumberTheory/Moduli/Equations.hs --- a/Math/NumberTheory/Moduli/Equations.hs +++ b/Math/NumberTheory/Moduli/Equations.hs @@ -82,7 +82,7 @@ combine :: [([Integer], Integer)] -> ([Integer], Integer) combine = foldl - (\(xs, xm) (ys, ym) -> ([ fromJust $ chinese (x, xm) (y, ym) | x <- xs, y <- ys ], xm * ym)) + (\(xs, xm) (ys, ym) -> ([ fst $ fromJust $ chinese (x, xm) (y, ym) | x <- xs, y <- ys ], xm * ym)) ([0], 1) solveQuadraticPrimePower @@ -130,10 +130,8 @@ | a `rem` p' == 0 = solveLinear' p' b c | otherwise - = map (\n -> n * recipModInteger (2 * a) p' `mod` p') - $ map (subtract b) + = map (\n -> (n - b) * recipModInteger (2 * a) p' `mod` p') $ sqrtsModPrime (b * b - 4 * a * c) p where p' :: Integer p' = unPrime p - diff --git a/Math/NumberTheory/Moduli/Internal.hs b/Math/NumberTheory/Moduli/Internal.hs --- a/Math/NumberTheory/Moduli/Internal.hs +++ b/Math/NumberTheory/Moduli/Internal.hs @@ -23,13 +23,12 @@ import GHC.Integer.GMP.Internals import Numeric.Natural -import Math.NumberTheory.ArithmeticFunctions import Math.NumberTheory.Moduli.Chinese import Math.NumberTheory.Moduli.Equations import Math.NumberTheory.Moduli.Singleton import Math.NumberTheory.Primes -import Math.NumberTheory.Powers.Modular import Math.NumberTheory.Roots +import Math.NumberTheory.Utils.FromIntegral -- https://en.wikipedia.org/wiki/Primitive_root_modulo_n#Finding_primitive_roots isPrimitiveRoot' @@ -44,14 +43,17 @@ CGOddPrimePower p k -> oddPrimePowerTest (unPrime p) k r CGDoubleOddPrimePower p k -> doubleOddPrimePowerTest (unPrime p) k r where - oddPrimeTest p g = let phi = totient p - pows = map (\pk -> phi `quot` unPrime (fst pk)) (factorise phi) - exps = map (\x -> powMod g x p) pows - in g /= 0 && gcd g p == 1 && notElem 1 exps oddPrimePowerTest p 1 g = oddPrimeTest p (g `mod` p) - oddPrimePowerTest p _ g = oddPrimeTest p (g `mod` p) && powMod g (p-1) (p*p) /= 1 + oddPrimePowerTest p _ g = oddPrimeTest p (g `mod` p) && case someNatVal (fromIntegral' (p * p)) of + SomeNat (_ :: Proxy pp) -> fromIntegral g ^ (p - 1) /= (1 :: Mod pp) + doubleOddPrimePowerTest p k g = odd g && oddPrimePowerTest p k g + oddPrimeTest p g = g /= 0 && gcd g p == 1 && case someNatVal (fromIntegral' p) of + SomeNat (_ :: Proxy p) -> all (\x -> fromIntegral g ^ x /= (1 :: Mod p)) pows + where + pows = map (\(q, _) -> (p - 1) `quot` unPrime q) (factorise (p - 1)) + -- Implementation of Bach reduction (https://www2.eecs.berkeley.edu/Pubs/TechRpts/1984/CSD-84-186.pdf) {-# INLINE discreteLogarithmPP #-} discreteLogarithmPP :: Integer -> Word -> Integer -> Integer -> Natural @@ -62,9 +64,8 @@ thetaA = theta p pkMinusOne a thetaB = theta p pkMinusOne b pkMinusOne = p^(k-1) - pkMinusPk1 = pkMinusOne * (p - 1) c = (recipModInteger thetaA pkMinusOne * thetaB) `rem` pkMinusOne - result = fromJust $ chinese (baseSol, p-1) (c, pkMinusOne) + (result, pkMinusPk1) = fromJust $ chinese (baseSol, p-1) (c, pkMinusOne) -- compute the homomorphism theta given in https://math.stackexchange.com/a/1864495/418148 {-# INLINE theta #-} @@ -83,7 +84,7 @@ -- made redundant, since n would be prime. discreteLogarithmPrime :: Integer -> Integer -> Integer -> Natural discreteLogarithmPrime p a b - | p < 100000000 = fromIntegral $ discreteLogarithmPrimeBSGS (fromInteger p) (fromInteger a) (fromInteger b) + | p < 100000000 = intToNatural $ discreteLogarithmPrimeBSGS (fromInteger p) (fromInteger a) (fromInteger b) | otherwise = discreteLogarithmPrimePollard p a b discreteLogarithmPrimeBSGS :: Int -> Int -> Int -> Int diff --git a/Math/NumberTheory/Moduli/Jacobi.hs b/Math/NumberTheory/Moduli/Jacobi.hs deleted file mode 100644 --- a/Math/NumberTheory/Moduli/Jacobi.hs +++ /dev/null @@ -1,18 +0,0 @@ --- | --- Module: Math.NumberTheory.Moduli.Jacobi --- Copyright: (c) 2011 Daniel Fischer, 2017-2018 Andrew Lelechenko --- Licence: MIT --- Maintainer: Andrew Lelechenko --- Description: Deprecated --- --- --- is a generalization of the Legendre symbol, useful for primality --- testing and integer factorization. --- - -module Math.NumberTheory.Moduli.Jacobi - {-# DEPRECATED "Use Math.NumberTheory.Moduli.Sqrt instead" #-} - ( module Math.NumberTheory.Moduli.JacobiSymbol - ) where - -import Math.NumberTheory.Moduli.JacobiSymbol 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 @@ -64,8 +64,8 @@ -- If arguments have a common factor, the result -- is 'Zero', otherwise it is 'MinusOne' or 'One'. -- --- >>> jacobi 1001 9911 --- Zero -- arguments have a common factor 11 +-- >>> jacobi 1001 9911 -- arguments have a common factor 11 +-- Zero -- >>> jacobi 1001 9907 -- MinusOne {-# SPECIALISE jacobi :: Integer -> Integer -> JacobiSymbol, diff --git a/Math/NumberTheory/Moduli/Multiplicative.hs b/Math/NumberTheory/Moduli/Multiplicative.hs --- a/Math/NumberTheory/Moduli/Multiplicative.hs +++ b/Math/NumberTheory/Moduli/Multiplicative.hs @@ -9,7 +9,6 @@ {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE ViewPatterns #-} -{-# LANGUAGE PatternSynonyms #-} module Math.NumberTheory.Moduli.Multiplicative ( -- * Multiplicative group @@ -36,7 +35,7 @@ import Math.NumberTheory.Primes -- | This type represents elements of the multiplicative group mod m, i.e. --- those elements which are coprime to m. Use @toMultElement@ to construct. +-- those elements which are coprime to m. Use @isMultElement@ to construct. newtype MultMod m = MultMod { multElement :: Mod m -- ^ Unwrap a residue. } deriving (Eq, Ord, Show) diff --git a/Math/NumberTheory/Moduli/PrimitiveRoot.hs b/Math/NumberTheory/Moduli/PrimitiveRoot.hs deleted file mode 100644 --- a/Math/NumberTheory/Moduli/PrimitiveRoot.hs +++ /dev/null @@ -1,19 +0,0 @@ --- | --- Module: Math.NumberTheory.Moduli.PrimitiveRoot --- Copyright: (c) 2017 Andrew Lelechenko, 2018 Bhavik Mehta --- Licence: MIT --- Maintainer: Andrew Lelechenko --- Description: Deprecated --- --- Primitive roots and cyclic groups. --- - -module Math.NumberTheory.Moduli.PrimitiveRoot - {-# DEPRECATED "Use Math.NumberTheory.Moduli.Multiplicative instead" #-} - ( -- * Primitive roots - PrimitiveRoot - , unPrimitiveRoot - , isPrimitiveRoot - ) where - -import Math.NumberTheory.Moduli.Multiplicative 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 @@ -12,13 +12,11 @@ {-# LANGUAGE DeriveGeneric #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE GADTs #-} -{-# LANGUAGE KindSignatures #-} {-# LANGUAGE LambdaCase #-} {-# LANGUAGE PatternSynonyms #-} {-# LANGUAGE PolyKinds #-} {-# LANGUAGE RankNTypes #-} {-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE StandaloneDeriving #-} {-# LANGUAGE TupleSections #-} {-# LANGUAGE TypeOperators #-} {-# LANGUAGE ViewPatterns #-} @@ -64,6 +62,7 @@ import Math.NumberTheory.Roots (highestPower) import Math.NumberTheory.Primes import Math.NumberTheory.Primes.Types +import Math.NumberTheory.Utils.FromIntegral -- | Wrapper to hide an unknown type-level natural. data Some (a :: Nat -> Type) where @@ -107,7 +106,7 @@ -- -- >>> :set -XDataKinds -- >>> sfactors :: SFactors Integer 13 --- SFactors {sfactorsFactors = [(Prime 13,1)]} +-- SFactors {unSFactors = [(Prime 13,1)]} sfactors :: forall a m. (Ord a, UniqueFactorisation a, KnownNat m) => SFactors a m sfactors = if m == 0 then error "sfactors: modulo must be positive" @@ -120,7 +119,7 @@ -- -- >>> import Math.NumberTheory.Primes -- >>> someSFactors (factorise 98) --- SFactors {sfactorsFactors = [(Prime 2,1),(Prime 7,2)]} +-- SFactors {unSFactors = [(Prime 2,1),(Prime 7,2)]} someSFactors :: (Ord a, Num a) => [(Prime a, Word)] -> Some (SFactors a) someSFactors = Some @@ -134,7 +133,7 @@ -- > toModulo :: SFactors Integer m -> Natural -- > toModulo t = case proofFromSFactors t of Sub Dict -> natVal t proofFromSFactors :: Integral a => SFactors a m -> (() :- KnownNat m) -proofFromSFactors (SFactors fs) = Sub $ unsafeCoerce (Magic Dict) (fromIntegral (factorBack fs) :: Natural) +proofFromSFactors (SFactors fs) = Sub $ unsafeCoerce (Magic Dict) (fromIntegral' (factorBack fs) :: Natural) -- | This singleton data type establishes a correspondence -- between a modulo @m@ on type level @@ -195,14 +194,16 @@ -- -- >>> :set -XDataKinds -- >>> import Data.Maybe --- >>> cyclicGroup :: CyclicGroup Integer 169 --- CGOddPrimePower' (Prime 13) 2 +-- >>> cyclicGroup :: Maybe (CyclicGroup Integer 169) +-- Just (CGOddPrimePower' (Prime 13) 2) -- --- >>> sfactorsToCyclicGroup (fromModulo 4) +-- >>> :set -XTypeOperators -XNoStarIsType +-- >>> import GHC.TypeNats +-- >>> sfactorsToCyclicGroup (sfactors :: SFactors Integer 4) -- Just CG4' --- >>> sfactorsToCyclicGroup (fromModulo (2 * 13 ^ 3)) +-- >>> sfactorsToCyclicGroup (sfactors :: SFactors Integer (2 * 13 ^ 3)) -- Just (CGDoubleOddPrimePower' (Prime 13) 3) --- >>> sfactorsToCyclicGroup (fromModulo (4 * 13)) +-- >>> sfactorsToCyclicGroup (sfactors :: SFactors Integer (4 * 13)) -- Nothing cyclicGroup :: forall a m. @@ -268,11 +269,13 @@ -- | Check whether a multiplicative group of residues, -- characterized by its modulo, is cyclic and, if yes, return its form. -- --- >>> sfactorsToCyclicGroup (fromModulo 4) +-- >>> :set -XTypeOperators -XNoStarIsType +-- >>> import GHC.TypeNats +-- >>> sfactorsToCyclicGroup (sfactors :: SFactors Integer 4) -- Just CG4' --- >>> sfactorsToCyclicGroup (fromModulo (2 * 13 ^ 3)) +-- >>> sfactorsToCyclicGroup (sfactors :: SFactors Integer (2 * 13 ^ 3)) -- Just (CGDoubleOddPrimePower' (Prime 13) 3) --- >>> sfactorsToCyclicGroup (fromModulo (4 * 13)) +-- >>> sfactorsToCyclicGroup (sfactors :: SFactors Integer (4 * 13)) -- Nothing sfactorsToCyclicGroup :: (Eq a, Num a) => SFactors a m -> Maybe (CyclicGroup a m) sfactorsToCyclicGroup (SFactors fs) = case fs of @@ -287,8 +290,8 @@ -- | Invert 'sfactorsToCyclicGroup'. -- -- >>> import Data.Maybe --- >>> cyclicGroupToSFactors (fromJust (sfactorsToCyclicGroup (fromModulo 4))) --- SFactors {sfactorsModulo = 4, sfactorsFactors = [(Prime 2,2)]} +-- >>> cyclicGroupToSFactors (fromJust (sfactorsToCyclicGroup (sfactors :: SFactors Integer 4))) +-- SFactors {unSFactors = [(Prime 2,2)]} cyclicGroupToSFactors :: Num a => CyclicGroup a m -> SFactors a m cyclicGroupToSFactors = SFactors . \case CG2' -> [(Prime 2, 1)] 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 @@ -39,6 +39,7 @@ -- (1 `modulo` 5) -- >>> (2 `modulo` 10) * (4 `modulo` 15) -- (3 `modulo` 5) +-- >>> import Data.Ratio -- >>> 2 `modulo` 10 + fromRational (3 % 7) -- (1 `modulo` 10) -- >>> 2 `modulo` 10 * fromRational (3 % 7) @@ -148,7 +149,7 @@ fromRational = InfMod {-# INLINE fromRational #-} recip x = case invertSomeMod x of - Nothing -> error $ "recip{SomeMod}: residue is not coprime with modulo" + Nothing -> error "recip{SomeMod}: residue is not coprime with modulo" Just y -> y -- | See the warning about division above. @@ -170,10 +171,11 @@ -- | Computes the inverse value, if it exists. -- --- >>> invertSomeMod (3 `modulo` 10) --- Just (7 `modulo` 10) -- because 3 * 7 = 1 :: Mod 10 +-- >>> invertSomeMod (3 `modulo` 10) -- because 3 * 7 = 1 :: Mod 10 +-- Just (7 `modulo` 10) -- >>> invertSomeMod (4 `modulo` 10) -- Nothing +-- >>> import Data.Ratio -- >>> invertSomeMod (fromRational (2 % 5)) -- Just 5 % 2 invertSomeMod :: SomeMod -> Maybe SomeMod diff --git a/Math/NumberTheory/Moduli/Sqrt.hs b/Math/NumberTheory/Moduli/Sqrt.hs --- a/Math/NumberTheory/Moduli/Sqrt.hs +++ b/Math/NumberTheory/Moduli/Sqrt.hs @@ -8,9 +8,11 @@ -- . -- -{-# LANGUAGE BangPatterns #-} -{-# LANGUAGE ViewPatterns #-} -{-# LANGUAGE CPP #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE TupleSections #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeApplications #-} +{-# LANGUAGE ViewPatterns #-} module Math.NumberTheory.Moduli.Sqrt ( -- * Modular square roots @@ -29,13 +31,14 @@ import Data.Constraint import Data.Maybe import Data.Mod +import Data.Proxy +import GHC.TypeNats (KnownNat, SomeNat(..), natVal, someNatVal) import Math.NumberTheory.Moduli.Chinese import Math.NumberTheory.Moduli.JacobiSymbol import Math.NumberTheory.Moduli.Singleton -import Math.NumberTheory.Powers.Modular (powMod) import Math.NumberTheory.Primes -import Math.NumberTheory.Utils (shiftToOddCount, splitOff, recipMod) +import Math.NumberTheory.Utils (shiftToOddCount, splitOff) import Math.NumberTheory.Utils.FromIntegral -- | List all modular square roots. @@ -60,15 +63,14 @@ ms = map (\(p, pow) -> unPrime p ^ pow) pps rs :: [[Integer]] - rs = map (\(p, pow) -> sqrtsModPrimePower n p pow) pps + rs = map (uncurry (sqrtsModPrimePower n)) pps cs :: [[(Integer, Integer)]] - cs = zipWith (\l m -> map (\x -> (x, m)) l) rs ms + cs = zipWith (\l m -> map (, m) l) rs ms - comb t1@(_, m1) t2@(_, m2) = (if ch < 0 then ch + m else ch, m) + comb t1 t2 = (if ch < 0 then ch + m else ch, m) where - ch = fromJust $ chinese t1 t2 - m = m1 * m2 + (ch, m) = fromJust $ chinese t1 t2 -- | List all square roots modulo the power of a prime. -- @@ -118,7 +120,8 @@ sqrtsModPrime n (unPrime -> prime) = case jacobi n prime of MinusOne -> [] Zero -> [0] - One -> let r = sqrtModP' (n `mod` prime) prime in [r, prime - r] + One -> case someNatVal (fromInteger prime) of + SomeNat (_ :: Proxy p) -> let r = toInteger (unMod (sqrtModP' @p (fromInteger n))) in [r, prime - r] ------------------------------------------------------------------------------- -- Internals @@ -126,75 +129,82 @@ -- | @sqrtModP' square prime@ finds a square root of @square@ modulo -- prime. @prime@ /must/ be a (positive) prime, and @square@ /must/ be a positive -- quadratic residue modulo @prime@, i.e. @'jacobi square prime == 1@. -sqrtModP' :: Integer -> Integer -> Integer -sqrtModP' square prime - | prime == 2 = square - | rem4 prime == 3 = powMod square ((prime + 1) `quot` 4) prime - | square `mod` prime == prime - 1 - = sqrtOfMinusOne prime - | otherwise = tonelliShanks square prime +sqrtModP' :: KnownNat p => Mod p -> Mod p +sqrtModP' square + | prime == 2 = square + | rem4 prime == 3 = square ^ ((prime + 1) `quot` 4) + | square == maxBound = sqrtOfMinusOne + | otherwise = tonelliShanks square + where + prime = natVal square -- | @p@ must be of form @4k + 1@ -sqrtOfMinusOne :: Integer -> Integer -sqrtOfMinusOne p - = head - $ filter (\n -> n /= 1 && n /= p - 1) - $ map (\n -> powMod n k p) - [2..p-2] +sqrtOfMinusOne :: KnownNat p => Mod p +sqrtOfMinusOne = res where + p = natVal res k = (p - 1) `quot` 4 + res = head + $ dropWhile (\n -> n == 1 || n == maxBound) + $ map (^ k) [2 .. maxBound - 1] -- | @tonelliShanks square prime@ calculates a square root of @square@ -- modulo @prime@, where @prime@ is a prime of the form @4*k + 1@ and -- @square@ is a positive quadratic residue modulo @prime@, using the -- Tonelli-Shanks algorithm. -tonelliShanks :: Integer -> Integer -> Integer -tonelliShanks square prime = loop rc t1 generator log2 +tonelliShanks :: forall p. KnownNat p => Mod p -> Mod p +tonelliShanks square = loop rc t1 generator log2 where - (wordToInt -> log2,q) = shiftToOddCount (prime-1) - nonSquare = findNonSquare prime - generator = powMod nonSquare q prime - rc = powMod square ((q+1) `quot` 2) prime - t1 = powMod square q prime - msqr x = (x*x) `rem` prime + prime = natVal square + (log2, q) = shiftToOddCount (prime - 1) + generator = findNonSquare ^ q + rc = square ^ ((q + 1) `quot` 2) + t1 = square ^ q + msquare 0 x = x - msquare k x = msquare (k-1) (msqr x) + msquare k x = msquare (k-1) (x * x) + findPeriod per 1 = per - findPeriod per x = findPeriod (per+1) (msqr x) + findPeriod per x = findPeriod (per + 1) (x * x) - loop :: Integer -> Integer -> Integer -> Int -> Integer + loop :: Mod p -> Mod p -> Mod p -> Word -> Mod p loop !r t c m | t == 1 = r | otherwise = loop nextR nextT nextC nextM where nextM = findPeriod 0 t b = msquare (m - 1 - nextM) c - nextR = (r*b) `rem` prime - nextC = msqr b - nextT = (t*nextC) `rem` prime + nextR = r * b + nextC = b * b + nextT = t * nextC -- | prime must be odd, n must be coprime with prime sqrtModPP' :: Integer -> Integer -> Word -> Maybe Integer sqrtModPP' n prime expo = case jacobi n prime of MinusOne -> Nothing Zero -> Nothing - One -> fixup $ sqrtModP' (n `mod` prime) prime + One -> case someNatVal (fromInteger prime) of + SomeNat (_ :: Proxy p) -> Just $ fixup $ sqrtModP' @p (fromInteger n) where - fixup r = let diff' = r*r-n - in if diff' == 0 - then Just r - else case splitOff prime diff' of - (e,q) | expo <= e -> Just r - | otherwise -> fmap (\inv -> hoist inv r (q `mod` prime) (prime^e)) (recipMod (2*r) prime) + fixup :: KnownNat p => Mod p -> Integer + fixup r + | diff' == 0 = r' + | expo <= e = r' + | otherwise = hoist (recip (2 * r)) r' (fromInteger q) (prime^e) + where + r' = toInteger (unMod r) + diff' = r' * r' - n + (e, q) = splitOff prime diff' + hoist :: KnownNat p => Mod p -> Integer -> Mod p -> Integer -> Integer hoist inv root elim pp - | diff' == 0 = root' - | expo <= ex = root' - | otherwise = hoist inv root' (nelim `mod` prime) (prime^ex) - where - root' = (root + (inv*(prime-elim))*pp) `mod` (prime*pp) - diff' = root'*root' - n - (ex, nelim) = splitOff prime diff' + | diff' == 0 = root' + | expo <= ex = root' + | otherwise = hoist inv root' (fromInteger nelim) (prime ^ ex) + where + root' = root + toInteger (unMod (inv * negate elim)) * pp + diff' = root' * root' - n + (ex, nelim) = splitOff prime diff' -- dirty, dirty sqM2P :: Integer -> Word -> Maybe Integer @@ -202,7 +212,7 @@ | e < 2 = Just (n `mod` 2) | n' == 0 = Just 0 | odd k = Nothing - | otherwise = fmap ((`mod` mdl) . (`shiftL` wordToInt k2)) $ solve s e2 + | otherwise = (`mod` mdl) . (`shiftL` wordToInt k2) <$> solve s e2 where mdl = 1 `shiftL` wordToInt e n' = n `mod` mdl @@ -233,17 +243,17 @@ rem8 :: Integral a => a -> Int rem8 n = fromIntegral n .&. 7 -findNonSquare :: Integer -> Integer -findNonSquare n - | rem8 n == 5 || rem8 n == 3 = 2 - | otherwise = search candidates - where - -- It is enough to consider only prime candidates, but - -- the probability that the smallest non-residue is > 67 - -- is small and 'jacobi' test is fast, - -- so we use [71..n] instead of filter isPrime [71..n]. - candidates = 3:5:7:11:13:17:19:23:29:31:37:41:43:47:53:59:61:67:[71..n] - search (p:ps) = case jacobi p n of - MinusOne -> p - _ -> search ps - search _ = error "Should never have happened, prime list exhausted." +findNonSquare :: KnownNat n => Mod n +findNonSquare = res + where + n = natVal res + res + | rem8 n == 3 || rem8 n == 5 = 2 + | otherwise = fromIntegral $ head $ + dropWhile (\p -> jacobi p n /= MinusOne) candidates + + -- It is enough to consider only prime candidates, but + -- the probability that the smallest non-residue is > 67 + -- is small and 'jacobi' test is fast, + -- so we use [71..n] instead of filter isPrime [71..n]. + candidates = 3:5:7:11:13:17:19:23:29:31:37:41:43:47:53:59:61:67:[71..n] diff --git a/Math/NumberTheory/MoebiusInversion.hs b/Math/NumberTheory/MoebiusInversion.hs --- a/Math/NumberTheory/MoebiusInversion.hs +++ b/Math/NumberTheory/MoebiusInversion.hs @@ -115,7 +115,7 @@ small <- MG.unsafeNew (mk0 + 1) :: ST s (G.Mutable v s t) MG.unsafeWrite small 0 0 - MG.unsafeWrite small 1 $! (fun 1) + MG.unsafeWrite small 1 $! fun 1 when (mk0 >= 2) $ MG.unsafeWrite small 2 $! (fun 2 - fun 1) diff --git a/Math/NumberTheory/Powers.hs b/Math/NumberTheory/Powers.hs deleted file mode 100644 --- a/Math/NumberTheory/Powers.hs +++ /dev/null @@ -1,47 +0,0 @@ --- | --- Module: Math.NumberTheory.Powers --- Copyright: (c) 2011 Daniel Fischer --- Licence: MIT --- Maintainer: Daniel Fischer --- Description: Deprecated --- --- Calculating integer roots, modular powers and related things. --- This module reexports the most needed functions from the implementation --- modules. The implementation modules provide some additional functions, --- in particular some unsafe functions which omit some tests for performance --- reasons. --- - -{-# OPTIONS_GHC -fno-warn-deprecations #-} - -module Math.NumberTheory.Powers - {-# DEPRECATED "Use Math.NumberTheory.Roots or Math.NumberTheory.Powers.Modular" #-} - ( -- * Integer Roots - -- ** Square roots - integerSquareRoot - , isSquare - , exactSquareRoot - -- ** Cube roots - , integerCubeRoot - , isCube - , exactCubeRoot - -- ** Fourth roots - , integerFourthRoot - , isFourthPower - , exactFourthRoot - -- ** General roots - , integerRoot - , isKthPower - , exactRoot - , isPerfectPower - , highestPower - -- * Modular powers - , powMod - ) where - -import Math.NumberTheory.Powers.Squares -import Math.NumberTheory.Powers.Cubes -import Math.NumberTheory.Powers.Fourth -import Math.NumberTheory.Powers.General - -import Math.NumberTheory.Powers.Modular diff --git a/Math/NumberTheory/Powers/Cubes.hs b/Math/NumberTheory/Powers/Cubes.hs deleted file mode 100644 --- a/Math/NumberTheory/Powers/Cubes.hs +++ /dev/null @@ -1,205 +0,0 @@ --- | --- Module: Math.NumberTheory.Powers.Cubes --- Copyright: (c) 2011 Daniel Fischer --- Licence: MIT --- Maintainer: Daniel Fischer --- Description: Deprecated --- --- Functions dealing with cubes. Moderately efficient calculation of integer --- cube roots and testing for cubeness. - -{-# LANGUAGE MagicHash, BangPatterns, CPP, FlexibleContexts #-} - -module Math.NumberTheory.Powers.Cubes - {-# DEPRECATED "Use Math.NumberTheory.Roots" #-} - ( integerCubeRoot - , integerCubeRoot' - , exactCubeRoot - , isCube - , isCube' - , isPossibleCube - ) where - -#include "MachDeps.h" - -import Control.Monad.ST -import Data.Bits -import qualified Data.Vector.Unboxed as V -import qualified Data.Vector.Unboxed.Mutable as MV - -import GHC.Base -import GHC.Integer -import GHC.Integer.GMP.Internals -import GHC.Integer.Logarithms (integerLog2#) - -import Numeric.Natural - -import Math.NumberTheory.Roots - --- | Calculate the integer cube root of a nonnegative integer @n@, --- that is, the largest integer @r@ such that @r^3 <= n@. --- The precondition @n >= 0@ is not checked. -{-# RULES -"integerCubeRoot'/Int" integerCubeRoot' = cubeRootInt' -"integerCubeRoot'/Word" integerCubeRoot' = cubeRootWord -"integerCubeRoot'/Integer" integerCubeRoot' = cubeRootIgr - #-} -{-# INLINE [1] integerCubeRoot' #-} -integerCubeRoot' :: Integral a => a -> a -integerCubeRoot' 0 = 0 -integerCubeRoot' n = newton3 n (approxCuRt n) - --- | Test whether a nonnegative integer is a cube. --- Before 'integerCubeRoot' is calculated, a few tests --- of remainders modulo small primes weed out most non-cubes. --- For testing many numbers, most of which aren't cubes, --- this is much faster than @let r = cubeRoot n in r*r*r == n@. --- The condition @n >= 0@ is /not/ checked. -{-# SPECIALISE isCube' :: Int -> Bool, - Word -> Bool, - Integer -> Bool, - Natural -> Bool - #-} -isCube' :: Integral a => a -> Bool -isCube' !n = isPossibleCube n - && (r*r*r == n) - where - r = integerCubeRoot' n - --- | Test whether a nonnegative number is possibly a cube. --- Only about 0.08% of all numbers pass this test. --- The precondition @n >= 0@ is /not/ checked. -{-# SPECIALISE isPossibleCube :: Int -> Bool, - Word -> Bool, - Integer -> Bool, - Natural -> Bool - #-} -isPossibleCube :: Integral a => a -> Bool -isPossibleCube !n - = V.unsafeIndex cr512 (fromIntegral n .&. 511) - && V.unsafeIndex cubeRes837 (fromIntegral (n `rem` 837)) - && V.unsafeIndex cubeRes637 (fromIntegral (n `rem` 637)) - && V.unsafeIndex cubeRes703 (fromIntegral (n `rem` 703)) - ----------------------------------------------------------------------- --- Utility Functions -- ----------------------------------------------------------------------- - --- Special case for 'Int', a little faster. --- For @n <= 2^64@, the truncated 'Double' is never --- more than one off. Things might overflow for @n@ --- close to @maxBound@, so check for overflow. -cubeRootInt' :: Int -> Int -cubeRootInt' 0 = 0 -cubeRootInt' n - | n < c || c < 0 = r-1 - | 0 < d && d < n = r+1 - | otherwise = r - where - x = fromIntegral n :: Double - r = truncate (x ** (1/3)) - c = r*r*r - d = c+3*r*(r+1) - -cubeRootWord :: Word -> Word -cubeRootWord 0 = 0 -cubeRootWord w -#if WORD_SIZE_IN_BITS == 64 - | r > 2642245 = 2642245 -#else - | r > 1625 = 1625 -#endif - | w < c = r-1 - | c < w && e < w && c < e = r+1 - | otherwise = r - where - r = truncate ((fromIntegral w) ** (1/3) :: Double) - c = r*r*r - d = 3*r*(r+1) - e = c+d - -cubeRootIgr :: Integer -> Integer -cubeRootIgr 0 = 0 -cubeRootIgr n = newton3 n (approxCuRt n) - -{-# SPECIALISE newton3 :: Integer -> Integer -> Integer #-} -newton3 :: Integral a => a -> a -> a -newton3 n a = go (step a) - where - step k = (2*k + n `quot` (k*k)) `quot` 3 - go k - | m < k = go m - | otherwise = k - where - m = step k - -{-# SPECIALISE approxCuRt :: Integer -> Integer #-} -approxCuRt :: Integral a => a -> a -approxCuRt 0 = 0 -approxCuRt n = fromInteger $ appCuRt (fromIntegral n) - --- threshold for shifting vs. direct fromInteger --- we shift when we expect more than 256 bits -#if WORD_SIZE_IN_BITS == 64 -#define THRESH 5 -#else -#define THRESH 9 -#endif - --- | approximate cube root, about 50 bits should be correct for large numbers -appCuRt :: Integer -> Integer -appCuRt (S# i#) = case double2Int# (int2Double# i# **## (1.0## /## 3.0##)) of - r# -> S# r# -appCuRt n@(Jp# bn#) - | isTrue# ((sizeofBigNat# bn#) <# THRESH#) = - floor (fromInteger n ** (1.0/3.0) :: Double) - | otherwise = case integerLog2# n of - l# -> case (l# `quotInt#` 3#) -# 51# of - h# -> case shiftRInteger n (3# *# h#) of - m -> case floor (fromInteger m ** (1.0/3.0) :: Double) of - r -> shiftLInteger r h# --- There's already handling for negative in integerCubeRoot, --- but integerCubeRoot' is exported directly too. -appCuRt _ = error "integerCubeRoot': negative argument" - --- not very discriminating, but cheap, so it's an overall gain -cr512 :: V.Vector Bool -cr512 = runST $ do - ar <- MV.replicate 512 True - let note s i - | i < 512 = MV.unsafeWrite ar i False >> note s (i+s) - | otherwise = return () - note 4 2 - note 8 4 - note 32 16 - note 64 32 - note 256 128 - MV.unsafeWrite ar 256 False - V.unsafeFreeze ar - --- Remainders modulo @3^3 * 31@ -cubeRes837 :: V.Vector Bool -cubeRes837 = runST $ do - ar <- MV.replicate 837 False - let note 837 = return () - note k = MV.unsafeWrite ar ((k*k*k) `rem` 837) True >> note (k+1) - note 0 - V.unsafeFreeze ar - --- Remainders modulo @7^2 * 13@ -cubeRes637 :: V.Vector Bool -cubeRes637 = runST $ do - ar <- MV.replicate 637 False - let note 637 = return () - note k = MV.unsafeWrite ar ((k*k*k) `rem` 637) True >> note (k+1) - note 0 - V.unsafeFreeze ar - --- Remainders modulo @19 * 37@ -cubeRes703 :: V.Vector Bool -cubeRes703 = runST $ do - ar <- MV.replicate 703 False - let note 703 = return () - note k = MV.unsafeWrite ar ((k*k*k) `rem` 703) True >> note (k+1) - note 0 - V.unsafeFreeze ar diff --git a/Math/NumberTheory/Powers/Fourth.hs b/Math/NumberTheory/Powers/Fourth.hs deleted file mode 100644 --- a/Math/NumberTheory/Powers/Fourth.hs +++ /dev/null @@ -1,216 +0,0 @@ --- | --- Module: Math.NumberTheory.Powers.Squares --- Copyright: (c) 2011 Daniel Fischer --- Licence: MIT --- Maintainer: Daniel Fischer --- Description: Deprecated --- --- Functions dealing with fourth powers. Efficient calculation of integer fourth --- roots and efficient testing for being a square's square. - -{-# LANGUAGE MagicHash, CPP, FlexibleContexts #-} - -module Math.NumberTheory.Powers.Fourth - {-# DEPRECATED "Use Math.NumberTheory.Roots" #-} - ( integerFourthRoot - , integerFourthRoot' - , exactFourthRoot - , isFourthPower - , isFourthPower' - , isPossibleFourthPower - ) where - -#include "MachDeps.h" - -import Control.Monad.ST -import Data.Bits -import qualified Data.Vector.Unboxed as V -import qualified Data.Vector.Unboxed.Mutable as MV - -import GHC.Base -import GHC.Integer -import GHC.Integer.GMP.Internals -import GHC.Integer.Logarithms (integerLog2#) - -import Numeric.Natural - -import Math.NumberTheory.Roots - --- | Calculate the integer fourth root of a nonnegative number, --- that is, the largest integer @r@ with @r^4 <= n@. --- Throws an error on negaitve input. -{-# SPECIALISE integerFourthRoot :: Int -> Int, - Word -> Word, - Integer -> Integer, - Natural -> Natural - #-} -integerFourthRoot :: Integral a => a -> a -integerFourthRoot = integerRoot (4 :: Word) - --- | Calculate the integer fourth root of a nonnegative number, --- that is, the largest integer @r@ with @r^4 <= n@. --- The condition is /not/ checked. -{-# RULES -"integerFourthRoot'/Int" integerFourthRoot' = biSqrtInt -"integerFourthRoot'/Word" integerFourthRoot' = biSqrtWord -"integerFourthRoot'/Integer" integerFourthRoot' = biSqrtIgr - #-} -{-# INLINE [1] integerFourthRoot' #-} -integerFourthRoot' :: Integral a => a -> a -integerFourthRoot' 0 = 0 -integerFourthRoot' n = newton4 n (approxBiSqrt n) - --- | Returns @Nothing@ if @n@ is not a fourth power, --- @Just r@ if @n == r^4@ and @r >= 0@. -{-# SPECIALISE exactFourthRoot :: Int -> Maybe Int, - Word -> Maybe Word, - Integer -> Maybe Integer, - Natural -> Maybe Natural - #-} -exactFourthRoot :: Integral a => a -> Maybe a -exactFourthRoot = exactRoot (4 :: Word) - --- | Test whether an integer is a fourth power. --- First nonnegativity is checked, then the unchecked --- test is called. -{-# SPECIALISE isFourthPower :: Int -> Bool, - Word -> Bool, - Integer -> Bool, - Natural -> Bool - #-} -isFourthPower :: Integral a => a -> Bool -isFourthPower = isKthPower (4 :: Word) - --- | Test whether a nonnegative number is a fourth power. --- The condition is /not/ checked. If a number passes the --- 'isPossibleFourthPower' test, its integer fourth root --- is calculated. -{-# SPECIALISE isFourthPower' :: Int -> Bool, - Word -> Bool, - Integer -> Bool, - Natural -> Bool - #-} -isFourthPower' :: Integral a => a -> Bool -isFourthPower' n = isPossibleFourthPower n && r2*r2 == n - where - r = integerFourthRoot' n - r2 = r*r - --- | Test whether a nonnegative number is a possible fourth power. --- The condition is /not/ checked. --- This eliminates about 99.958% of numbers. -{-# SPECIALISE isPossibleFourthPower :: Int -> Bool, - Word -> Bool, - Integer -> Bool, - Natural -> Bool - #-} -isPossibleFourthPower :: Integral a => a -> Bool -isPossibleFourthPower n - = V.unsafeIndex biSqRes256 (fromIntegral n .&. 255) - && V.unsafeIndex biSqRes425 (fromIntegral (n `rem` 425)) - && V.unsafeIndex biSqRes377 (fromIntegral (n `rem` 377)) - -{-# SPECIALISE newton4 :: Integer -> Integer -> Integer #-} -newton4 :: Integral a => a -> a -> a -newton4 n a = go (step a) - where - step k = (3*k + n `quot` (k*k*k)) `quot` 4 - go k - | m < k = go m - | otherwise = k - where - m = step k - -{-# SPECIALISE approxBiSqrt :: Integer -> Integer #-} -approxBiSqrt :: Integral a => a -> a -approxBiSqrt = fromInteger . appBiSqrt . fromIntegral - --- threshold for shifting vs. direct fromInteger --- we shift when we expect more than 384 bits -#if WORD_SIZE_IN_BITS == 64 -#define THRESH 7 -#else -#define THRESH 13 -#endif - --- Find a fairly good approximation to the fourth root. --- About 48 bits should be correct for large Integers. -appBiSqrt :: Integer -> Integer -appBiSqrt (S# i#) = S# (double2Int# (sqrtDouble# (sqrtDouble# (int2Double# i#)))) -appBiSqrt n@(Jp# bn#) - | isTrue# ((sizeofBigNat# bn#) <# THRESH#) = - floor (sqrt . sqrt $ fromInteger n :: Double) - | otherwise = case integerLog2# n of - l# -> case uncheckedIShiftRA# l# 2# -# 47# of - h# -> case shiftRInteger n (4# *# h#) of - m -> case floor (sqrt $ sqrt $ fromInteger m :: Double) of - r -> shiftLInteger r h# --- There's already a check for negative in integerFourthRoot, --- but integerFourthRoot' is exported directly too. -appBiSqrt _ = error "integerFourthRoot': negative argument" - - -biSqRes256 :: V.Vector Bool -biSqRes256 = runST $ do - ar <- MV.replicate 256 False - let note 257 = return () - note i = MV.unsafeWrite ar i True >> note (i+16) - MV.unsafeWrite ar 0 True - MV.unsafeWrite ar 16 True - note 1 - V.unsafeFreeze ar - -biSqRes425 :: V.Vector Bool -biSqRes425 = runST $ do - ar <- MV.replicate 425 False - let note 154 = return () - note i = MV.unsafeWrite ar ((i*i*i*i) `rem` 425) True >> note (i+1) - note 0 - V.unsafeFreeze ar - -biSqRes377 :: V.Vector Bool -biSqRes377 = runST $ do - ar <- MV.replicate 377 False - let note 144 = return () - note i = MV.unsafeWrite ar ((i*i*i*i) `rem` 377) True >> note (i+1) - note 0 - V.unsafeFreeze ar - -biSqrtInt :: Int -> Int -biSqrtInt 0 = 0 -biSqrtInt n -#if WORD_SIZE_IN_BITS == 64 - | r > 55108 = 55108 -#else - | r > 215 = 215 -#endif - | n < r4 = r-1 - | otherwise = r - where - x :: Double - x = fromIntegral n - -- timed faster than x**0.25, never too small - r = truncate (sqrt (sqrt x)) - r2 = r*r - r4 = r2*r2 - -biSqrtWord :: Word -> Word -biSqrtWord 0 = 0 -biSqrtWord n -#if WORD_SIZE_IN_BITS == 64 - | r > 65535 = 65535 -#else - | r > 255 = 255 -#endif - | n < r4 = r-1 - | otherwise = r - where - x :: Double - x = fromIntegral n - r = truncate (sqrt (sqrt x)) - r2 = r*r - r4 = r2*r2 - -biSqrtIgr :: Integer -> Integer -biSqrtIgr 0 = 0 -biSqrtIgr n = newton4 n (approxBiSqrt n) diff --git a/Math/NumberTheory/Powers/General.hs b/Math/NumberTheory/Powers/General.hs deleted file mode 100644 --- a/Math/NumberTheory/Powers/General.hs +++ /dev/null @@ -1,89 +0,0 @@ --- | --- Module: Math.NumberTheory.Powers.General --- Copyright: (c) 2011 Daniel Fischer --- Licence: MIT --- Maintainer: Daniel Fischer --- Description: Deprecated --- --- Calculating integer roots and determining perfect powers. --- The algorithms are moderately efficient. --- - -{-# LANGUAGE BangPatterns #-} -{-# LANGUAGE CPP #-} -{-# LANGUAGE MagicHash #-} -{-# LANGUAGE UnboxedTuples #-} - -{-# OPTIONS_GHC -fno-warn-deprecations #-} - -module Math.NumberTheory.Powers.General - {-# DEPRECATED "Use Math.NumberTheory.Roots" #-} - ( integerRoot - , exactRoot - , isKthPower - , isPerfectPower - , highestPower - , largePFPower - ) where - -#include "MachDeps.h" - -import Math.NumberTheory.Logarithms (integerLogBase') -import qualified Math.NumberTheory.Powers.Squares as P2 -import qualified Math.NumberTheory.Powers.Cubes as P3 -import qualified Math.NumberTheory.Powers.Fourth as P4 -import Math.NumberTheory.Utils.FromIntegral (intToWord) - -import Math.NumberTheory.Roots - --- | @'largePFPower' bd n@ produces the pair @(b,k)@ with the largest --- exponent @k@ such that @n == b^k@, where @bd > 1@ (it is expected --- that @bd@ is much larger, at least @1000@ or so), @n > bd^2@ and @n@ --- has no prime factors @p <= bd@, skipping the trial division phase --- of @'highestPower'@ when that is a priori known to be superfluous. --- It is only present to avoid duplication of work in factorisation --- and primality testing, it is not expected to be generally useful. --- The assumptions are not checked, if they are not satisfied, wrong --- results and wasted work may be the consequence. -largePFPower :: Integer -> Integer -> (Integer, Word) -largePFPower bd n = rawPower ln n - where - ln = intToWord (integerLogBase' (bd+1) n) - ------------------------------------------------------------------------------------------- --- Auxiliary functions -- ------------------------------------------------------------------------------------------- - -rawPower :: Word -> Integer -> (Integer, Word) -rawPower mx n - | mx < 2 = (n,1) - | mx == 2 = case P2.exactSquareRoot n of - Just r -> (r,2) - Nothing -> (n,1) -rawPower mx n = case P4.exactFourthRoot n of - Just r -> case rawPower (mx `quot` 4) r of - (m,e) -> (m, 4*e) - Nothing -> case P2.exactSquareRoot n of - Just r -> case rawOddPower (mx `quot` 2) r of - (m,e) -> (m, 2*e) - Nothing -> rawOddPower mx n - -rawOddPower :: Word -> Integer -> (Integer, Word) -rawOddPower mx n - | mx < 3 = (n,1) -rawOddPower mx n = case P3.exactCubeRoot n of - Just r -> case rawOddPower (mx `quot` 3) r of - (m,e) -> (m, 3*e) - Nothing -> badPower mx n - -badPower :: Word -> Integer -> (Integer, Word) -badPower mx n - | mx < 5 = (n,1) - | otherwise = go 1 mx n (takeWhile (<= mx) $ scanl (+) 5 $ cycle [2,4]) - where - go !e b m (k:ks) - | b < k = (m,e) - | otherwise = case exactRoot k m of - Just r -> go (e*k) (b `quot` k) r (k:ks) - Nothing -> go e b m ks - go e _ m [] = (m,e) diff --git a/Math/NumberTheory/Powers/Modular.hs b/Math/NumberTheory/Powers/Modular.hs --- a/Math/NumberTheory/Powers/Modular.hs +++ b/Math/NumberTheory/Powers/Modular.hs @@ -7,20 +7,15 @@ -- Modular powers (a. k. a. modular exponentiation). -- -{-# LANGUAGE CPP #-} -{-# LANGUAGE MagicHash #-} - module Math.NumberTheory.Powers.Modular + {-# DEPRECATED "Use Data.Mod or Data.Mod.Word instead" #-} ( powMod , powModWord , powModInt ) where -import qualified GHC.Integer.GMP.Internals as GMP (powModInteger) - -import GHC.Exts (Word(..)) import GHC.Natural (powModNatural) -import qualified GHC.Integer.GMP.Internals as GMP (powModWord) +import qualified GHC.Integer.GMP.Internals as GMP (powModInteger) import Math.NumberTheory.Utils.FromIntegral -- | @powMod@ @b@ @e@ @m@ computes (@b^e@) \`mod\` @m@ in effective way. @@ -40,12 +35,12 @@ -- need both to fit into machine word and to handle large moduli, -- take a look at 'powModInt' and 'powModWord'. -- --- >>> powMod 3 101 (2^60-1 :: Integer) --- 1018105167100379328 -- correct --- >>> powMod 3 101 (2^60-1 :: Int) --- 1115647832265427613 -- incorrect due to overflow --- >>> powModInt 3 101 (2^60-1 :: Int) --- 1018105167100379328 -- correct +-- >>> powMod 3 101 (2^60-1 :: Integer) -- correct +-- 1018105167100379328 +-- >>> powMod 3 101 (2^60-1 :: Int) -- incorrect due to overflow +-- 1115647832265427613 +-- >>> powModInt 3 101 (2^60-1 :: Int) -- correct +-- 1018105167100379328 powMod :: (Integral a, Integral b) => a -> b -> a -> a powMod x y m | m <= 0 = error "powModInt: non-positive modulo" @@ -54,7 +49,7 @@ where f _ 0 acc = acc f b e acc = f (b * b `rem` m) (e `quot` 2) - (if odd e then (b * acc `rem` m) else acc) + (if odd e then b * acc `rem` m else acc) {-# INLINE [1] powMod #-} {-# RULES @@ -76,7 +71,7 @@ -- >>> powModWord 3 101 (2^60-1) -- 1018105167100379328 powModWord :: Word -> Word -> Word -> Word -powModWord (W# x) (W# y) (W# m) = W# (GMP.powModWord x y m) +powModWord b e m = fromInteger $ GMP.powModInteger (toInteger b) (toInteger e) (toInteger m) -- | Specialised version of 'powMod', able to handle large moduli correctly. -- diff --git a/Math/NumberTheory/Powers/Squares.hs b/Math/NumberTheory/Powers/Squares.hs deleted file mode 100644 --- a/Math/NumberTheory/Powers/Squares.hs +++ /dev/null @@ -1,227 +0,0 @@ --- | --- Module: Math.NumberTheory.Powers.Squares --- Copyright: (c) 2011 Daniel Fischer --- Licence: MIT --- Maintainer: Daniel Fischer --- Description: Deprecated --- --- Functions dealing with squares. Efficient calculation of integer square roots --- and efficient testing for squareness. - -{-# LANGUAGE MagicHash, BangPatterns, PatternGuards, CPP, FlexibleContexts #-} - -{-# OPTIONS_GHC -fno-warn-deprecations #-} - -module Math.NumberTheory.Powers.Squares - {-# DEPRECATED "Use Math.NumberTheory.Roots" #-} - ( -- * Square root calculation - integerSquareRoot - , integerSquareRoot' - , integerSquareRootRem - , integerSquareRootRem' - , exactSquareRoot - -- * Tests for squares - , isSquare - , isSquare' - , isPossibleSquare - , isPossibleSquare2 - ) where - -#include "MachDeps.h" - -import Control.Monad.ST -import Data.Bits -import qualified Data.Vector.Unboxed as V -import qualified Data.Vector.Unboxed.Mutable as MV - -import Numeric.Natural - -import Math.NumberTheory.Powers.Squares.Internal - -import Math.NumberTheory.Roots - --- | Calculate the integer square root of a nonnegative number @n@, --- that is, the largest integer @r@ with @r*r <= n@. --- The precondition @n >= 0@ is not checked. -{-# RULES -"integerSquareRoot'/Int" integerSquareRoot' = isqrtInt' -"integerSquareRoot'/Word" integerSquareRoot' = isqrtWord -"integerSquareRoot'/Integer" integerSquareRoot' = isqrtInteger - #-} -{-# INLINE [1] integerSquareRoot' #-} -integerSquareRoot' :: Integral a => a -> a -integerSquareRoot' = isqrtA - --- | Calculate the integer square root of a nonnegative number as well as --- the difference of that number with the square of that root, that is if --- @(s,r) = integerSquareRootRem n@ then @s^2 <= n == s^2+r < (s+1)^2@. -{-# SPECIALISE integerSquareRootRem :: - Int -> (Int, Int), - Word -> (Word, Word), - Integer -> (Integer, Integer), - Natural -> (Natural, Natural) - #-} -integerSquareRootRem :: Integral a => a -> (a, a) -integerSquareRootRem n - | n < 0 = error "integerSquareRootRem: negative argument" - | otherwise = integerSquareRootRem' n - --- | Calculate the integer square root of a nonnegative number as well as --- the difference of that number with the square of that root, that is if --- @(s,r) = integerSquareRootRem' n@ then @s^2 <= n == s^2+r < (s+1)^2@. --- The precondition @n >= 0@ is not checked. -{-# RULES -"integerSquareRootRem'/Integer" integerSquareRootRem' = karatsubaSqrt - #-} -{-# INLINE [1] integerSquareRootRem' #-} -integerSquareRootRem' :: Integral a => a -> (a, a) -integerSquareRootRem' n = (s, n - s * s) - where - s = integerSquareRoot' n - --- | Test whether the input (a nonnegative number) @n@ is a square. --- The same as 'isSquare', but without the negativity test. --- Faster if many known positive numbers are tested. --- --- The precondition @n >= 0@ is not tested, passing negative --- arguments may cause any kind of havoc. -{-# SPECIALISE isSquare' :: Int -> Bool, - Word -> Bool, - Integer -> Bool, - Natural -> Bool - #-} -isSquare' :: Integral a => a -> Bool -isSquare' n - | isPossibleSquare n - , (_, 0) <- integerSquareRootRem' n = True - | otherwise = False - --- | Test whether a non-negative number may be a square. --- Non-negativity is not checked, passing negative arguments may --- cause any kind of havoc. --- --- First the remainder modulo 256 is checked (that can be calculated --- easily without division and eliminates about 82% of all numbers). --- After that, the remainders modulo 9, 25, 7, 11 and 13 are tested --- to eliminate altogether about 99.436% of all numbers. --- --- This is the test used by 'exactSquareRoot'. For large numbers, --- the slower but more discriminating test 'isPossibleSquare2' is --- faster. -{-# SPECIALISE isPossibleSquare :: Int -> Bool, - Word -> Bool, - Integer -> Bool, - Natural -> Bool - #-} -isPossibleSquare :: Integral a => a -> Bool -isPossibleSquare n - = V.unsafeIndex sr256 ((fromIntegral n) .&. 255) - && V.unsafeIndex sr693 (fromIntegral (n `rem` 693)) - && V.unsafeIndex sr325 (fromIntegral (n `rem` 325)) - --- | Test whether a non-negative number may be a square. --- Non-negativity is not checked, passing negative arguments may --- cause any kind of havoc. --- --- First the remainder modulo 256 is checked (that can be calculated --- easily without division and eliminates about 82% of all numbers). --- After that, the remainders modulo several small primes are tested --- to eliminate altogether about 99.98954% of all numbers. --- --- For smallish to medium sized numbers, this hardly performs better --- than 'isPossibleSquare', which uses smaller arrays, but for large --- numbers, where calculating the square root becomes more expensive, --- it is much faster (if the vast majority of tested numbers aren't squares). -{-# SPECIALISE isPossibleSquare2 :: Int -> Bool, - Word -> Bool, - Integer -> Bool, - Natural -> Bool - #-} -isPossibleSquare2 :: Integral a => a -> Bool -isPossibleSquare2 n - = V.unsafeIndex sr256 ((fromIntegral n) .&. 255) - && V.unsafeIndex sr819 (fromIntegral (n `rem` 819)) - && V.unsafeIndex sr1025 (fromIntegral (n `rem` 1025)) - && V.unsafeIndex sr2047 (fromIntegral (n `rem` 2047)) - && V.unsafeIndex sr4097 (fromIntegral (n `rem` 4097)) - && V.unsafeIndex sr341 (fromIntegral (n `rem` 341)) - ------------------------------------------------------------------------------ --- Auxiliary Stuff - --- Make an array indicating whether a remainder is a square remainder. -sqRemArray :: Int -> V.Vector Bool -sqRemArray md = runST $ do - ar <- MV.replicate md False - let !stop = (md `quot` 2) + 1 - fill k - | k < stop = MV.unsafeWrite ar ((k*k) `rem` md) True >> fill (k+1) - | otherwise = return () - MV.unsafeWrite ar 0 True - MV.unsafeWrite ar 1 True - fill 2 - V.unsafeFreeze ar - -sr256 :: V.Vector Bool -sr256 = sqRemArray 256 - -sr819 :: V.Vector Bool -sr819 = sqRemArray 819 - -sr4097 :: V.Vector Bool -sr4097 = sqRemArray 4097 - -sr341 :: V.Vector Bool -sr341 = sqRemArray 341 - -sr1025 :: V.Vector Bool -sr1025 = sqRemArray 1025 - -sr2047 :: V.Vector Bool -sr2047 = sqRemArray 2047 - -sr693 :: V.Vector Bool -sr693 = sqRemArray 693 - -sr325 :: V.Vector Bool -sr325 = sqRemArray 325 - --- Specialisations for Int, Word, and Integer - --- For @n <= 2^64@, the result of --- --- > truncate (sqrt $ fromIntegral n) --- --- is never too small and never more than one too large. --- The multiplication doesn't overflow for 32 or 64 bit Ints. -isqrtInt' :: Int -> Int -isqrtInt' n - | n < r*r = r-1 - | otherwise = r - where - !r = (truncate :: Double -> Int) . sqrt $ fromIntegral n --- With -O2, that should be translated to the below -{- -isqrtInt' n@(I# i#) - | r# *# r# ># i# = I# (r# -# 1#) - | otherwise = I# r# - where - !r# = double2Int# (sqrtDouble# (int2Double# i#)) --} - --- Same for Word. -isqrtWord :: Word -> Word -isqrtWord n - | n < (r*r) -#if WORD_SIZE_IN_BITS == 64 - || r == 4294967296 --- Double interprets values near maxBound as 2^64, we don't have that problem for 32 bits -#endif - = r-1 - | otherwise = r - where - !r = (fromIntegral :: Int -> Word) . (truncate :: Double -> Int) . sqrt $ fromIntegral n - -{-# INLINE isqrtInteger #-} -isqrtInteger :: Integer -> Integer -isqrtInteger = fst . karatsubaSqrt diff --git a/Math/NumberTheory/Powers/Squares/Internal.hs b/Math/NumberTheory/Powers/Squares/Internal.hs deleted file mode 100644 --- a/Math/NumberTheory/Powers/Squares/Internal.hs +++ /dev/null @@ -1,136 +0,0 @@ --- | --- Module: Math.NumberTheory.Powers.Squares.Internal --- Copyright: (c) 2016 Andrew Lelechenko --- Licence: MIT --- Maintainer: Andrew Lelechenko --- Description: Deprecated --- --- Internal functions dealing with square roots. End-users should not import this module. - -{-# LANGUAGE MagicHash #-} -{-# LANGUAGE BangPatterns #-} -{-# LANGUAGE PatternGuards #-} -{-# LANGUAGE CPP #-} -{-# LANGUAGE FlexibleContexts #-} - -{-# OPTIONS_HADDOCK hide #-} - -module Math.NumberTheory.Powers.Squares.Internal - {-# DEPRECATED "Use Math.NumberTheory.Roots" #-} - ( karatsubaSqrt - , isqrtA - ) where - -#include "MachDeps.h" - -import Data.Bits - -import GHC.Base -import GHC.Integer -import GHC.Integer.GMP.Internals -import GHC.Integer.Logarithms (integerLog2#) - -import Math.NumberTheory.Logarithms (integerLog2) - --- Find approximation to square root in 'Integer', then --- find the integer square root by the integer variant --- of Heron's method. Takes only a handful of steps --- unless the input is really large. -{-# SPECIALISE isqrtA :: Integer -> Integer #-} -isqrtA :: Integral a => a -> a -isqrtA 0 = 0 -isqrtA n = heron n (fromInteger . appSqrt . fromIntegral $ n) - --- Heron's method for integers. First make one step to ensure --- the value we're working on is @>= r@, then we have --- @k == r@ iff @k <= step k@. -{-# SPECIALISE heron :: Integer -> Integer -> Integer #-} -heron :: Integral a => a -> a -> a -heron n a = go (step a) - where - step k = (k + n `quot` k) `quot` 2 - go k - | m < k = go m - | otherwise = k - where - m = step k - --- threshold for shifting vs. direct fromInteger --- we shift when we expect more than 256 bits -#if WORD_SIZE_IN_BITS == 64 -#define THRESH 5 -#else -#define THRESH 9 -#endif - --- Find a fairly good approximation to the square root. --- At most one off for small Integers, about 48 bits should be correct --- for large Integers. -appSqrt :: Integer -> Integer -appSqrt (S# i#) = S# (double2Int# (sqrtDouble# (int2Double# i#))) -appSqrt n@(Jp# bn#) - | isTrue# ((sizeofBigNat# bn#) <# THRESH#) = - floor (sqrt $ fromInteger n :: Double) - | otherwise = case integerLog2# n of - l# -> case uncheckedIShiftRA# l# 1# -# 47# of - h# -> case shiftRInteger n (2# *# h#) of - m -> case floor (sqrt $ fromInteger m :: Double) of - r -> shiftLInteger r h# --- There's already a check for negative in integerSquareRoot, --- but integerSquareRoot' is exported directly too. -appSqrt _ = error "integerSquareRoot': negative argument" - - --- Integer square root with remainder, using the Karatsuba Square Root --- algorithm from --- Paul Zimmermann. Karatsuba Square Root. [Research Report] RR-3805, 1999, --- pp.8. - -karatsubaSqrt :: Integer -> (Integer, Integer) -karatsubaSqrt 0 = (0, 0) -karatsubaSqrt n - | lgN < 2300 = - let s = isqrtA n in (s, n - s * s) - | otherwise = - if lgN .&. 2 /= 0 then - karatsubaStep k (karatsubaSplit k n) - else - -- before we split n into 4 part we must ensure that the first part - -- is at least 2^k/4, since this doesn't happen here we scale n by - -- multiplying it by 4 - let n' = n `unsafeShiftL` 2 - (s, r) = karatsubaStep k (karatsubaSplit k n') - r' | s .&. 1 == 0 = r - | otherwise = r + double s - 1 - in (s `unsafeShiftR` 1, r' `unsafeShiftR` 2) - where - k = lgN `unsafeShiftR` 2 + 1 - lgN = integerLog2 n - -karatsubaStep :: Int -> (Integer, Integer, Integer, Integer) -> (Integer, Integer) -karatsubaStep k (a3, a2, a1, a0) - | r >= 0 = (s, r) - | otherwise = (s - 1, r + double s - 1) - where - r = cat u a0 - q * q - s = s' `unsafeShiftL` k + q - (q, u) = cat r' a1 `quotRem` double s' - (s', r') = karatsubaSqrt (cat a3 a2) - cat x y = x `unsafeShiftL` k .|. y - {-# INLINE cat #-} - -karatsubaSplit :: Int -> Integer -> (Integer, Integer, Integer, Integer) -karatsubaSplit k n0 = (a3, a2, a1, a0) - where - a3 = n3 - n3 = n2 `unsafeShiftR` k - a2 = n2 .&. m - n2 = n1 `unsafeShiftR` k - a1 = n1 .&. m - n1 = n0 `unsafeShiftR` k - a0 = n0 .&. m - m = 1 `unsafeShiftL` k - 1 - -double :: Bits a => a -> a -double x = x `unsafeShiftL` 1 -{-# INLINE double #-} diff --git a/Math/NumberTheory/Prefactored.hs b/Math/NumberTheory/Prefactored.hs --- a/Math/NumberTheory/Prefactored.hs +++ b/Math/NumberTheory/Prefactored.hs @@ -7,7 +7,6 @@ -- Type for numbers, accompanied by their factorisation. -- -{-# LANGUAGE CPP #-} {-# LANGUAGE TypeFamilies #-} {-# OPTIONS_GHC -fno-warn-unused-imports #-} @@ -66,7 +65,7 @@ -- >>> prefFactors $ totient (p^2 * q^3) -- Coprimes {unCoprimes = [(1000000000000000000000000000057,1),(41666666666666666666666666669,1),(2000000000000000000000000000071,2),(111111111111111111111111111115,1),(2,4),(3,3)]} -- >>> prefFactors $ totient $ totient (p^2 * q^3) --- Coprimes {unCoprimes = [(39521,1),(6046667,1),(22222222222222222222222222223,1),(2000000000000000000000000000071,1),(361696272343,1),(85331809838489,1),(227098769,1),(199937,1),(5,3),(41666666666666666666666666669,1),(2,22),(3,8)]} +-- Coprimes {unCoprimes = [(39521,1),(227098769,1),(22222222222222222222222222223,1),(2000000000000000000000000000071,1),(361696272343,1),(85331809838489,1),(6046667,1),(199937,1),(5,3),(41666666666666666666666666669,1),(2,22),(3,8)]} -- -- Pairwise coprimality of factors is crucial, because it allows -- us to process them independently, possibly even @@ -99,6 +98,13 @@ -- Prefactored {prefValue = 88045650000, prefFactors = Coprimes {unCoprimes = [(28,2),(33,3),(5,5)]}} fromFactors :: Semiring a => Coprimes a Word -> Prefactored a fromFactors as = Prefactored (getMul $ foldMap (\(a, k) -> Mul $ a ^ k) (unCoprimes as)) as + +instance (Eq a, GcdDomain a) => Semiring (Prefactored a) where + Prefactored v1 _ `plus` Prefactored v2 _ + = fromValue (v1 `plus` v2) + Prefactored v1 f1 `times` Prefactored v2 f2 + = Prefactored (v1 `times` v2) (f1 <> f2) + fromNatural n = fromValue (fromNatural n) instance (Eq a, Num a, GcdDomain a) => Num (Prefactored a) where Prefactored v1 _ + Prefactored v2 _ diff --git a/Math/NumberTheory/Primes.hs b/Math/NumberTheory/Primes.hs --- a/Math/NumberTheory/Primes.hs +++ b/Math/NumberTheory/Primes.hs @@ -5,7 +5,6 @@ -- Maintainer: Andrew Lelechenko -- -{-# LANGUAGE CPP #-} {-# LANGUAGE FlexibleInstances #-} {-# LANGUAGE LambdaCase #-} @@ -14,6 +13,7 @@ module Math.NumberTheory.Primes ( Prime , unPrime + , toPrimeIntegral , nextPrime , precPrime , UniqueFactorisation(..) @@ -123,7 +123,7 @@ -- > precPrime 97 == 97 precPrime :: (Bits a, Integral a, UniqueFactorisation a) => a -> Prime a precPrime n - | n < 2 = error $ "precPrime: tried to take `precPrime` of an argument less than 2" + | n < 2 = error "precPrime: tried to take `precPrime` of an argument less than 2" | n < 3 = Prime 2 | n < 5 = Prime 3 | n < 7 = Prime 5 @@ -153,7 +153,7 @@ Prime 2 -> Prime 3 Prime 3 -> Prime 5 Prime 5 -> Prime 7 - Prime p -> head $ mapMaybe isPrime $ map fromWheel30 [toWheel30 p + 1 ..] + Prime p -> head $ mapMaybe (isPrime . fromWheel30) [toWheel30 p + 1 ..] succGenericBounded :: (Bits a, Integral a, UniqueFactorisation a, Bounded a) @@ -163,7 +163,7 @@ Prime 2 -> Prime 3 Prime 3 -> Prime 5 Prime 5 -> Prime 7 - Prime p -> case mapMaybe isPrime $ map fromWheel30 [toWheel30 p + 1 .. toWheel30 maxBound] of + Prime p -> case mapMaybe (isPrime . fromWheel30) [toWheel30 p + 1 .. toWheel30 maxBound] of [] -> error "Enum.succ{Prime}: tried to take `succ' near `maxBound'" q : _ -> q @@ -173,7 +173,7 @@ Prime 3 -> Prime 2 Prime 5 -> Prime 3 Prime 7 -> Prime 5 - Prime p -> head $ mapMaybe isPrime $ map fromWheel30 [toWheel30 p - 1, toWheel30 p - 2 ..] + Prime p -> head $ mapMaybe (isPrime . fromWheel30) [toWheel30 p - 1, toWheel30 p - 2 ..] -- 'dropWhile' is important, because 'psieveFrom' can actually contain primes less than p. enumFromGeneric :: Integral a => Prime a -> [Prime a] @@ -206,7 +206,7 @@ -> [Prime a] enumFromToGeneric' p@(Prime p') q@(Prime q') = takeWhile (<= q) $ dropWhile (< p) $ case chooseAlgorithm p' q' of - IsPrime -> Prime 2 : Prime 3 : Prime 5 : mapMaybe isPrime (map fromWheel30 [toWheel30 p' .. toWheel30 q']) + IsPrime -> Prime 2 : Prime 3 : Prime 5 : mapMaybe (isPrime . fromWheel30) [toWheel30 p' .. toWheel30 q'] Sieve -> if q' < fromIntegral sieveRange then primeList $ primeSieve $ toInteger q' 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 @@ -7,7 +7,6 @@ -- Number of primes not exceeding @n@, @π(n)@, and @n@-th prime. -- {-# LANGUAGE BangPatterns #-} -{-# LANGUAGE CPP #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE ScopedTypeVariables #-} @@ -18,14 +17,13 @@ , nthPrime ) where -#include "MachDeps.h" - import Math.NumberTheory.Primes.Sieve.Eratosthenes (PrimeSieve(..), primeList, primeSieve, psieveFrom, sieveTo, sieveBits, sieveRange) import Math.NumberTheory.Primes.Sieve.Indexing (toPrim, idxPr) import Math.NumberTheory.Primes.Counting.Approximate (nthPrimeApprox, approxPrimeCount) import Math.NumberTheory.Primes.Types import Math.NumberTheory.Roots +import Math.NumberTheory.Utils.FromIntegral import Control.Monad.ST import Data.Array.Base @@ -52,12 +50,12 @@ primeCount n | n > primeCountMaxArg = error $ "primeCount: can't handle bound " ++ show n | n < 2 = 0 - | n < 1000 = fromIntegral . length . takeWhile (<= n) . map unPrime . primeList . primeSieve $ max 242 n + | n < 1000 = intToInteger . length . takeWhile (<= n) . map unPrime . primeList . primeSieve $ max 242 n | n < 30000 = runST $ do ba <- sieveTo n (s,e) <- getBounds ba ct <- countFromTo s e ba - return (fromIntegral $ ct+3) + return (intToInteger $ ct+3) | otherwise = let !ub = cop $ fromInteger n !sr = integerSquareRoot ub @@ -114,7 +112,7 @@ | otherwise = lowSieve p0 shortage -- a third count wouldn't make it faster, I think where - gap = truncate (log (fromIntegral p0 :: Double)) + gap = truncate (log (intToDouble p0 :: Double)) est = toInteger shortage * gap p1 = toInteger p0 + est goodEnough = 3*est*est*est < 2*p1*p1 -- a second counting would be more work than sieving @@ -127,7 +125,7 @@ | otherwise = tooHigh n b (c-n) where - gap = truncate (log (fromIntegral p0 :: Double)) + gap = truncate (log (intToDouble p0 :: Double)) b = p0 - (surplus * gap * 11) `quot` 10 c = fromInteger (primeCount (toInteger b)) @@ -152,11 +150,11 @@ sieveCountST :: forall s. Int64 -> Int64 -> Int64 -> ST s Integer sieveCountST ub cr sr = do - let psieves = psieveFrom (fromIntegral cr) + let psieves = psieveFrom (int64ToInteger cr) pisr = approxPrimeCount sr picr = approxPrimeCount cr diff = pisr - picr - size = fromIntegral (diff + diff `quot` 50) + 30 + 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 feed voff !wi !ri uar sves @@ -184,28 +182,28 @@ | otherwise = do qb <- unsafeRead store wi let dist = qb - voff - 7 - if dist < fromIntegral sieveRange + if dist < intToInt64 sieveRange then do let (b,j) = idxPr (dist+7) !li = (b `shiftL` 3) .|. j new <- if li < si then return 0 else countFromTo si li stu - let nbtw = btw + fromIntegral new + 1 + let nbtw = btw + intToInteger new + 1 eat (acc+nbtw) nbtw voff (wi-1) (li+1) stu sves else do - let (cpl,fds) = dist `quotRem` fromIntegral sieveRange + let (cpl,fds) = dist `quotRem` intToInt64 sieveRange (b,j) = idxPr (fds+7) !li = (b `shiftL` 3) .|. j ctLoop !lac 0 (PS vO ba : more) = do nstu <- unsafeThaw ba new <- countFromTo 0 li nstu - let nbtw = btw + lac + 1 + fromIntegral new - eat (acc+nbtw) nbtw (fromIntegral vO) (wi-1) (li+1) nstu more + let nbtw = btw + lac + 1 + intToInteger new + eat (acc+nbtw) nbtw (integerToInt64 vO) (wi-1) (li+1) nstu more ctLoop lac s (ps : more) = do let !new = countAll ps - ctLoop (lac + fromIntegral new) (s-1) more + ctLoop (lac + intToInteger new) (s-1) more ctLoop _ _ [] = error "Primes ended" new <- countFromTo si (sieveBits-1) stu - ctLoop (fromIntegral new) (cpl-1) sves + ctLoop (intToInteger new) (cpl-1) sves case psieves of (PS vO ba : more) -> feed (fromInteger vO) 0 0 ba more _ -> error "No primes sieved" @@ -215,7 +213,7 @@ calcST :: forall s. Int64 -> Int64 -> ST s Integer calcST lim plim = do - !parr <- sieveTo (fromIntegral plim) + !parr <- sieveTo (int64ToInteger plim) (plo,phi) <- getBounds parr !pct <- countFromTo plo phi parr !ar1 <- unsafeNewArray_ (0,end-1) @@ -239,14 +237,14 @@ | i < stop = do !k <- unsafeRead ar i !v <- unsafeRead ar (i+1) - cgo (acc + fromIntegral v*cp6 k) (i+2) - | otherwise = return (acc+fromIntegral pct+2) + cgo (acc + int64ToInteger v*cp6 k) (i+2) + | otherwise = return (acc+intToInteger pct+2) in cgo 0 0 go 2 start ar1 ar2 where (bt,ri) = idxPr plim !start = 8*bt + ri - !size = fromIntegral $ (integerSquareRoot lim) `quot` 4 + !size = int64ToInt $ integerSquareRoot lim `quot` 4 !end = 2*size treat :: Int -> Int64 -> STUArray s Int Int64 -> STUArray s Int Int64 -> ST s Int @@ -266,9 +264,9 @@ !key <- unsafeRead old qi !val <- unsafeRead old (qi+1) let !q0 = key `quot` n - !r0 = fromIntegral (q0 `rem` 30030) - !nkey = q0 - fromIntegral (cpDfAr `unsafeAt` r0) - nk0 = q0 + fromIntegral (cpGpAr `unsafeAt` (r0+1) + 1) + !r0 = int64ToInt (q0 `rem` 30030) + !nkey = q0 - int8ToInt64 (cpDfAr `unsafeAt` r0) + nk0 = q0 + int8ToInt64 (cpGpAr `unsafeAt` (r0+1) + 1) !nlim = n*nk0 (wi1,ci1) <- copyTo end nkey old ci new wi ckey <- unsafeRead old ci1 @@ -332,11 +330,11 @@ cp6 :: Int64 -> Integer cp6 k = case k `quotRem` 30030 of - (q,r) -> 5760*fromIntegral q + - fromIntegral (cpCtAr `unsafeAt` fromIntegral r) + (q,r) -> 5760*int64ToInteger q + + int16ToInteger (cpCtAr `unsafeAt` int64ToInt r) cop :: Int64 -> Int64 -cop m = m - fromIntegral (cpDfAr `unsafeAt` fromIntegral (m `rem` 30030)) +cop m = m - int8ToInt64 (cpDfAr `unsafeAt` int64ToInt (m `rem` 30030)) -------------------------------------------------------------------------------- @@ -413,21 +411,17 @@ ------------------------------------------------------------------------------- -- Prime counting -#if SIZEOF_HSWORD == 8 - -#define RMASK 63 -#define WSHFT 6 -#define TOPB 32 -#define TOPM 0xFFFFFFFF +rMASK :: Int +rMASK = finiteBitSize (0 :: Word) - 1 -#else +wSHFT :: (Bits a, Num a) => a +wSHFT = if finiteBitSize (0 :: Word) == 64 then 6 else 5 -#define RMASK 31 -#define WSHFT 5 -#define TOPB 16 -#define TOPM 0xFFFF +tOPB :: Int +tOPB = finiteBitSize (0 :: Word) `shiftR` 1 -#endif +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 @@ -448,7 +442,7 @@ then go (k-bc) (i+1) else let j = bc - k px = top w j bc - in v0 + toPrim (px + (i `shiftL` WSHFT)) + in v0 + toPrim (px + (i `shiftL` wSHFT)) -- count all set bits in a chunk, do it wordwise for speed. countAll :: PrimeSieve -> Int @@ -465,7 +459,7 @@ -- 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 +top w j bc = go 0 tOPB tOPM bn w where !bn = bc-j go !_ _ !_ !_ 0 = error "Too few bits set" @@ -482,14 +476,14 @@ 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 + 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))) + return (acc + popCount (w `shiftL` (rMASK - ei))) | otherwise = do w <- unsafeRead wa i count (acc + popCount w) (i+1) @@ -500,4 +494,4 @@ else do w <- unsafeRead wa sb let !w1 = w `shiftR` si - return (popCount (w1 `shiftL` (RMASK - ei + si))) + return (popCount (w1 `shiftL` (rMASK - ei + si))) 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 @@ -69,16 +69,13 @@ -- an arbitrary manner from the bit-pattern of @n@. -- -- __Warning:__ there are no guarantees of any particular --- order of prime factors, do not expect them to be ascending. E. g., --- --- >>> factorise 10251562501 --- [(101701,1),(100801,1)] +-- order of prime factors, do not expect them to be ascending. factorise :: Integral a => a -> [(a, Word)] factorise 0 = error "0 has no prime factorisation" factorise n' = map (first fromIntegral) sfs <> map (first fromInteger) rest where n = abs n' - (sfs, mb) = smallFactors (fromIntegral n) + (sfs, mb) = smallFactors (fromIntegral' n) sg = mkStdGen (fromIntegral n `xor` 0xdeadbeef) rest = case mb of Nothing -> [] @@ -97,8 +94,8 @@ -> Maybe Int -- ^ Estimated number of digits of smallest prime factor -> Integer -- ^ The number to factorise -> [(Integer, Word)] -- ^ List of prime factors and exponents -stdGenFactorisation primeBound sg digits n - = curveFactorisation primeBound bailliePSW (\m -> randomR (6,m-2)) sg digits n +stdGenFactorisation primeBound = + curveFactorisation primeBound bailliePSW (\m -> randomR (6, m - 2)) -- | 'curveFactorisation' is the driver for the factorisation. Its performance (and success) -- can be influenced by passing appropriate arguments. If you know that @n@ has no prime divisors @@ -181,7 +178,7 @@ -- Since all @cs@ are coprime, we can factor each of -- them and just concat results, without summing up -- powers of the same primes in different elements. - fmap mconcat $ flip mapM cs $ + fmap mconcat $ forM cs $ \(x, xm) -> if ptest x then pure $ singlePrimeFactor x xm else repFact x b1 b2 (count - 1) @@ -307,7 +304,7 @@ n = pointN q b0 = b1 - b1 `rem` wheel - qks = zip [0..] $ map (\k -> multiply k q) wheelCoprimes + qks = zip [0..] $ map (`multiply` q) wheelCoprimes qs = enumAndMultiplyFromThenTo q b0 (b0 + wheel) b2 rs = foldl' (\ts (_cHi, p) -> foldl' (\us (_cLo, pq) -> @@ -337,7 +334,7 @@ pThen = multiply thn p pStep = multiply step p - progression = pFrom : pThen : zipWith (\x0 x1 -> add x0 pStep x1) progression (tail progression) + progression = pFrom : pThen : zipWith (`add` pStep) progression (tail progression) -- primes, compactly stored as a bit sieve primeStore :: [PrimeSieve] @@ -365,7 +362,7 @@ !(Ptr smallPrimesAddr#) = smallPrimesPtr goBigNat :: BigNat -> Int -> ([(Natural, Word)], Maybe Natural) - goBigNat !m !i@(I# i#) + goBigNat !m i@(I# i#) | isTrue# (sizeofBigNat# m ==# 1#) = goWord (bigNatToWord m) i | i >= smallPrimesLength @@ -389,7 +386,7 @@ = if isTrue# (m# `leWord#` 4294967295##) -- 65536 * 65536 - 1 then ([(NatS# m#, 1)], Nothing) else ([], Just (NatS# m#)) - goWord m# !i@(I# i#) = let p# = indexWord16OffAddr# smallPrimesAddr# i# in + goWord m# i@(I# i#) = let p# = indexWord16OffAddr# smallPrimesAddr# i# in if isTrue# (m# `ltWord#` (p# `timesWord#` p#)) then ([(NatS# m#, 1)], Nothing) else case m# `quotRemWord#` p# of diff --git a/Math/NumberTheory/Primes/Factorisation/TrialDivision.hs b/Math/NumberTheory/Primes/Factorisation/TrialDivision.hs --- a/Math/NumberTheory/Primes/Factorisation/TrialDivision.hs +++ b/Math/NumberTheory/Primes/Factorisation/TrialDivision.hs @@ -49,7 +49,7 @@ trialDivisionTo bd | bd < 100 = trialDivisionTo 100 | bd < 10000000 = trialDivisionWith (map unPrime $ primeList $ primeSieve bd) - | otherwise = trialDivisionWith (takeWhile (<= bd) $ map unPrime $ (psieveList >>= primeList)) + | otherwise = trialDivisionWith (takeWhile (<= bd) $ map unPrime $ psieveList >>= primeList) -- | Check whether a number is coprime to all of the numbers in the list -- (assuming that list contains only numbers > 1 and is ascending). @@ -68,4 +68,4 @@ trialDivisionPrimeTo bd | bd < 100 = trialDivisionPrimeTo 100 | bd < 10000000 = trialDivisionPrimeWith (map unPrime $ primeList $ primeSieve bd) - | otherwise = trialDivisionPrimeWith (takeWhile (<= bd) $ map unPrime $ (psieveList >>= primeList)) + | otherwise = trialDivisionPrimeWith (takeWhile (<= bd) $ map unPrime $ psieveList >>= primeList) diff --git a/Math/NumberTheory/Primes/IntSet.hs b/Math/NumberTheory/Primes/IntSet.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Primes/IntSet.hs @@ -0,0 +1,338 @@ +-- | +-- Module: Math.NumberTheory.Primes.IntSet +-- Copyright: (c) 2020 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- +-- A newtype wrapper around 'IntSet'. +-- +-- This module is intended to be imported qualified, e. g., +-- +-- > import Math.NumberTheory.Primes.IntSet (PrimeIntSet) +-- > import qualified Math.NumberTheory.Primes.IntSet as PrimeIntSet +-- + +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE DeriveDataTypeable #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeApplications #-} +{-# LANGUAGE TypeFamilies #-} + +module Math.NumberTheory.Primes.IntSet + ( -- * Set type + PrimeIntSet + , unPrimeIntSet + -- * Construction + -- | Use 'Data.Monoid.mempty' to create an empty set. + , singleton + , fromList + , fromAscList + , fromDistinctAscList + -- * Insertion + , insert + -- * Deletion + , delete + -- * Query + , member + , notMember + , lookupEQ + , lookupLT + , lookupGT + , lookupLE + , lookupGE + , null + , size + , isSubsetOf + , isProperSubsetOf + , disjoint + -- * Combine + -- | Use 'Data.Semigroup.<>' for unions. + , difference + , (\\) + , symmetricDifference + , intersection + -- * Filter + , filter + , partition + , split + , splitMember + , splitLookupEQ + , splitRoot + -- * Folds + , foldr + , foldl + , foldr' + , foldl' + -- * Min/Max + , deleteMin + , deleteMax + , minView + , maxView + -- * Conversion + , toAscList + , toDescList + ) where + +import Prelude ((>), (/=), (==), (-), Eq, Ord, Show, Monoid, Bool, Maybe(..), Int, Word, otherwise) +import Control.DeepSeq (NFData) +import Data.Coerce (coerce) +import Data.Data (Data) +import Data.Function (on) +import Data.IntSet (IntSet) +import qualified Data.IntSet.Internal as IS +import Data.Semigroup (Semigroup) +import qualified GHC.Exts (IsList(..)) + +import Math.NumberTheory.Primes.Types (Prime(..)) +import Math.NumberTheory.Utils.FromIntegral (wordToInt, intToWord) +import Data.Bits (Bits(..)) +import Utils.Containers.Internal.BitUtil (highestBitMask) + +-- | A set of 'Prime' integers. +newtype PrimeIntSet = PrimeIntSet { + -- | Convert to a set of integers. + unPrimeIntSet :: IntSet + } + deriving (Eq, Ord, Data, Show, Semigroup, Monoid, NFData) + +instance GHC.Exts.IsList PrimeIntSet where + type Item PrimeIntSet = Prime Int + fromList = coerce IS.fromList + toList = coerce IS.toList + +-- | Build a singleton set. +singleton :: Prime Int -> PrimeIntSet +singleton = coerce IS.singleton + +-- | Build a set from a list of primes. +fromList :: [Prime Int] -> PrimeIntSet +fromList = coerce IS.fromList + +-- | Build a set from an ascending list of primes +-- (the precondition is not checked). +fromAscList :: [Prime Int] -> PrimeIntSet +fromAscList = coerce IS.fromAscList + +-- | Build a set from an ascending list of distinct primes +-- (the precondition is not checked). +fromDistinctAscList :: [Prime Int] -> PrimeIntSet +fromDistinctAscList = coerce IS.fromDistinctAscList + +-- | Insert a prime into the set. +insert :: Prime Int -> PrimeIntSet -> PrimeIntSet +insert = coerce IS.insert + +-- | Delete an integer from the set. +delete :: Int -> PrimeIntSet -> PrimeIntSet +delete = coerce IS.delete + +-- | Check whether the given prime is a member of the set. +member :: Prime Int -> PrimeIntSet -> Bool +member = coerce IS.member + +-- | Check whether the given prime is not a member of the set. +notMember :: Prime Int -> PrimeIntSet -> Bool +notMember = coerce IS.notMember + +-- | Find a prime in the set, +-- equal to the given integer, if any exists. +lookupEQ :: Int -> PrimeIntSet -> Maybe (Prime Int) +lookupEQ x xs + | coerce member x xs = Just (Prime x) + | otherwise = Nothing + +-- | Find the largest prime in the set, +-- smaller than the given integer, if any exists. +lookupLT :: Int -> PrimeIntSet -> Maybe (Prime Int) +lookupLT = coerce IS.lookupLT + +-- | Find the smallest prime in the set, +-- greater than the given integer, if any exists. +lookupGT :: Int -> PrimeIntSet -> Maybe (Prime Int) +lookupGT = coerce IS.lookupGT + +-- | Find the largest prime in the set, +-- smaller or equal to the given integer, if any exists. +lookupLE :: Int -> PrimeIntSet -> Maybe (Prime Int) +lookupLE = coerce IS.lookupLE + +-- | Find the smallest prime in the set, +-- greater or equal to the given integer, if any exists. +lookupGE :: Int -> PrimeIntSet -> Maybe (Prime Int) +lookupGE = coerce IS.lookupGE + +-- | Check whether the set is empty. +null :: PrimeIntSet -> Bool +null = coerce IS.null + +-- | Cardinality of the set. +size :: PrimeIntSet -> Int +size = coerce IS.size + +-- | Check whether the first argument is a subset of the second one. +isSubsetOf :: PrimeIntSet -> PrimeIntSet -> Bool +isSubsetOf = coerce IS.isSubsetOf + +-- | Check whether the first argument is a proper subset of the second one. +isProperSubsetOf :: PrimeIntSet -> PrimeIntSet -> Bool +isProperSubsetOf = coerce IS.isProperSubsetOf + +#if MIN_VERSION_containers(0,5,11) +-- | Check whether two sets are disjoint. +disjoint :: PrimeIntSet -> PrimeIntSet -> Bool +disjoint = coerce IS.disjoint +#else +-- | Check whether two sets are disjoint. +disjoint :: PrimeIntSet -> PrimeIntSet -> Bool +disjoint (PrimeIntSet x) (PrimeIntSet y) = IS.null (IS.intersection x y) +#endif + +-- | Difference between a set of primes and a set of integers. +difference :: PrimeIntSet -> IntSet -> PrimeIntSet +difference = coerce IS.difference + +-- | An alias to 'difference'. +(\\) :: PrimeIntSet -> IntSet -> PrimeIntSet +(\\) = coerce (IS.\\) + +infixl 9 \\{- -} + +-- | Symmetric difference of two sets of primes. +symmetricDifference :: PrimeIntSet -> PrimeIntSet -> PrimeIntSet +symmetricDifference = coerce symmDiff + +-- | Intersection of a set of primes and a set of integers. +intersection :: PrimeIntSet -> IntSet -> PrimeIntSet +intersection = coerce IS.intersection + +-- | Filter primes satisfying a predicate. +filter :: (Prime Int -> Bool) -> PrimeIntSet -> PrimeIntSet +filter = coerce IS.filter + +-- | Partition primes according to a predicate. +partition :: (Prime Int -> Bool) -> PrimeIntSet -> (PrimeIntSet, PrimeIntSet) +partition = coerce IS.partition + +-- | Split into primes strictly less and strictly greater +-- than the first argument. +split :: Int -> PrimeIntSet -> (PrimeIntSet, PrimeIntSet) +split = coerce IS.split + +-- | Simultaneous 'split' and 'member'. +splitMember :: Prime Int -> PrimeIntSet -> (PrimeIntSet, Bool, PrimeIntSet) +splitMember = coerce IS.splitMember + +-- | Simultaneous 'split' and 'lookupEQ'. +splitLookupEQ :: Int -> PrimeIntSet -> (PrimeIntSet, Maybe (Prime Int), PrimeIntSet) +splitLookupEQ x xs = (lt, if eq then Just (Prime x) else Nothing, gt) + where + (lt, eq, gt) = coerce IS.splitMember x xs + +-- | Decompose a set into pieces based on the structure of the underlying tree. +splitRoot :: PrimeIntSet -> [PrimeIntSet] +splitRoot = coerce IS.splitRoot + +-- | Fold a set using the given right-associative operator. +foldr :: forall b. (Prime Int -> b -> b) -> b -> PrimeIntSet -> b +foldr = coerce (IS.foldr @b) + +-- | Fold a set using the given left-associative operator. +foldl :: forall a. (a -> Prime Int -> a) -> a -> PrimeIntSet -> a +foldl = coerce (IS.foldl @a) + +-- | A strict version of 'foldr'. +foldr' :: forall b. (Prime Int -> b -> b) -> b -> PrimeIntSet -> b +foldr' = coerce (IS.foldr' @b) + +-- | A strict version of 'foldl'. +foldl' :: forall a. (a -> Prime Int -> a) -> a -> PrimeIntSet -> a +foldl' = coerce (IS.foldl' @a) + +-- | Delete the smallest prime in the set. +deleteMin :: PrimeIntSet -> PrimeIntSet +deleteMin = coerce IS.deleteMin + +-- | Delete the largest prime in the set. +deleteMax :: PrimeIntSet -> PrimeIntSet +deleteMax = coerce IS.deleteMax + +-- | Split a set into the smallest prime and the rest, if non-empty. +minView :: PrimeIntSet -> Maybe (Prime Int, PrimeIntSet) +minView = coerce IS.minView + +-- | Split a set into the largest prime and the rest, if non-empty. +maxView :: PrimeIntSet -> Maybe (Prime Int, PrimeIntSet) +maxView = coerce IS.maxView + +-- | Convert the set to a list of ascending primes. +toAscList :: PrimeIntSet -> [Prime Int] +toAscList = coerce IS.toAscList + +-- | Convert the set to a list of descending primes. +toDescList :: PrimeIntSet -> [Prime Int] +toDescList = coerce IS.toDescList + +------------------------------------------------------------------------------- +-- IntSet helpers + +-- | Symmetric difference of two sets. +-- Implementation is inspired by 'Data.IntSet.union' +-- and 'Data.IntSet.difference'. +symmDiff :: IntSet -> IntSet -> IntSet +symmDiff t1 t2 = case t1 of + IS.Bin p1 m1 l1 r1 -> case t2 of + IS.Bin p2 m2 l2 r2 + | shorter m1 m2 -> symmDiff1 + | shorter m2 m1 -> symmDiff2 + | p1 == p2 -> bin p1 m1 (symmDiff l1 l2) (symmDiff r1 r2) + | otherwise -> link p1 t1 p2 t2 + where + symmDiff1 + | mask p2 m1 /= p1 = link p1 t1 p2 t2 + | p2 .&. m1 == 0 = bin p1 m1 (symmDiff l1 t2) r1 + | otherwise = bin p1 m1 l1 (symmDiff r1 t2) + symmDiff2 + | mask p1 m2 /= p2 = link p1 t1 p2 t2 + | p1 .&. m2 == 0 = bin p2 m2 (symmDiff t1 l2) r2 + | otherwise = bin p2 m2 l2 (symmDiff t1 r2) + IS.Tip kx bm -> symmDiffBM kx bm t1 + IS.Nil -> t1 + IS.Tip kx bm -> symmDiffBM kx bm t2 + IS.Nil -> t2 + +shorter :: Int -> Int -> Bool +shorter = (>) `on` intToWord + +symmDiffBM :: Int -> Word -> IntSet -> IntSet +symmDiffBM !kx !bm t = case t of + IS.Bin p m l r + | mask kx m /= p -> link kx (IS.Tip kx bm) p t + | kx .&. m == 0 -> bin p m (symmDiffBM kx bm l) r + | otherwise -> bin p m l (symmDiffBM kx bm r) + IS.Tip kx' bm' + | kx' == kx -> if bm' == bm then IS.Nil else IS.Tip kx (bm' `xor` bm) + | otherwise -> link kx (IS.Tip kx bm) kx' t + IS.Nil -> IS.Tip kx bm + +link :: Int -> IntSet -> Int -> IntSet -> IntSet +link p1 t1 p2 t2 + | p1 .&. m == 0 = IS.Bin p m t1 t2 + | otherwise = IS.Bin p m t2 t1 + where + m = wordToInt (highestBitMask (intToWord p1 `xor` intToWord p2)) + p = mask p1 m +{-# INLINE link #-} + +bin :: Int -> Int -> IntSet -> IntSet -> IntSet +bin p m l r = case r of + IS.Nil -> l + _ -> case l of + IS.Nil -> r + _ -> IS.Bin p m l r +{-# INLINE bin #-} + +mask :: Int -> Int -> Int +mask i m = i .&. (complement (m - 1) `xor` m) +{-# INLINE mask #-} 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 @@ -7,7 +7,6 @@ -- Sieve -- {-# LANGUAGE BangPatterns #-} -{-# LANGUAGE CPP #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE LambdaCase #-} {-# LANGUAGE ScopedTypeVariables #-} @@ -25,8 +24,6 @@ , sieveTo ) where -#include "MachDeps.h" - import Control.Monad (when) import Control.Monad.ST import Data.Array.Base @@ -39,24 +36,35 @@ import Math.NumberTheory.Primes.Sieve.Indexing import Math.NumberTheory.Primes.Types import Math.NumberTheory.Roots +import Math.NumberTheory.Utils.FromIntegral -#define IX_MASK 0xFFFFF -#define IX_BITS 20 -#define IX_J_MASK 0x7FFFFF -#define IX_J_BITS 23 -#define J_MASK 7 -#define J_BITS 3 -#define SIEVE_KB 128 +iXMASK :: Num a => a +iXMASK = 0xFFFFF +iXBITS :: Int +iXBITS = 20 + +iXJMASK :: Num a => a +iXJMASK = 0x7FFFFF + +iXJBITS :: Int +iXJBITS = 23 + +jMASK :: Int +jMASK = 7 + +jBITS :: Int +jBITS = 3 + -- Sieve in 128K chunks. -- Large enough to get something done per chunk -- and hopefully small enough to fit in the cache. sieveBytes :: Int -sieveBytes = SIEVE_KB*1024 +sieveBytes = 128 * 1024 -- Number of bits per chunk. sieveBits :: Int -sieveBits = 8*sieveBytes +sieveBits = 8 * sieveBytes -- Last index of chunk. lastIndex :: Int @@ -64,21 +72,10 @@ -- Range of a chunk. sieveRange :: Int -sieveRange = 30*sieveBytes - -type CacheWord = Word64 +sieveRange = 30 * sieveBytes -#if SIZEOF_HSWORD == 8 -#define RMASK 63 -#define WSHFT 6 -#define TOPB 32 -#define TOPM 0xFFFFFFFF -#else -#define RMASK 31 -#define WSHFT 5 -#define TOPB 16 -#define TOPM 0xFFFF -#endif +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) @@ -135,23 +132,18 @@ -- 'primes' is a polymorphic list, so the results of computations are not retained in memory. -- Make it monomorphic to take advantages of memoization. Compare -- --- >>> :set +s --- >>> primes !! 1000000 :: Prime Int +-- >>> primes !! 1000000 :: Prime Int -- (5.32 secs, 6,945,267,496 bytes) -- Prime 15485867 --- (5.32 secs, 6,945,267,496 bytes) --- >>> primes !! 1000000 :: Prime Int +-- >>> primes !! 1000000 :: Prime Int -- (5.19 secs, 6,945,267,496 bytes) -- Prime 15485867 --- (5.19 secs, 6,945,267,496 bytes) -- -- against -- -- >>> let primes' = primes :: [Prime Int] --- >>> primes' !! 1000000 :: Prime Int +-- >>> primes' !! 1000000 :: Prime Int -- (5.29 secs, 6,945,269,856 bytes) -- Prime 15485867 --- (5.29 secs, 6,945,269,856 bytes) --- >>> primes' !! 1000000 :: Prime Int +-- >>> primes' !! 1000000 :: Prime Int -- (0.02 secs, 336,232 bytes) -- Prime 15485867 --- (0.02 secs, 336,232 bytes) primes :: Integral a => [Prime a] primes = (coerce :: [a] -> [Prime a]) @@ -167,30 +159,30 @@ sqlim = plim*plim cache = runSTUArray $ do sieve <- sieveTo (4801 :: Integer) - new <- unsafeNewArray_ (0,1287) :: ST s (STUArray s Int CacheWord) + new <- unsafeNewArray_ (0,1287) :: ST s (STUArray s Int Word64) let fill j indx | 1279 < indx = return new -- index of 4801 = 159*30 + 31 ~> 159*8+7 | otherwise = do p <- unsafeRead sieve indx if p then do - let !i = indx .&. J_MASK - k = indx `shiftR` J_BITS - strt1 = (k*(30*k + 2*rho i) + byte i) `shiftL` J_BITS + idx i - !strt = fromIntegral (strt1 .&. IX_MASK) - !skip = fromIntegral (strt1 `shiftR` IX_BITS) - !ixes = fromIntegral indx `shiftL` IX_J_BITS + strt `shiftL` J_BITS + fromIntegral i + let !i = indx .&. jMASK + k = indx `shiftR` jBITS + strt1 = (k*(30*k + 2*rho i) + byte i) `shiftL` jBITS + idx i + !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 fill (j+2) (indx+1) else fill j (indx+1) fill 0 0 -makeSieves :: Integer -> Integer -> Integer -> Integer -> UArray Int CacheWord -> [PrimeSieve] +makeSieves :: Integer -> Integer -> Integer -> Integer -> UArray Int Word64 -> [PrimeSieve] makeSieves plim sqlim bitOff valOff cache | valOff' < sqlim = let (nc, bs) = runST $ do - cch <- unsafeThaw cache :: ST s (STUArray s Int CacheWord) + cch <- unsafeThaw cache :: ST s (STUArray s Int Word64) bs0 <- slice cch fcch <- unsafeFreeze cch fbs0 <- unsafeFreeze bs0 @@ -207,10 +199,10 @@ return (fcch, fbs0) in PS valOff bs : makeSieves plim' sqlim' bitOff' valOff' nc where - valOff' = valOff + fromIntegral sieveRange - bitOff' = bitOff + fromIntegral sieveBits + valOff' = valOff + intToInteger sieveRange + bitOff' = bitOff + intToInteger sieveBits -slice :: STUArray s Int CacheWord -> ST s (STUArray s Int Bool) +slice :: STUArray s Int Word64 -> ST s (STUArray s Int Bool) slice cache = do hi <- snd `fmap` getBounds cache sieve <- newArray (0,lastIndex) True @@ -222,25 +214,25 @@ then unsafeWrite cache pr (w-1) else do ixes <- unsafeRead cache (pr+1) - let !stj = fromIntegral ixes .&. IX_J_MASK -- position of multiple and index of cofactor - !ixw = fromIntegral (ixes `shiftR` IX_J_BITS) -- prime data, up to 41 bits - !i = ixw .&. J_MASK + 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 !k = ixw - i -- On 32-bits, k > 44717396 means overflow is possible in tick - !o = i `shiftL` J_BITS - !j = stj .&. J_MASK -- index of cofactor - !s = stj `shiftR` J_BITS -- index of first multiple to tick off + !o = i `shiftL` jBITS + !j = stj .&. jMASK -- index of cofactor + !s = stj `shiftR` jBITS -- index of first multiple to tick off (n, u) <- tick k o j s - let !skip = fromIntegral (n `shiftR` IX_BITS) - !strt = fromIntegral (n .&. IX_MASK) + let !skip = intToWord64 (n `shiftR` iXBITS) + !strt = intToWord64 (n .&. iXMASK) unsafeWrite cache pr skip - unsafeWrite cache (pr+1) ((ixes .&. complement IX_J_MASK) .|. strt `shiftL` J_BITS .|. fromIntegral u) + 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) - tick stp off ((j+1) .&. J_MASK) (ix + stp*delta j + tau (off+j)) + tick stp off ((j+1) .&. jMASK) (ix + stp*delta j + tau (off+j)) treat 0 -- | Sieve up to bound in one go. @@ -250,7 +242,7 @@ (bytes,lidx) = idxPr bound !mxidx = 8*bytes+lidx mxval :: Integer - mxval = 30*fromIntegral bytes + fromIntegral (rho lidx) + mxval = 30*intToInteger bytes + intToInteger (rho lidx) !mxsve = integerSquareRoot mxval (kr,r) = idxPr mxsve !svbd = 8*kr+r @@ -262,20 +254,20 @@ | otherwise = do p <- unsafeRead ar ix when p (unsafeWrite ar ix False) - tick stp off ((j+1) .&. J_MASK) (ix + stp*delta j + tau (off+j)) + 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 - when p (do let i = ix .&. J_MASK - k = ix `shiftR` J_BITS - !off = i `shiftL` J_BITS + when p (do let i = ix .&. jMASK + k = ix `shiftR` jBITS + !off = i `shiftL` jBITS !stp = ix - i tick stp off i (start k i)) sift (ix+1) sift 0 -growCache :: Integer -> Integer -> UArray Int CacheWord -> ST s (STUArray s Int CacheWord) +growCache :: Integer -> Integer -> UArray Int Word64 -> ST s (STUArray s Int Word64) growCache offset plim old = do let (_,num) = bounds old (bt,ix) = idxPr plim @@ -284,7 +276,7 @@ sieve <- sieveTo nlim -- Implement SieveFromTo for this, it's pretty wasteful when nlim isn't (_,hi) <- getBounds sieve -- very small anymore more <- countFromToWd start hi sieve - new <- unsafeNewArray_ (0,num+2*more) :: ST s (STUArray s Int CacheWord) + new <- unsafeNewArray_ (0,num+2*more) :: ST s (STUArray s Int Word64) let copy i | num < i = return () | otherwise = do @@ -297,16 +289,16 @@ p <- unsafeRead sieve indx if p then do - let !i = indx .&. J_MASK + let !i = indx .&. jMASK k :: Integer - k = fromIntegral (indx `shiftR` J_BITS) - strt0 = ((k*(30*k + fromIntegral (2*rho i)) - + fromIntegral (byte i)) `shiftL` J_BITS) - + fromIntegral (idx i) + k = intToInteger (indx `shiftR` jBITS) + strt0 = ((k*(30*k + intToInteger (2*rho i)) + + intToInteger (byte i)) `shiftL` jBITS) + + intToInteger (idx i) strt1 = strt0 - offset - !strt = fromIntegral strt1 .&. IX_MASK - !skip = fromIntegral (strt1 `shiftR` IX_BITS) - !ixes = fromIntegral indx `shiftL` IX_J_BITS .|. strt `shiftL` J_BITS .|. fromIntegral i + !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 fill (j+2) (indx+1) @@ -320,8 +312,8 @@ 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 + let !sb = start `shiftR` wSHFT + !eb = end `shiftR` wSHFT count !acc i | eb < i = return acc | otherwise = do @@ -342,7 +334,7 @@ bitOff = 8*k0 start = valOff+7 ssr = integerSquareRoot (start-1) + 1 - end1 = start - 6 + fromIntegral sieveRange + end1 = start - 6 + intToInteger sieveRange plim0 = integerSquareRoot end1 plim = plim0 + 4801 - (plim0 `rem` 4800) sqlim = plim*plim @@ -350,40 +342,40 @@ sieve <- sieveTo plim (lo,hi) <- getBounds sieve pct <- countFromToWd lo hi sieve - new <- unsafeNewArray_ (0,2*pct-1) :: ST s (STUArray s Int CacheWord) + new <- unsafeNewArray_ (0,2*pct-1) :: ST s (STUArray s Int Word64) let fill j indx | hi < indx = return new | otherwise = do isPr <- unsafeRead sieve indx if isPr then do - let !i = indx .&. J_MASK - !moff = i `shiftL` J_BITS + let !i = indx .&. jMASK + !moff = i `shiftL` jBITS k :: Integer - k = fromIntegral (indx `shiftR` J_BITS) - p = 30*k+fromIntegral (rho i) + k = intToInteger (indx `shiftR` jBITS) + p = 30*k+intToInteger (rho i) q0 = (start-1) `quot` p - (skp0,q1) = q0 `quotRem` fromIntegral sieveRange + (skp0,q1) = q0 `quotRem` intToInteger sieveRange (b0,r0) | q1 == 0 = (-1,6) | q1 < 7 = (-1,7) - | otherwise = idxPr (fromIntegral q1 :: Int) + | otherwise = idxPr (integerToInt q1 :: Int) (b1,r1) | r0 == 7 = (b0+1,0) | otherwise = (b0,r0+1) - b2 = skp0*fromIntegral sieveBytes + fromIntegral b1 - strt0 = ((k*(30*b2 + fromIntegral (rho r1)) - + b2 * fromIntegral (rho i) - + fromIntegral (mu (moff + r1))) `shiftL` J_BITS) - + fromIntegral (nu (moff + r1)) - strt1 = ((k*(30*k + fromIntegral (2*rho i)) - + fromIntegral (byte i)) `shiftL` J_BITS) - + fromIntegral (idx i) + b2 = skp0*intToInteger sieveBytes + intToInteger b1 + strt0 = ((k*(30*b2 + intToInteger (rho r1)) + + b2 * intToInteger (rho i) + + intToInteger (mu (moff + r1))) `shiftL` jBITS) + + intToInteger (nu (moff + r1)) + strt1 = ((k*(30*k + intToInteger (2*rho i)) + + intToInteger (byte i)) `shiftL` jBITS) + + intToInteger (idx i) (strt2,r2) | p < ssr = (strt0 - bitOff,r1) | otherwise = (strt1 - bitOff, i) - !strt = fromIntegral strt2 .&. IX_MASK - !skip = fromIntegral (strt2 `shiftR` IX_BITS) - !ixes = fromIntegral indx `shiftL` IX_J_BITS .|. strt `shiftL` J_BITS .|. fromIntegral r2 + !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 fill (j+2) (indx+1) @@ -393,14 +385,14 @@ {-# INLINE delta #-} delta :: Int -> Int -delta i = unsafeAt deltas i +delta = unsafeAt deltas deltas :: UArray Int Int deltas = listArray (0,7) [4,2,4,2,4,6,2,6] {-# INLINE tau #-} tau :: Int -> Int -tau i = unsafeAt taus i +tau = unsafeAt taus taus :: UArray Int Int taus = listArray (0,63) @@ -416,25 +408,25 @@ {-# INLINE byte #-} byte :: Int -> Int -byte i = unsafeAt startByte i +byte = unsafeAt startByte startByte :: UArray Int Int startByte = listArray (0,7) [1,3,5,9,11,17,27,31] {-# INLINE idx #-} idx :: Int -> Int -idx i = unsafeAt startIdx i +idx = unsafeAt startIdx startIdx :: UArray Int Int startIdx = listArray (0,7) [4,7,4,4,7,4,7,7] {-# INLINE mu #-} mu :: Int -> Int -mu i = unsafeAt mArr i +mu = unsafeAt mArr {-# INLINE nu #-} nu :: Int -> Int -nu i = unsafeAt nArr i +nu = unsafeAt nArr mArr :: UArray Int Int mArr = listArray (0,63) 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 @@ -21,8 +21,8 @@ | n0 < 7 = (0, 0) | otherwise = (fromIntegral bytes0, rm3) where - n = if (fromIntegral n0 .&. 1 == (1 :: Int)) - then n0 else (n0-1) + n = if fromIntegral n0 .&. 1 == (1 :: Int) + then n0 else n0 - 1 (bytes0,rm0) = (n-7) `quotRem` 30 rm1 = fromIntegral rm0 rm2 = rm1 `quot` 3 @@ -37,7 +37,7 @@ {-# INLINE rho #-} rho :: Int -> Int -rho i = unsafeAt residues i +rho = unsafeAt residues residues :: UArray Int Int residues = listArray (0,7) [7,11,13,17,19,23,29,31] diff --git a/Math/NumberTheory/Primes/Small.hs b/Math/NumberTheory/Primes/Small.hs --- a/Math/NumberTheory/Primes/Small.hs +++ b/Math/NumberTheory/Primes/Small.hs @@ -21,7 +21,7 @@ import GHC.Word smallPrimesFromTo :: Word16 -> Word16 -> [Word16] -smallPrimesFromTo !(W16# from#) !(W16# to#) = go k0# +smallPrimesFromTo (W16# from#) (W16# to#) = go k0# where !(Ptr smallPrimesAddr#) = smallPrimesPtr fromD# = word2Double# from# diff --git a/Math/NumberTheory/Primes/Testing/Certified.hs b/Math/NumberTheory/Primes/Testing/Certified.hs --- a/Math/NumberTheory/Primes/Testing/Certified.hs +++ b/Math/NumberTheory/Primes/Testing/Certified.hs @@ -155,7 +155,7 @@ where sr = integerSquareRoot n pbd = min 1000000 (sr+20) - prms = map unPrime $ primeList (primeSieve $ pbd) + prms = map unPrime $ primeList (primeSieve pbd) go a b afs (p:ps) | a > b = (a,afs,b) | otherwise = case splitOff p b of diff --git a/Math/NumberTheory/Primes/Testing/Probabilistic.hs b/Math/NumberTheory/Primes/Testing/Probabilistic.hs --- a/Math/NumberTheory/Primes/Testing/Probabilistic.hs +++ b/Math/NumberTheory/Primes/Testing/Probabilistic.hs @@ -7,7 +7,6 @@ -- Probabilistic primality tests, Miller-Rabin and Baillie-PSW. {-# LANGUAGE BangPatterns #-} -{-# LANGUAGE CPP #-} {-# LANGUAGE MagicHash #-} {-# LANGUAGE ScopedTypeVariables #-} @@ -20,8 +19,6 @@ , lucasTest ) where -#include "MachDeps.h" - import Data.Bits import Data.Mod import Data.Proxy @@ -30,8 +27,9 @@ import GHC.TypeNats (KnownNat, SomeNat(..), someNatVal) import Math.NumberTheory.Moduli.JacobiSymbol -import Math.NumberTheory.Utils +import Math.NumberTheory.Primes.Small import Math.NumberTheory.Roots +import Math.NumberTheory.Utils -- | @isPrime n@ tests whether @n@ is a prime (negative or positive). -- It is a combination of trial division and Baillie-PSW test. @@ -60,9 +58,17 @@ -- If @millerRabinV k n@ returns @False@ then @n@ is definitely composite. -- Otherwise @n@ may appear composite with probability @1/4^k@. millerRabinV :: Int -> Integer -> Bool -millerRabinV (I# k) n = case testPrimeInteger n k of - 0# -> False - _ -> True +millerRabinV k n + | n < 0 = millerRabinV k (-n) + | n < 2 = False + | n < 4 = True + | otherwise = go smallPrimes + where + go (p:ps) + | p*p > n = True + | otherwise = (n `rem` p /= 0) && go ps + go [] = all (isStrongFermatPP n) (take k smallPrimes) + smallPrimes = map toInteger $ smallPrimesFromTo minBound maxBound -- | @'isStrongFermatPP' n b@ tests whether non-negative @n@ is -- a strong Fermat probable prime for base @b@. @@ -168,7 +174,7 @@ -- n odd positive, n > abs q, index odd testLucas :: Integer -> Integer -> Integer -> (Integer, Integer, Integer) -testLucas n q (S# i#) = look (WORD_SIZE_IN_BITS - 2) +testLucas n q (S# i#) = look (finiteBitSize (0 :: Word) - 2) where j = I# i# look k @@ -191,14 +197,14 @@ s# = sizeofBigNat# bn# test j# = case indexBigNat# bn# j# of 0## -> test (j# -# 1#) - w# -> look (j# -# 1#) (W# w#) (WORD_SIZE_IN_BITS - 1) + w# -> look (j# -# 1#) (W# w#) (finiteBitSize (0 :: Word) - 1) look j# w i | testBit w i = go j# w (i - 1) 1 1 1 q | otherwise = look j# w (i-1) go k# w i un un1 vn qn | i < 0 = if isTrue# (k# <# 0#) then (un,vn,qn) - else go (k# -# 1#) (W# (indexBigNat# bn# k#)) (WORD_SIZE_IN_BITS - 1) un un1 vn qn + else go (k# -# 1#) (W# (indexBigNat# bn# k#)) (finiteBitSize (0 :: Word) - 1) un un1 vn qn | testBit w i = go k# w (i-1) u2n1 u2n2 v2n1 q2n1 | otherwise = go k# w (i-1) u2n u2n1 v2n q2n where diff --git a/Math/NumberTheory/Primes/Types.hs b/Math/NumberTheory/Primes/Types.hs --- a/Math/NumberTheory/Primes/Types.hs +++ b/Math/NumberTheory/Primes/Types.hs @@ -8,17 +8,28 @@ -- Should not be exposed to users. -- -{-# LANGUAGE CPP #-} -{-# LANGUAGE TypeFamilies #-} -{-# LANGUAGE DeriveGeneric #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE DeriveGeneric #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MultiParamTypeClasses #-} module Math.NumberTheory.Primes.Types ( Prime(..) + , toPrimeIntegral ) where +import Data.Bits import GHC.Generics import Control.DeepSeq +import qualified Data.Vector.Generic as G +import qualified Data.Vector.Generic.Mutable as M +import qualified Data.Vector.Unboxed as U +import Math.NumberTheory.Utils.FromIntegral + +-- $setup +-- >>> import Math.NumberTheory.Primes (nextPrime, precPrime) + -- | Wrapper for prime elements of @a@. It is supposed to be constructed -- by 'Math.NumberTheory.Primes.nextPrime' / 'Math.NumberTheory.Primes.precPrime'. -- and eliminated by 'unPrime'. @@ -28,13 +39,14 @@ -- -- * Generate primes from the given interval: -- +-- >>> :set -XFlexibleContexts -- >>> [nextPrime 101 .. precPrime 130] -- [Prime 101,Prime 103,Prime 107,Prime 109,Prime 113,Prime 127] -- -- * Generate an infinite list of primes: -- --- >>> [nextPrime 101 ..] --- [Prime 101,Prime 103,Prime 107,Prime 109,Prime 113,Prime 127... +-- > [nextPrime 101 ..] +-- > [Prime 101,Prime 103,Prime 107,Prime 109,Prime 113,Prime 127... -- -- * Generate primes from the given interval of form p = 6k+5: -- @@ -48,7 +60,7 @@ -- -- * Get previous prime: -- --- >>> prec (nextPrime 101) +-- >>> pred (nextPrime 101) -- Prime 97 -- -- * Count primes less than a given number (cf. 'Math.NumberTheory.Primes.Counting.approxPrimeCount'): @@ -72,3 +84,86 @@ showsPrec d (Prime p) r = (if d > 10 then "(" ++ s ++ ")" else s) ++ r where s = "Prime " ++ show p + +newtype instance U.MVector s (Prime a) = MV_Prime (U.MVector s a) +newtype instance U.Vector (Prime a) = V_Prime (U.Vector a) + +instance U.Unbox a => U.Unbox (Prime a) + +instance M.MVector U.MVector a => M.MVector U.MVector (Prime a) where + {-# INLINE basicLength #-} + {-# INLINE basicUnsafeSlice #-} + {-# INLINE basicOverlaps #-} + {-# INLINE basicUnsafeNew #-} + {-# INLINE basicInitialize #-} + {-# INLINE basicUnsafeReplicate #-} + {-# INLINE basicUnsafeRead #-} + {-# INLINE basicUnsafeWrite #-} + {-# INLINE basicClear #-} + {-# INLINE basicSet #-} + {-# INLINE basicUnsafeCopy #-} + {-# INLINE basicUnsafeGrow #-} + basicLength (MV_Prime v) = M.basicLength v + basicUnsafeSlice i n (MV_Prime v) = MV_Prime $ M.basicUnsafeSlice i n v + basicOverlaps (MV_Prime v1) (MV_Prime v2) = M.basicOverlaps v1 v2 + basicUnsafeNew n = MV_Prime <$> M.basicUnsafeNew n + basicInitialize (MV_Prime v) = M.basicInitialize v + basicUnsafeReplicate n x = MV_Prime <$> M.basicUnsafeReplicate n (unPrime x) + basicUnsafeRead (MV_Prime v) i = Prime <$> M.basicUnsafeRead v i + basicUnsafeWrite (MV_Prime v) i x = M.basicUnsafeWrite v i (unPrime x) + basicClear (MV_Prime v) = M.basicClear v + basicSet (MV_Prime v) x = M.basicSet v (unPrime x) + basicUnsafeCopy (MV_Prime v1) (MV_Prime v2) = M.basicUnsafeCopy v1 v2 + basicUnsafeMove (MV_Prime v1) (MV_Prime v2) = M.basicUnsafeMove v1 v2 + basicUnsafeGrow (MV_Prime v) n = MV_Prime <$> M.basicUnsafeGrow v n + +instance G.Vector U.Vector a => G.Vector U.Vector (Prime a) where + {-# INLINE basicUnsafeFreeze #-} + {-# INLINE basicUnsafeThaw #-} + {-# INLINE basicLength #-} + {-# INLINE basicUnsafeSlice #-} + {-# INLINE basicUnsafeIndexM #-} + {-# INLINE elemseq #-} + basicUnsafeFreeze (MV_Prime v) = V_Prime <$> G.basicUnsafeFreeze v + basicUnsafeThaw (V_Prime v) = MV_Prime <$> G.basicUnsafeThaw v + basicLength (V_Prime v) = G.basicLength v + basicUnsafeSlice i n (V_Prime v) = V_Prime $ G.basicUnsafeSlice i n v + basicUnsafeIndexM (V_Prime v) i = Prime <$> G.basicUnsafeIndexM v i + basicUnsafeCopy (MV_Prime mv) (V_Prime v) = G.basicUnsafeCopy mv v + elemseq _ = seq + +-- | Convert between primes of different types, similar in spirit to 'toIntegralSized'. +-- +-- A simpler version of this function is: +-- +-- > toPrimeIntegral :: (Integral a, Integral b) => a -> Maybe b +-- > toPrimeIntegral (Prime a) +-- > | toInteger a == b = Just (Prime (fromInteger b)) +-- > | otherwise = Nothing +-- > where +-- > b = toInteger a +-- +-- The point of 'toPrimeIntegral' is to avoid redundant conversions and conditions, +-- when it is safe to do so, determining type sizes statically with 'bitSizeMaybe'. +-- For example, 'toPrimeIntegral' from 'Prime' 'Int' to 'Prime' 'Word' boils down to +-- 'Just' . 'fromIntegral'. +-- +toPrimeIntegral :: (Integral a, Integral b, Bits a, Bits b) => Prime a -> Maybe (Prime b) +toPrimeIntegral (Prime a) = case unsignedWidth b of + Nothing -> res + Just bW -> case unsignedWidth a of + Just aW + | aW <= bW -> res + _ + | a <= bit bW - 1 -> res + | otherwise -> Nothing + where + b = fromIntegral' a + res = Just (Prime b) +{-# INLINE toPrimeIntegral #-} + +unsignedWidth :: Bits a => a -> Maybe Int +unsignedWidth t + | isSigned t = subtract 1 <$> bitSizeMaybe t + | otherwise = bitSizeMaybe t +{-# INLINE unsignedWidth #-} diff --git a/Math/NumberTheory/Quadratic/EisensteinIntegers.hs b/Math/NumberTheory/Quadratic/EisensteinIntegers.hs --- a/Math/NumberTheory/Quadratic/EisensteinIntegers.hs +++ b/Math/NumberTheory/Quadratic/EisensteinIntegers.hs @@ -45,7 +45,7 @@ infix 6 :+ -- | An Eisenstein integer is @a + bω@, where @a@ and @b@ are both integers. -data EisensteinInteger = (:+) { real :: !Integer, imag :: !Integer } +data EisensteinInteger = !Integer :+ !Integer deriving (Eq, Ord, Generic) instance NFData EisensteinInteger @@ -80,7 +80,7 @@ times = (*) zero = 0 :+ 0 one = 1 :+ 0 - fromNatural n = fromIntegral n :+ 0 + fromNatural n = naturalToInteger n :+ 0 instance S.Ring EisensteinInteger where negate = negate @@ -88,14 +88,20 @@ -- | Returns an @EisensteinInteger@'s sign, and its associate in the first -- sextant. absSignum :: EisensteinInteger -> (EisensteinInteger, EisensteinInteger) +absSignum 0 = (0, 0) absSignum z@(a :+ b) - | a == 0 && b == 0 = (z, 0) -- origin - | a > b && b >= 0 = (z, 1) -- first sextant: 0 ≤ Arg(η) < π/3 - | b >= a && a > 0 = ((-ω) * z, 1 + ω) -- second sextant: π/3 ≤ Arg(η) < 2π/3 - | b > 0 && 0 >= a = ((-1 - ω) * z, ω) -- third sextant: 2π/3 ≤ Arg(η) < π - | a < b && b <= 0 = (- z, -1) -- fourth sextant: -π < Arg(η) < -2π/3 or Arg(η) = π - | b <= a && a < 0 = (ω * z, -1 - ω) -- fifth sextant: -2π/3 ≤ Arg(η) < -π/3 - | otherwise = ((1 + ω) * z, -ω) -- sixth sextant: -π/3 ≤ Arg(η) < 0 + -- first sextant: 0 ≤ Arg(z) < π/3 + | a > b && b >= 0 = (z, 1) + -- second sextant: π/3 ≤ Arg(z) < 2π/3 + | b >= a && a > 0 = (b :+ (b - a), 1 :+ 1) + -- third sextant: 2π/3 ≤ Arg(z) < π + | b > 0 && 0 >= a = ((b - a) :+ (-a), 0 :+ 1) + -- fourth sextant: -π ≤ Arg(z) < -2π/3 + | a < b && b <= 0 = (-z, -1) + -- fifth sextant: -2π/3 ≤ Arg(η) < -π/3 + | b <= a && a < 0 = ((-b) :+ (a - b), (-1) :+ (-1)) + -- sixth sextant: -π/3 ≤ Arg(η) < 0 + | otherwise = ((a - b) :+ a, 0 :+ (-1)) -- | List of all Eisenstein units, counterclockwise across all sextants, -- starting with @1@. @@ -109,24 +115,23 @@ instance GcdDomain EisensteinInteger instance Euclidean EisensteinInteger where - degree = fromInteger . norm - quotRem = divHelper - --- | Function that does most of the underlying work for @divMod@ and --- @quotRem@, apart from choosing the specific integer division algorithm. --- This is instead done by the calling function (either @divMod@ which uses --- @div@, or @quotRem@, which uses @quot@.) -divHelper - :: EisensteinInteger - -> EisensteinInteger - -> (EisensteinInteger, EisensteinInteger) -divHelper g h = (q, r) + degree = fromInteger . norm + quotRem x (d :+ 0) = quotRemInt x d + quotRem x y = (q, x - q * y) where - nr :+ ni = g * conjugate h - denom = norm h - q = ((nr + signum nr * denom `quot` 2) `quot` denom) :+ ((ni + signum ni * denom `quot` 2) `quot` denom) - r = g - h * q + (q, _) = quotRemInt (x * conjugate y) (norm y) +quotRemInt :: EisensteinInteger -> Integer -> (EisensteinInteger, EisensteinInteger) +quotRemInt z 1 = ( z, 0) +quotRemInt z (-1) = (-z, 0) +quotRemInt (a :+ b) c = (qa :+ qb, (ra - bumpA) :+ (rb - bumpB)) + where + halfC = abs c `quot` 2 + bumpA = signum a * halfC + bumpB = signum b * halfC + (qa, ra) = (a + bumpA) `quotRem` c + (qb, rb) = (b + bumpB) `quotRem` c + -- | Conjugate a Eisenstein integer. conjugate :: EisensteinInteger -> EisensteinInteger conjugate (a :+ b) = (a - b) :+ (-b) @@ -167,6 +172,7 @@ -- | Find an Eisenstein integer whose norm is the given prime number -- in the form @3k + 1@. -- +-- >>> import Math.NumberTheory.Primes (nextPrime) -- >>> findPrime (nextPrime 7) -- Prime 3+2*ω findPrime :: Prime Integer -> U.Prime EisensteinInteger @@ -263,7 +269,7 @@ gp' = Prime (x :+ (x - y)) (k, k', z') = divideByPrime gp gp' (unPrime p) e z - quotI (a :+ b) n = (a `quot` n :+ b `quot` n) + quotI (a :+ b) n = a `quot` n :+ b `quot` n -- | Remove @p@ and @conjugate p@ factors from the argument, where -- @p@ is an Eisenstein prime. diff --git a/Math/NumberTheory/Quadratic/GaussianIntegers.hs b/Math/NumberTheory/Quadratic/GaussianIntegers.hs --- a/Math/NumberTheory/Quadratic/GaussianIntegers.hs +++ b/Math/NumberTheory/Quadratic/GaussianIntegers.hs @@ -8,7 +8,6 @@ -- computing their prime factorisations. -- -{-# LANGUAGE BangPatterns #-} {-# LANGUAGE DeriveGeneric #-} {-# LANGUAGE TypeFamilies #-} @@ -74,36 +73,43 @@ times = (*) zero = 0 :+ 0 one = 1 :+ 0 - fromNatural n = fromIntegral n :+ 0 + fromNatural n = naturalToInteger n :+ 0 instance S.Ring GaussianInteger where negate = negate absSignum :: GaussianInteger -> (GaussianInteger, GaussianInteger) +absSignum 0 = (0, 0) absSignum z@(a :+ b) - | a == 0 && b == 0 = (z, 0) -- origin - | a > 0 && b >= 0 = (z, 1) -- first quadrant: (0, inf) x [0, inf)i - | a <= 0 && b > 0 = (b :+ (-a), ι) -- second quadrant: (-inf, 0] x (0, inf)i - | a < 0 && b <= 0 = ((-a) :+ (-b), -1) -- third quadrant: (-inf, 0) x (-inf, 0]i - | otherwise = ((-b) :+ a, -ι) -- fourth quadrant: [0, inf) x (-inf, 0)i + -- first quadrant: (0, inf) x [0, inf)i + | a > 0 && b >= 0 = (z, 1) + -- second quadrant: (-inf, 0] x (0, inf)i + | a <= 0 && b > 0 = (b :+ (-a), ι) + -- third quadrant: (-inf, 0) x (-inf, 0]i + | a < 0 && b <= 0 = (-z, -1) + -- fourth quadrant: [0, inf) x (-inf, 0)i + | otherwise = ((-b) :+ a, -ι) instance GcdDomain GaussianInteger instance Euclidean GaussianInteger where - degree = fromInteger . norm - quotRem = divHelper - -divHelper - :: GaussianInteger - -> GaussianInteger - -> (GaussianInteger, GaussianInteger) -divHelper g h = (q, r) + degree = fromInteger . norm + quotRem x (d :+ 0) = quotRemInt x d + quotRem x y = (q, x - q * y) where - nr :+ ni = g * conjugate h - denom = norm h - q = ((nr + signum nr * denom `quot` 2) `quot` denom) :+ ((ni + signum ni * denom `quot` 2) `quot` denom) - r = g - h * q + (q, _) = quotRemInt (x * conjugate y) (norm y) +quotRemInt :: GaussianInteger -> Integer -> (GaussianInteger, GaussianInteger) +quotRemInt z 1 = ( z, 0) +quotRemInt z (-1) = (-z, 0) +quotRemInt (a :+ b) c = (qa :+ qb, (ra - bumpA) :+ (rb - bumpB)) + where + halfC = abs c `quot` 2 + bumpA = signum a * halfC + bumpB = signum b * halfC + (qa, ra) = (a + bumpA) `quotRem` c + (qb, rb) = (b + bumpB) `quotRem` c + -- |Conjugate a Gaussian integer. conjugate :: GaussianInteger -> GaussianInteger conjugate (r :+ i) = r :+ (-i) @@ -138,6 +144,7 @@ -- of form 4k + 1 using -- . -- +-- >>> import Math.NumberTheory.Primes (nextPrime) -- >>> findPrime (nextPrime 5) -- Prime 2+ι findPrime :: Prime Integer -> U.Prime GaussianInteger @@ -205,8 +212,7 @@ where (d1, z') = go1 c 0 z d2 = c - d1 - z'' = head $ drop (wordToInt d2) - $ iterate (\g -> fromMaybe err $ (g * unPrime p) `quotEvenI` np) z' + z'' = iterate (\g -> fromMaybe err $ (g * unPrime p) `quotEvenI` np) z' !! wordToInt d2 go1 :: Word -> Word -> GaussianInteger -> (Word, GaussianInteger) go1 0 d z = (d, z) @@ -219,7 +225,7 @@ err = error $ "divideByPrime: malformed arguments" ++ show (p, np, k) quotI :: GaussianInteger -> Integer -> GaussianInteger -quotI (x :+ y) n = (x `quot` n :+ y `quot` n) +quotI (x :+ y) n = x `quot` n :+ y `quot` n quotEvenI :: GaussianInteger -> Integer -> Maybe GaussianInteger quotEvenI (x :+ y) n 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 @@ -12,26 +12,20 @@ -- Top-level definitions in this module are polymorphic, so the results of computations are not retained in memory. -- Make them monomorphic to take advantages of memoization. Compare -- --- >>> :set +s --- >>> binomial !! 1000 !! 1000 :: Integer +-- >>> binomial !! 1000 !! 1000 :: Integer -- (0.01 secs, 1,385,512 bytes) -- 1 --- (0.01 secs, 1,385,512 bytes) --- >>> binomial !! 1000 !! 1000 :: Integer +-- >>> binomial !! 1000 !! 1000 :: Integer -- (0.01 secs, 1,381,616 bytes) -- 1 --- (0.01 secs, 1,381,616 bytes) -- -- against -- -- >>> let binomial' = binomial :: [[Integer]] --- >>> binomial' !! 1000 !! 1000 :: Integer +-- >>> binomial' !! 1000 !! 1000 :: Integer -- (0.01 secs, 1,381,696 bytes) -- 1 --- (0.01 secs, 1,381,696 bytes) --- >>> binomial' !! 1000 !! 1000 :: Integer +-- >>> binomial' !! 1000 !! 1000 :: Integer -- (0.01 secs, 391,152 bytes) -- 1 --- (0.01 secs, 391,152 bytes) {-# LANGUAGE BangPatterns #-} -{-# LANGUAGE CPP #-} {-# LANGUAGE ScopedTypeVariables #-} module Math.NumberTheory.Recurrences.Bilinear @@ -65,7 +59,7 @@ -- | Infinite zero-based table of binomial coefficients (also known as Pascal triangle). -- --- prop> binomial !! n !! k == n! / k! / (n - k)! +-- > binomial !! n !! k == n! / k! / (n - k)! -- -- Note that 'binomial' !! n !! k is asymptotically slower -- than 'binomialLine' n !! k, @@ -82,7 +76,7 @@ -- | Pascal triangle, rotated by 45 degrees. -- --- prop> binomialRotated !! n !! k == (n + k)! / n! / k! == binomial !! (n + k) !! k +-- > binomialRotated !! n !! k == (n + k)! / n! / k! == binomial !! (n + k) !! k -- -- Note that 'binomialRotated' !! n !! k is asymptotically slower -- than 'binomialDiagonal' n !! k, @@ -119,7 +113,7 @@ -- [1,6,21,56,126,252] binomialDiagonal :: (Enum a, GcdDomain a) => a -> [a] binomialDiagonal n = scanl' - (\x k -> fromJust $ (x `times` (n `plus` k) `divide` k)) + (\x k -> fromJust (x `times` (n `plus` k) `divide` k)) one [one..] {-# SPECIALIZE binomialDiagonal :: Int -> [Int] #-} @@ -129,7 +123,7 @@ -- | Prime factors of a binomial coefficient. -- --- prop> binomialFactors n k == factorise (binomial !! n !! k) +-- > binomialFactors n k == factorise (binomial !! n !! k) -- -- >>> binomialFactors 10 4 -- [(Prime 2,1),(Prime 3,1),(Prime 5,1),(Prime 7,1)] @@ -140,7 +134,7 @@ | otherwise = filter ((/= 0) . snd) $ map (\p -> (p, mult (unPrime p) n - mult (unPrime p) (n - k) - mult (unPrime p) k)) - $ [minBound .. precPrime n] + [minBound .. precPrime n] where mult :: Word -> Word -> Word mult p m = go mp mp @@ -158,7 +152,7 @@ -- Complexity: @stirling1 !! n !! k@ is O(n ln n) bits long, its computation -- takes O(k n^2 ln n) time and forces thunks @stirling1 !! i !! j@ for @0 <= i <= n@ and @max(0, k - n + i) <= j <= k@. -- --- One could also consider 'Math.Combinat.Numbers.unsignedStirling1st' to compute stand-alone values. +-- One could also consider 'Math.Combinat.Numbers.unsignedStirling1st' from package to compute stand-alone values. stirling1 :: (Num a, Enum a) => [[a]] stirling1 = scanl f [1] [0..] where @@ -176,7 +170,7 @@ -- Complexity: @stirling2 !! n !! k@ is O(n ln n) bits long, its computation -- takes O(k n^2 ln n) time and forces thunks @stirling2 !! i !! j@ for @0 <= i <= n@ and @max(0, k - n + i) <= j <= k@. -- --- One could also consider 'Math.Combinat.Numbers.stirling2nd' to compute stand-alone values. +-- One could also consider 'Math.Combinat.Numbers.stirling2nd' from package to compute stand-alone values. stirling2 :: (Num a, Enum a) => [[a]] stirling2 = iterate f [1] where @@ -248,9 +242,9 @@ -- Complexity: @bernoulli !! n@ is O(n ln n) bits long, its computation -- takes O(n^3 ln n) time and forces thunks @stirling2 !! i !! j@ for @0 <= i <= n@ and @0 <= j <= i@. -- --- One could also consider 'Math.Combinat.Numbers.bernoulli' to compute stand-alone values. +-- One could also consider 'Math.Combinat.Numbers.bernoulli' from package to compute stand-alone values. bernoulli :: Integral a => [Ratio a] -bernoulli = helperForB_E_EP id (map recip [1..]) +bernoulli = helperForBEEP id (map recip [1..]) {-# SPECIALIZE bernoulli :: [Ratio Int] #-} {-# SPECIALIZE bernoulli :: [Rational] #-} @@ -265,7 +259,7 @@ faulhaberPoly p = zipWith (*) ((0:) $ reverse - $ take (p+1) $ bernoulli) + $ take (p + 1) bernoulli) $ map (% (fromIntegral p+1)) $ zipWith (*) (iterate negate (if odd p then 1 else -1)) $ binomial !! (p+1) @@ -278,7 +272,7 @@ -- >>> take 10 euler' :: [Rational] -- [1 % 1,0 % 1,(-1) % 1,0 % 1,5 % 1,0 % 1,(-61) % 1,0 % 1,1385 % 1,0 % 1] euler' :: forall a . Integral a => [Ratio a] -euler' = tail $ helperForB_E_EP tail as +euler' = tail $ helperForBEEP tail as where as :: [Ratio a] as = zipWith3 @@ -309,7 +303,7 @@ -- >>> take 10 eulerPolyAt1 :: [Rational] -- [1 % 1,1 % 2,0 % 1,(-1) % 4,0 % 1,1 % 2,0 % 1,(-17) % 8,0 % 1,31 % 2] eulerPolyAt1 :: forall a . Integral a => [Ratio a] -eulerPolyAt1 = tail $ helperForB_E_EP tail (map recip (iterate (2 *) 1)) +eulerPolyAt1 = tail $ helperForBEEP tail (map recip (iterate (2 *) 1)) {-# SPECIALIZE eulerPolyAt1 :: [Ratio Int] #-} {-# SPECIALIZE eulerPolyAt1 :: [Rational] #-} @@ -333,16 +327,16 @@ -- zipping four lists together with multiplication, with one of those lists -- being the sublists in @stirling2@, and two of them being the factorial -- sequence and @cycle [1, -1]@. The remaining list is passed to --- @helperForB_E_EP@ as an argument. +-- @helperForBEEP@ as an argument. -- -- Note: This function has a @([Ratio a] -> [Ratio a])@ argument because -- @bernoulli !! n@ will use, for all nonnegative @n@, every element in -- @stirling2 !! n@, while @euler, eulerPolyAt1@ only use -- @tail $ stirling2 !! n@. As such, this argument serves to pass @id@ -- in the former case, and @tail@ in the latter. -helperForB_E_EP :: Integral a => ([Ratio a] -> [Ratio a]) -> [Ratio a] -> [Ratio a] -helperForB_E_EP g xs = map (f . g) stirling2 +helperForBEEP :: Integral a => ([Ratio a] -> [Ratio a]) -> [Ratio a] -> [Ratio a] +helperForBEEP g xs = map (f . g) stirling2 where f = sum . zipWith4 (\sgn fact x stir -> sgn * fact * x * stir) (cycle [1, -1]) factorial xs -{-# SPECIALIZE helperForB_E_EP :: ([Ratio Int] -> [Ratio Int]) -> [Ratio Int] -> [Ratio Int] #-} -{-# SPECIALIZE helperForB_E_EP :: ([Rational] -> [Rational]) -> [Rational] -> [Rational] #-} +{-# SPECIALIZE helperForBEEP :: ([Ratio Int] -> [Ratio Int]) -> [Ratio Int] -> [Ratio Int] #-} +{-# SPECIALIZE helperForBEEP :: ([Rational] -> [Rational]) -> [Rational] -> [Rational] #-} diff --git a/Math/NumberTheory/Recurrences/Linear.hs b/Math/NumberTheory/Recurrences/Linear.hs --- a/Math/NumberTheory/Recurrences/Linear.hs +++ b/Math/NumberTheory/Recurrences/Linear.hs @@ -7,7 +7,6 @@ -- Efficient calculation of linear recurrent sequences, including Fibonacci and Lucas sequences. {-# LANGUAGE BangPatterns #-} -{-# LANGUAGE CPP #-} module Math.NumberTheory.Recurrences.Linear ( factorial @@ -19,8 +18,6 @@ , generalLucas ) where -#include "MachDeps.h" - import Data.Bits import Numeric.Natural import Math.NumberTheory.Primes @@ -41,7 +38,7 @@ -- | Prime factors of a factorial. -- --- prop> factorialFactors n == factorise (factorial !! n) +-- > factorialFactors n == factorise (factorial !! n) -- -- >>> factorialFactors 10 -- [(Prime 2,8),(Prime 3,4),(Prime 5,2),(Prime 7,1)] @@ -83,7 +80,7 @@ fibonacciPair n | n < 0 = let (f,g) = fibonacciPair (-(n+1)) in if testBit n 0 then (g, -f) else (-g, f) | n == 0 = (0, 1) - | otherwise = look (WORD_SIZE_IN_BITS - 2) + | otherwise = look (finiteBitSize (0 :: Word) - 2) where look k | testBit n k = go (k-1) 0 1 @@ -112,7 +109,7 @@ lucasPair n | n < 0 = let (f,g) = lucasPair (-(n+1)) in if testBit n 0 then (-g, f) else (g, -f) | n == 0 = (2, 1) - | otherwise = look (WORD_SIZE_IN_BITS - 2) + | otherwise = look (finiteBitSize (0 :: Word) - 2) where look k | testBit n k = go (k-1) 0 1 @@ -143,7 +140,7 @@ generalLucas p q k | k < 0 = error "generalLucas: negative index" | k == 0 = (0,1,2,p) - | otherwise = look (WORD_SIZE_IN_BITS - 2) + | otherwise = look (finiteBitSize (0 :: Word) - 2) where look i | testBit k i = go (i-1) 1 p p q diff --git a/Math/NumberTheory/Recurrences/Pentagonal.hs b/Math/NumberTheory/Recurrences/Pentagonal.hs --- a/Math/NumberTheory/Recurrences/Pentagonal.hs +++ b/Math/NumberTheory/Recurrences/Pentagonal.hs @@ -30,45 +30,17 @@ interleave (n : ns) (m : ms) = n : m : interleave ns ms interleave _ _ = [] --- | When calculating the @n@-th partition number @p(n)@ using the sum --- @p(n) = p(n-1) + p(n-2) - p(n-5) - p(n-7) + p(n-11) + ...@, the signs of each --- term alternate every two terms, starting with a positive sign. --- @pentagonalSigns@ takes a list of numbers and produces such an alternated --- sequence. --- Examples: --- --- >>> pentagonalSigns [1..5] --- [1,2,-3,-4,5] --- --- >>> pentagonalSigns [1..6] --- [1,2,-3,-4,5,6] -pentagonalSigns :: Num a => [a] -> [a] -pentagonalSigns = zipWith (*) (cycle [1, 1, -1, -1]) - --- [Implementation notes for partition function] --- --- @p(n) = p(n-1) + p(n-2) - p(n-5) - p(n-7) + p(n-11) + ...@, where @p(0) = 1@ +-- | @p(n) = p(n-1) + p(n-2) - p(n-5) - p(n-7) + p(n-11) + ...@, where @p(0) = 1@ -- and @p(k) = 0@ for a negative integer @k@. Uses a @Chimera@ from the -- @chimera@ package to memoize previous results. --- --- Example: calculating @partition !! 10@, assuming the memoization map is --- filled and called @dict@. --- --- * @tail [0, 1, 2, 5, 7, 12 ,15, 22, 26, 35, ..] == [1, 2, 5, 7, 12 ,15, 22, 26, 35, 40, ..]@. --- * @takeWhile (\m -> 10 - m >= 0) [1, 2, 5, 7, 12 ,15, 22, 26, 35, 40, ..] == [1, 2, 5, 7]@. --- * @map (\m -> dict ! fromIntegral (10 - m)) [1, 2, 5, 7] == [dict ! 9, dict ! 8, dict ! 5, dict ! 3] == [30, 22, 7, 3]@ --- * @pentagonalSigns [30, 22, 7, 3] == [30, 22, 7, 3] == [30, 22, -7, -3]@ --- * @sum [30, 22, -7, -3] == 42@ --- --- Notes: --- 1. @tail@ is applied to @pents@ because otherwise the calculation of --- @p(n)@ would involve a duplicated @p(n-1)@ term (see the above example). --- 2. Calculating @partition !! k@, where @k@ is any index equal or higher --- than @maxBound :: Int@ results in undefined behavior. - partitionF :: Num a => (Word -> a) -> Word -> a partitionF _ 0 = 1 -partitionF f n = sum $ pentagonalSigns $ map (f . (n -)) $ takeWhile (<= n) $ tail pents +partitionF f n + = sum + $ zipWith (*) (cycle [1, 1, -1, -1]) + $ map (f . (n -)) + $ takeWhile (<= n) + $ tail pents -- | Infinite zero-based table of . -- diff --git a/Math/NumberTheory/RootsOfUnity.hs b/Math/NumberTheory/RootsOfUnity.hs --- a/Math/NumberTheory/RootsOfUnity.hs +++ b/Math/NumberTheory/RootsOfUnity.hs @@ -50,7 +50,7 @@ -- | This Semigroup is in fact a group, so @'stimes'@ can be called with a negative first argument. instance Semigroup RootOfUnity where RootOfUnity q1 <> RootOfUnity q2 = toRootOfUnity (q1 + q2) - stimes k (RootOfUnity q) = toRootOfUnity (q * fromIntegral k) + stimes k (RootOfUnity q) = toRootOfUnity (q * (toInteger k % 1)) instance Monoid RootOfUnity where mappend = (<>) diff --git a/Math/NumberTheory/SmoothNumbers.hs b/Math/NumberTheory/SmoothNumbers.hs --- a/Math/NumberTheory/SmoothNumbers.hs +++ b/Math/NumberTheory/SmoothNumbers.hs @@ -12,7 +12,6 @@ {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE TypeApplications #-} module Math.NumberTheory.SmoothNumbers ( SmoothBasis @@ -23,7 +22,7 @@ , smoothOver' ) where -import Prelude hiding (div, mod, gcd) +import Prelude hiding (div, mod, gcd, (+)) import Data.Euclidean import Data.List (nub) import Data.Maybe diff --git a/Math/NumberTheory/Utils.hs b/Math/NumberTheory/Utils.hs --- a/Math/NumberTheory/Utils.hs +++ b/Math/NumberTheory/Utils.hs @@ -6,9 +6,18 @@ -- -- Some utilities, mostly for bit twiddling. -- -{-# LANGUAGE CPP, MagicHash, UnboxedTuples, BangPatterns #-} + +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE MagicHash #-} +{-# LANGUAGE UnboxedTuples #-} +{-# LANGUAGE RankNTypes #-} +{-# LANGUAGE KindSignatures #-} +{-# LANGUAGE DataKinds #-} +{-# LANGUAGE GADTs #-} + module Math.NumberTheory.Utils - ( shiftToOddCount + ( SomeKnown(..) + , shiftToOddCount , shiftToOdd , shiftToOdd# , shiftToOddCount# @@ -22,10 +31,10 @@ , toWheel30 , fromWheel30 + , withSomeKnown + , intVal ) where -#include "MachDeps.h" - import Prelude hiding (mod, quotRem) import qualified Prelude as P @@ -34,7 +43,10 @@ import Data.Semiring (Semiring(..), isZero) import GHC.Base import GHC.Integer.GMP.Internals +import qualified Math.NumberTheory.Utils.FromIntegral as UT import GHC.Natural +import GHC.TypeNats +import Math.NumberTheory.Utils.FromIntegral (intToWord) -- | Remove factors of @2@ and count them. If -- @n = 2^k*m@ with @m@ odd, the result is @(k, m)@. @@ -47,7 +59,7 @@ #-} {-# INLINE [1] shiftToOddCount #-} shiftToOddCount :: Integral a => a -> (Word, a) -shiftToOddCount n = case shiftOCInteger (fromIntegral n) of +shiftToOddCount n = case shiftOCInteger (toInteger n) of (z, o) -> (z, fromInteger o) -- | Specialised version for @'Word'@. @@ -71,10 +83,10 @@ (# z#, w# #) -> (W# z#, wordToInteger w#) shiftOCInteger n@(Jp# bn#) = case bigNatZeroCount bn# of 0## -> (0, n) - z# -> (W# z#, bigNatToInteger (bn# `shiftRBigNat` (word2Int# z#))) + z# -> (W# z#, bigNatToInteger (bn# `shiftRBigNat` word2Int# z#)) shiftOCInteger n@(Jn# bn#) = case bigNatZeroCount bn# of 0## -> (0, n) - z# -> (W# z#, bigNatToNegInteger (bn# `shiftRBigNat` (word2Int# z#))) + z# -> (W# z#, bigNatToNegInteger (bn# `shiftRBigNat` word2Int# z#)) -- | Specialised version for @'Natural'@. -- Precondition: argument nonzero (not checked). @@ -85,21 +97,22 @@ (# z#, w# #) -> (W# z#, NatS# w#) shiftOCNatural n@(NatJ# bn#) = case bigNatZeroCount bn# of 0## -> (0, n) - z# -> (W# z#, bigNatToNatural (bn# `shiftRBigNat` (word2Int# z#))) + z# -> (W# z#, bigNatToNatural (bn# `shiftRBigNat` word2Int# z#)) shiftToOddCountBigNat :: BigNat -> (Word, BigNat) shiftToOddCountBigNat bn# = case bigNatZeroCount bn# of 0## -> (0, bn#) - z# -> (W# z#, bn# `shiftRBigNat` (word2Int# z#)) + z# -> (W# z#, bn# `shiftRBigNat` word2Int# z#) -- | Count trailing zeros in a @'BigNat'@. -- Precondition: argument nonzero (not checked, Integer invariant). bigNatZeroCount :: BigNat -> Word# bigNatZeroCount bn# = count 0## 0# where + !(W# bitSize#) = intToWord (finiteBitSize (0 :: Word)) count a# i# = case indexBigNat# bn# i# of - 0## -> count (a# `plusWord#` WORD_SIZE_IN_BITS##) (i# +# 1#) + 0## -> count (a# `plusWord#` bitSize#) (i# +# 1#) w# -> a# `plusWord#` ctz# w# -- | Remove factors of @2@. If @n = 2^k*m@ with @m@ odd, the result is @m@. @@ -111,7 +124,7 @@ #-} {-# INLINE [1] shiftToOdd #-} shiftToOdd :: Integral a => a -> a -shiftToOdd n = fromInteger (shiftOInteger (fromIntegral n)) +shiftToOdd n = fromInteger (shiftOInteger (toInteger n)) -- | Specialised version for @'Int'@. -- Precondition: argument nonzero (not checked). @@ -129,10 +142,10 @@ shiftOInteger (S# i#) = wordToInteger (shiftToOdd# (int2Word# i#)) shiftOInteger n@(Jp# bn#) = case bigNatZeroCount bn# of 0## -> n - z# -> bigNatToInteger (bn# `shiftRBigNat` (word2Int# z#)) + z# -> bigNatToInteger (bn# `shiftRBigNat` word2Int# z#) shiftOInteger n@(Jn# bn#) = case bigNatZeroCount bn# of 0## -> n - z# -> bigNatToNegInteger (bn# `shiftRBigNat` (word2Int# z#)) + z# -> bigNatToNegInteger (bn# `shiftRBigNat` word2Int# z#) -- | Shift argument right until the result is odd. -- Precondition: argument not @0@, not checked. @@ -200,3 +213,15 @@ fromWheel30 :: (Num a, Bits a) => a -> a fromWheel30 i = ((i `shiftL` 2 - i `shiftR` 2) .|. 1) + ((i `shiftL` 1 - i `shiftR` 1) .&. 2) + +------------------------------------------------------------------------------- +-- Helpers for dealing with data types parametrised by natural numbers. + +data SomeKnown (f :: Nat -> Type) where + SomeKnown :: KnownNat k => f k -> SomeKnown f + +withSomeKnown :: (forall k. KnownNat k => f k -> a) -> SomeKnown f -> a +withSomeKnown f (SomeKnown x) = f x + +intVal :: KnownNat k => a k -> Int +intVal = UT.naturalToInt . natVal 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 @@ -9,7 +9,6 @@ {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeApplications #-} -{-# LANGUAGE ViewPatterns #-} module Math.NumberTheory.Utils.DirichletSeries ( DirichletSeries diff --git a/Math/NumberTheory/Utils/FromIntegral.hs b/Math/NumberTheory/Utils/FromIntegral.hs --- a/Math/NumberTheory/Utils/FromIntegral.hs +++ b/Math/NumberTheory/Utils/FromIntegral.hs @@ -7,19 +7,39 @@ -- Monomorphic `fromIntegral`. -- -{-# LANGUAGE CPP #-} +{-# LANGUAGE CPP #-} module Math.NumberTheory.Utils.FromIntegral ( wordToInt , wordToInteger , intToWord + , intToInt8 + , intToInt64 + , int8ToInt64 + , intToWord8 + , intToWord64 + , int8ToInt + , int64ToInt + , word8ToInt + , word64ToInt , intToInteger + , int16ToInteger + , int64ToInteger + , word64ToInteger , naturalToInteger , integerToNatural , integerToWord + , integerToWord64 , integerToInt + , integerToInt64 + , intToNatural + , naturalToInt + , intToDouble + , fromIntegral' ) where +import Data.Int +import Data.Word import Numeric.Natural wordToInt :: Word -> Int @@ -34,22 +54,99 @@ intToWord = fromIntegral {-# INLINE intToWord #-} +intToInt8 :: Int -> Int8 +intToInt8 = fromIntegral +{-# INLINE intToInt8 #-} + +intToInt64 :: Int -> Int64 +intToInt64 = fromIntegral +{-# INLINE intToInt64 #-} + +int8ToInt64 :: Int8 -> Int64 +int8ToInt64 = fromIntegral +{-# INLINE int8ToInt64 #-} + +intToWord8 :: Int -> Word8 +intToWord8 = fromIntegral +{-# INLINE intToWord8 #-} + +intToWord64 :: Int -> Word64 +intToWord64 = fromIntegral +{-# INLINE intToWord64 #-} + +int8ToInt :: Int8 -> Int +int8ToInt = fromIntegral +{-# INLINE int8ToInt #-} + +int64ToInt :: Int64 -> Int +int64ToInt = fromIntegral +{-# INLINE int64ToInt #-} + +word8ToInt :: Word8 -> Int +word8ToInt = fromIntegral +{-# INLINE word8ToInt #-} + +word64ToInt :: Word64 -> Int +word64ToInt = fromIntegral +{-# INLINE word64ToInt #-} + intToInteger :: Int -> Integer intToInteger = fromIntegral {-# INLINE intToInteger #-} +int16ToInteger :: Int16 -> Integer +int16ToInteger = fromIntegral +{-# INLINE int16ToInteger #-} + +int64ToInteger :: Int64 -> Integer +int64ToInteger = fromIntegral +{-# INLINE int64ToInteger #-} + +word64ToInteger :: Word64 -> Integer +word64ToInteger = fromIntegral +{-# INLINE word64ToInteger #-} + naturalToInteger :: Natural -> Integer naturalToInteger = fromIntegral {-# INLINE naturalToInteger #-} integerToNatural :: Integer -> Natural -integerToNatural = fromIntegral +integerToNatural = fromIntegral' {-# INLINE integerToNatural #-} integerToWord :: Integer -> Word integerToWord = fromIntegral {-# INLINE integerToWord #-} +integerToWord64 :: Integer -> Word64 +integerToWord64 = fromIntegral +{-# INLINE integerToWord64 #-} + integerToInt :: Integer -> Int integerToInt = fromIntegral {-# INLINE integerToInt #-} + +integerToInt64 :: Integer -> Int64 +integerToInt64 = fromIntegral +{-# INLINE integerToInt64 #-} + +intToNatural :: Int -> Natural +intToNatural = fromIntegral +{-# INLINE intToNatural #-} + +naturalToInt :: Natural -> Int +naturalToInt = fromIntegral +{-# INLINE naturalToInt #-} + +intToDouble :: Int -> Double +intToDouble = fromIntegral +{-# INLINE intToDouble #-} + +fromIntegral' :: (Integral a, Num b) => a -> b +#if __GLASGOW_HASKELL__ == 900 && __GLASGOW_HASKELL_PATCHLEVEL1__ == 1 +-- Cannot use fromIntegral because of https://gitlab.haskell.org/ghc/ghc/-/issues/19411 +fromIntegral' = fromInteger . toInteger +#else +fromIntegral' = fromIntegral +#endif +{-# INLINE fromIntegral' #-} diff --git a/Math/NumberTheory/Utils/Hyperbola.hs b/Math/NumberTheory/Utils/Hyperbola.hs --- a/Math/NumberTheory/Utils/Hyperbola.hs +++ b/Math/NumberTheory/Utils/Hyperbola.hs @@ -7,9 +7,6 @@ -- Highest points under hyperbola. -- -{-# LANGUAGE CPP #-} -{-# LANGUAGE LambdaCase #-} - module Math.NumberTheory.Utils.Hyperbola ( pointsUnderHyperbola ) where @@ -50,19 +47,15 @@ -- | bresenham(x+1) -> bresenham(x) for x >= (2n)^1/3 stepBack :: Bresenham -> Bresenham -stepBack (Bresenham x' beta' gamma' delta1' epsilon') = - if eps >= x - then (if eps >= x `shiftL` 1 - then {- delta2 = 2 -} - let delta1 = delta1' + 2 in (Bresenham x (beta' + delta1) (gamma' + delta1 `shiftL` 1 - x `shiftL` 1) delta1 (eps - x `shiftL` 1)) - else {- delta1 = 1 -} - let delta1 = delta1' + 1 in (Bresenham x (beta' + delta1) (gamma' + delta1 `shiftL` 1 - x) delta1 (eps - x)) - ) - else (if eps >= 0 - then {- delta2 = 0 -} - (Bresenham x (beta' + delta1') (gamma' + delta1' `shiftL` 1) delta1' eps) - else {- delta2 = -1 -} - let delta1 = delta1' - 1 in (Bresenham x (beta' + delta1) (gamma' + delta1 `shiftL` 1 + x) delta1 (eps + x))) +stepBack (Bresenham x' beta' gamma' delta1' epsilon') + | eps >= x `shiftL` 1 {- delta2 = 2 -} + = let delta1 = delta1' + 2 in Bresenham x (beta' + delta1) (gamma' + delta1 `shiftL` 1 - x `shiftL` 1) delta1 (eps - x `shiftL` 1) + | eps >= x {- delta1 = 1 -} + = let delta1 = delta1' + 1 in Bresenham x (beta' + delta1) (gamma' + delta1 `shiftL` 1 - x) delta1 (eps - x) + | eps >= 0 {- delta2 = 0 -} + = Bresenham x (beta' + delta1') (gamma' + delta1' `shiftL` 1) delta1' eps + | otherwise {- delta2 = -1 -} + = let delta1 = delta1' - 1 in Bresenham x (beta' + delta1) (gamma' + delta1 `shiftL` 1 + x) delta1 (eps + x) where x = x' - 1 eps = epsilon' + gamma' diff --git a/Math/NumberTheory/Zeta/Dirichlet.hs b/Math/NumberTheory/Zeta/Dirichlet.hs --- a/Math/NumberTheory/Zeta/Dirichlet.hs +++ b/Math/NumberTheory/Zeta/Dirichlet.hs @@ -24,6 +24,7 @@ -- | Infinite sequence of exact values of Dirichlet beta-function at odd arguments, starting with @β(1)@. -- +-- >>> import Data.ExactPi -- >>> approximateValue (betasOdd !! 25) :: Double -- 0.9999999999999987 -- >>> import Data.Number.Fixed diff --git a/Math/NumberTheory/Zeta/Hurwitz.hs b/Math/NumberTheory/Zeta/Hurwitz.hs --- a/Math/NumberTheory/Zeta/Hurwitz.hs +++ b/Math/NumberTheory/Zeta/Hurwitz.hs @@ -50,7 +50,7 @@ -- @a + n@ aPlusN :: a - aPlusN = a + fromIntegral digitsOfPrecision + aPlusN = a + fromInteger digitsOfPrecision -- @[(a + n)^s | s <- [0, 1, 2 ..]]@ powsOfAPlusN :: [a] @@ -75,7 +75,7 @@ (\powOfA int -> powOfA * fromInteger int) powsOfAPlusN [-1, 0..] - in map ((/) aPlusN) denoms + in map (aPlusN /) denoms -- [ 1 | ] -- [ ----------- | s <- [0 ..] ] @@ -110,7 +110,7 @@ fracs :: [a] fracs = map - (\pochh -> sum $ zipWith (\s p -> s * fromInteger p) second pochh) + (sum . zipWith (\s p -> s * fromInteger p) second) pochhammers -- Infinite list of @T@ values in 4.8.5 formula, for every @s@ in diff --git a/Math/NumberTheory/Zeta/Utils.hs b/Math/NumberTheory/Zeta/Utils.hs --- a/Math/NumberTheory/Zeta/Utils.hs +++ b/Math/NumberTheory/Zeta/Utils.hs @@ -15,7 +15,7 @@ -- | Joins two lists element-by-element together into one, starting with the -- first one provided as argument. -- --- >>> take 10 $ intertwine [0, 2 ..] [1, 3 ..] +-- >>> take 10 (intertwine [0, 2 ..] [1, 3 ..]) -- [0,1,2,3,4,5,6,7,8,9] intertwine :: [a] -> [a] -> [a] intertwine [] ys = ys diff --git a/arithmoi.cabal b/arithmoi.cabal --- a/arithmoi.cabal +++ b/arithmoi.cabal @@ -1,12 +1,11 @@ name: arithmoi -version: 0.11.0.1 -cabal-version: >=1.10 +version: 0.12.0.0 +cabal-version: 2.0 build-type: Simple license: MIT license-file: LICENSE -copyright: (c) 2016-2020 Andrew Lelechenko, 2016-2019 Carter Schonwald, 2011 Daniel Fischer -maintainer: Andrew Lelechenko andrew dot lelechenko at gmail dot com, - Carter Schonwald carter at wellposed dot com +copyright: (c) 2016-2021 Andrew Lelechenko, 2016-2019 Carter Schonwald, 2011 Daniel Fischer +maintainer: Andrew Lelechenko homepage: https://github.com/Bodigrim/arithmoi bug-reports: https://github.com/Bodigrim/arithmoi/issues synopsis: Efficient basic number-theoretic functions. @@ -18,7 +17,7 @@ powers (integer roots and tests, modular exponentiation). category: Math, Algorithms, Number Theory author: Andrew Lelechenko, Daniel Fischer -tested-with: GHC ==8.2.2 GHC ==8.4.4 GHC ==8.6.5 GHC ==8.8.3 GHC ==8.10.1 +tested-with: GHC ==8.2.2 GHC ==8.4.4 GHC ==8.6.5 GHC ==8.8.4 GHC ==8.10.4 GHC ==9.0.1 extra-source-files: changelog.md @@ -35,7 +34,7 @@ constraints, deepseq, exact-pi >=0.5, - integer-gmp <1.1, + integer-gmp <1.2, integer-logarithms >=1.0, integer-roots >=1.0, mod, @@ -51,30 +50,23 @@ Math.NumberTheory.ArithmeticFunctions.Moebius Math.NumberTheory.ArithmeticFunctions.SieveBlock Math.NumberTheory.Curves.Montgomery + Math.NumberTheory.Diophantine Math.NumberTheory.DirichletCharacters - Math.NumberTheory.Euclidean Math.NumberTheory.Euclidean.Coprimes Math.NumberTheory.Moduli Math.NumberTheory.Moduli.Chinese Math.NumberTheory.Moduli.Class - Math.NumberTheory.Moduli.DiscreteLogarithm + Math.NumberTheory.Moduli.Cbrt Math.NumberTheory.Moduli.Equations - Math.NumberTheory.Moduli.Jacobi Math.NumberTheory.Moduli.Multiplicative - Math.NumberTheory.Moduli.PrimitiveRoot Math.NumberTheory.Moduli.Singleton Math.NumberTheory.Moduli.Sqrt Math.NumberTheory.MoebiusInversion - Math.NumberTheory.Powers - Math.NumberTheory.Powers.Cubes - Math.NumberTheory.Powers.Fourth - Math.NumberTheory.Powers.General Math.NumberTheory.Powers.Modular - Math.NumberTheory.Powers.Squares - Math.NumberTheory.Powers.Squares.Internal Math.NumberTheory.Prefactored Math.NumberTheory.Primes Math.NumberTheory.Primes.Counting + Math.NumberTheory.Primes.IntSet Math.NumberTheory.Primes.Testing Math.NumberTheory.Quadratic.GaussianIntegers Math.NumberTheory.Quadratic.EisensteinIntegers @@ -110,7 +102,7 @@ Math.NumberTheory.Zeta.Riemann Math.NumberTheory.Zeta.Utils default-language: Haskell2010 - ghc-options: -Wall -Widentities -Wcompat + ghc-options: -Wall -Widentities -Wcompat -Wno-deprecations test-suite arithmoi-tests build-depends: @@ -118,11 +110,12 @@ arithmoi, containers, exact-pi >=0.4.1.1, - integer-gmp <1.1, + integer-gmp <1.2, integer-roots >=1.0, mod, QuickCheck >=2.10, quickcheck-classes >=0.6.3, + random >=1.0 && <1.3, semirings >=0.2, smallcheck >=1.2 && <1.3, tasty >=0.10, @@ -131,13 +124,15 @@ tasty-rerun >=1.1.17, tasty-smallcheck >=0.8 && <0.9, transformers >=0.5, - vector + vector, + vector-sized other-modules: Math.NumberTheory.ArithmeticFunctionsTests Math.NumberTheory.ArithmeticFunctions.InverseTests Math.NumberTheory.ArithmeticFunctions.MertensTests Math.NumberTheory.ArithmeticFunctions.SieveBlockTests Math.NumberTheory.CurvesTests + Math.NumberTheory.DiophantineTests Math.NumberTheory.DirichletCharactersTests Math.NumberTheory.EisensteinIntegersTests Math.NumberTheory.GaussianIntegersTests @@ -145,16 +140,18 @@ Math.NumberTheory.Moduli.ChineseTests Math.NumberTheory.Moduli.DiscreteLogarithmTests Math.NumberTheory.Moduli.ClassTests + Math.NumberTheory.Moduli.CbrtTests Math.NumberTheory.Moduli.EquationsTests Math.NumberTheory.Moduli.JacobiTests Math.NumberTheory.Moduli.PrimitiveRootTests Math.NumberTheory.Moduli.SingletonTests Math.NumberTheory.Moduli.SqrtTests Math.NumberTheory.MoebiusInversionTests - Math.NumberTheory.Powers.ModularTests Math.NumberTheory.PrefactoredTests Math.NumberTheory.Primes.CountingTests Math.NumberTheory.Primes.FactorisationTests + -- Math.NumberTheory.Primes.LinearAlgebraTests + -- Math.NumberTheory.Primes.QuadraticSieveTests Math.NumberTheory.Primes.SequenceTests Math.NumberTheory.Primes.SieveTests Math.NumberTheory.Primes.TestingTests @@ -174,9 +171,9 @@ main-is: Test.hs default-language: Haskell2010 hs-source-dirs: test-suite - ghc-options: -Wall -Widentities -Wcompat + ghc-options: -Wall -Widentities -Wcompat -threaded -benchmark arithmoi-gauge +benchmark arithmoi-bench build-depends: base, arithmoi, @@ -184,12 +181,15 @@ constraints, containers, deepseq, - gauge, integer-logarithms, mod, random, semirings, vector + build-depends: + tasty-bench + mixins: + tasty-bench (Test.Tasty.Bench as Gauge.Main) other-modules: Math.NumberTheory.ArithmeticFunctionsBench Math.NumberTheory.DiscreteLogarithmBench diff --git a/benchmark/Math/NumberTheory/MertensBench.hs b/benchmark/Math/NumberTheory/MertensBench.hs --- a/benchmark/Math/NumberTheory/MertensBench.hs +++ b/benchmark/Math/NumberTheory/MertensBench.hs @@ -1,6 +1,3 @@ -{-# LANGUAGE CPP #-} -{-# LANGUAGE LambdaCase #-} - {-# OPTIONS_GHC -fno-warn-type-defaults #-} module Math.NumberTheory.MertensBench diff --git a/benchmark/Math/NumberTheory/PrimitiveRootsBench.hs b/benchmark/Math/NumberTheory/PrimitiveRootsBench.hs --- a/benchmark/Math/NumberTheory/PrimitiveRootsBench.hs +++ b/benchmark/Math/NumberTheory/PrimitiveRootsBench.hs @@ -39,7 +39,7 @@ , bench "10000000000000061" $ nf cyclicWrap (10^16 + 61) -- large prime , bench "2*3^20000" $ nf cyclicWrap (2*3^20000) -- twice prime to large power , bench "10000000000000046" $ nf cyclicWrap (10^16 + 46) -- twice large prime - , bench "224403121196654400" $ nf cyclicWrap (224403121196654400) -- highly composite + , bench "224403121196654400" $ nf cyclicWrap 224403121196654400 -- highly composite ] , bgroup "check prim roots" [ bench "3^20000" $ nf (primRootWrap 3 20000) 2 -- prime to large power diff --git a/benchmark/Math/NumberTheory/RecurrencesBench.hs b/benchmark/Math/NumberTheory/RecurrencesBench.hs --- a/benchmark/Math/NumberTheory/RecurrencesBench.hs +++ b/benchmark/Math/NumberTheory/RecurrencesBench.hs @@ -32,7 +32,7 @@ benchSuite :: Benchmark benchSuite = bgroup "Recurrences" [ bgroup "Bilinear" - [ benchTriangle "binomial" binomial 1000 + [ benchTriangle "binomial" binomial 100 , benchTriangle "stirling1" stirling1 100 , benchTriangle "stirling2" stirling2 100 , benchTriangle "eulerian1" eulerian1 100 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 @@ -20,16 +20,14 @@ filterIsPrimeBench :: Benchmark filterIsPrimeBench = bgroup "filterIsPrime" $ - map (\(x, y) -> bench (show (x, y)) $ nf filterIsPrime (x, x + y)) - [ (10 ^ x, 10 ^ y) + [ bench (show (10^x, 10^y)) $ nf filterIsPrime (10^x, 10^x + 10^y) | x <- [5..8] , y <- [3..x-1] ] eratosthenesBench :: Benchmark eratosthenesBench = bgroup "eratosthenes" $ - map (\(x, y) -> bench (show (x, y)) $ nf eratosthenes (x, x + y)) - [ (10 ^ x, 10 ^ y) + [ bench (show (10^x, 10^y)) $ nf eratosthenes (10^x, 10^x + 10^y) | x <- [10..17] , y <- [6..x-1] , x == 10 || y == 7 diff --git a/benchmark/Math/NumberTheory/SieveBlockBench.hs b/benchmark/Math/NumberTheory/SieveBlockBench.hs --- a/benchmark/Math/NumberTheory/SieveBlockBench.hs +++ b/benchmark/Math/NumberTheory/SieveBlockBench.hs @@ -45,7 +45,7 @@ carmichaelBlockConfig = SieveBlockConfig { sbcEmpty = 1 -- There is a specialized 'gcd' for Word, but not 'lcm'. - , sbcAppend = (\x y -> (x `quot` (gcd x y)) * y) + , sbcAppend = \x y -> (x `quot` gcd x y) * y , sbcFunctionOnPrimePower = carmichaelHelper . unPrime } diff --git a/benchmark/Math/NumberTheory/ZetaBench.hs b/benchmark/Math/NumberTheory/ZetaBench.hs --- a/benchmark/Math/NumberTheory/ZetaBench.hs +++ b/benchmark/Math/NumberTheory/ZetaBench.hs @@ -10,6 +10,6 @@ benchSuite :: Benchmark benchSuite = bgroup "Zeta" - [ bench "riemann zeta" $ nf (\eps -> sum $ take 20 $ zetas eps) (1e-15 :: Double) - , bench "dirichlet beta" $ nf (\eps -> sum $ take 20 $ betas eps) (1e-15 :: Double) + [ bench "riemann zeta" $ nf (sum . take 20 . zetas) (1e-15 :: Double) + , bench "dirichlet beta" $ nf (sum . take 20 . betas) (1e-15 :: Double) ] diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -1,5 +1,22 @@ # Changelog +## 0.12.0.0 + +### Added + +* Define cubic symbol ([#194](https://github.com/Bodigrim/arithmoi/pull/194)). +* Add `instance Unbox (Prime a)` and `toPrimeIntegral` helper ([#201](https://github.com/Bodigrim/arithmoi/pull/201)). +* Implement Cornacchia's algorithm for Diophantine equations ([#195](https://github.com/Bodigrim/arithmoi/pull/195)). +* Define a wrapper `PrimeIntSet` for sets of primes ([#205](https://github.com/Bodigrim/arithmoi/pull/205)). + +### Deprecated + +* Deprecate `Math.NumberTheory.Powers.Modular`, use `Data.Mod` or `Data.Mod.Word` instead. + +### Removed + +* Remove modules and functions, deprecated in the previous release. + ## 0.11.0.1 ### Changed @@ -44,6 +61,10 @@ * Deprecate `Math.NumberTheory.Moduli.{DiscreteLogarithm,PrimitiveRoot}`, use `Math.NumberTheory.Moduli.Multiplicative` instead. + +### Removed + +* Remove modules and functions, deprecated in the previous release. ### Fixed diff --git a/test-suite/Math/NumberTheory/ArithmeticFunctionsTests.hs b/test-suite/Math/NumberTheory/ArithmeticFunctionsTests.hs --- a/test-suite/Math/NumberTheory/ArithmeticFunctionsTests.hs +++ b/test-suite/Math/NumberTheory/ArithmeticFunctionsTests.hs @@ -292,7 +292,7 @@ nFreedomProperty3 (Power n) (Positive m) = let n' | n == maxBound = n | otherwise = n + 1 - zet = 1 / zetas 1e-14 !! (fromIntegral n') :: Double + zet = 1 / zetas 1e-14 !! fromIntegral n' :: Double m' = 100 * m nfree = fromIntegral m' / fromIntegral (head (drop (m' - 1) $ nFrees n' :: [Integer])) diff --git a/test-suite/Math/NumberTheory/CurvesTests.hs b/test-suite/Math/NumberTheory/CurvesTests.hs --- a/test-suite/Math/NumberTheory/CurvesTests.hs +++ b/test-suite/Math/NumberTheory/CurvesTests.hs @@ -22,9 +22,7 @@ import Math.NumberTheory.TestUtils (==>?) :: Maybe a -> (a -> Property) -> Property -x ==>? f = case x of - Nothing -> discard - Just y -> f y +x ==>? f = maybe discard f x isValid :: KnownNat n => Point a24 n -> Property isValid p diff --git a/test-suite/Math/NumberTheory/DiophantineTests.hs b/test-suite/Math/NumberTheory/DiophantineTests.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/NumberTheory/DiophantineTests.hs @@ -0,0 +1,40 @@ +-- Tests for Math.NumberTheory.Diophantine + +{-# LANGUAGE CPP #-} + +{-# OPTIONS_GHC -fno-warn-type-defaults #-} + +module Math.NumberTheory.DiophantineTests + ( testSuite + ) where + +import Data.List (sort) + +import Test.Tasty + +import Math.NumberTheory.Diophantine +import Math.NumberTheory.Roots (integerSquareRoot) +import Math.NumberTheory.TestUtils + +cornacchiaTest :: Positive Integer -> Positive Integer -> Bool +cornacchiaTest (Positive d) (Positive a) = gcd d m /= 1 || all checkSoln (cornacchia d m) + where m = d + a + checkSoln (x, y) = x*x + d*y*y == m + +-- Testing against a slower reference implementation on coprime inputs +cornacchiaBruteForce :: Positive Integer -> Positive Integer -> Bool +cornacchiaBruteForce (Positive d) (Positive a) = gcd d m /= 1 || findSolutions [] 1 == sort (cornacchia d m) + where m = d + a + -- Simple O(sqrt (m/d)) brute force by considering all possible y values + findSolutions acc y + | x2 <= 0 = acc + | x*x == x2 = findSolutions ((x,y) : acc) (y+1) + | otherwise = findSolutions acc (y+1) + where x2 = m - d*y*y + x = integerSquareRoot x2 + +testSuite :: TestTree +testSuite = testGroup "Diophantine" + [ testSmallAndQuick "Cornacchia correct" cornacchiaTest + , testSmallAndQuick "Cornacchia same solutions as brute force" cornacchiaBruteForce + ] 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 @@ -8,7 +8,6 @@ -- {-# LANGUAGE GADTs #-} -{-# LANGUAGE PatternSynonyms #-} {-# LANGUAGE Rank2Types #-} {-# LANGUAGE ScopedTypeVariables #-} {-# LANGUAGE TypeApplications #-} @@ -196,13 +195,8 @@ succValid = validChar . succ inducedValid :: forall d. KnownNat d => DirichletCharacter d -> Positive Natural -> Bool -inducedValid chi (Positive k) = - case someNatVal (d*k) of - SomeNat (Proxy :: Proxy n) -> - case induced @n chi of - Just chi2 -> validChar chi2 - Nothing -> False - where d = natVal @d Proxy +inducedValid chi (Positive k) = case someNatVal (natVal @d Proxy * k) of + SomeNat (Proxy :: Proxy n) -> maybe False validChar (induced @n chi) jacobiValid :: Positive Natural -> Bool jacobiValid (Positive n) = diff --git a/test-suite/Math/NumberTheory/EisensteinIntegersTests.hs b/test-suite/Math/NumberTheory/EisensteinIntegersTests.hs --- a/test-suite/Math/NumberTheory/EisensteinIntegersTests.hs +++ b/test-suite/Math/NumberTheory/EisensteinIntegersTests.hs @@ -48,7 +48,7 @@ -- | Verify that @rem@ produces a remainder smaller than the divisor with -- regards to the Euclidean domain's function. remProperty1 :: E.EisensteinInteger -> E.EisensteinInteger -> Bool -remProperty1 x y = (y == 0) || (E.norm $ x `rem` y) < (E.norm y) +remProperty1 x y = (y == 0) || E.norm (x `rem` y) < E.norm y -- | Verify that @quot@ and @rem@ are what `quotRem` produces. quotRemProperty1 :: E.EisensteinInteger -> E.EisensteinInteger -> Bool @@ -74,7 +74,7 @@ gcdEProperty2 :: E.EisensteinInteger -> E.EisensteinInteger -> E.EisensteinInteger -> Bool gcdEProperty2 z z1 z2 = z == 0 - || (gcd z1' z2') `rem` z == 0 + || gcd z1' z2' `rem` z == 0 where z1' = z * z1 z2' = z * z2 @@ -99,7 +99,7 @@ -- | Checks that the numbers produced by @primes@ are actually Eisenstein -- primes. primesProperty1 :: Positive Int -> Bool -primesProperty1 (Positive index) = all (isJust . isPrime . (unPrime :: Prime E.EisensteinInteger -> E.EisensteinInteger)) $ take index $ E.primes +primesProperty1 (Positive index) = all (isJust . isPrime . (unPrime :: Prime E.EisensteinInteger -> E.EisensteinInteger)) $ take index E.primes -- | Checks that the infinite list of Eisenstein primes @primes@ is ordered -- by the numbers' norm. @@ -113,7 +113,7 @@ -- sextant. primesProperty3 :: Positive Int -> Bool primesProperty3 (Positive index) = - all (\e -> abs (unPrime e) == (unPrime e :: E.EisensteinInteger)) $ take index $ E.primes + all (\e -> abs (unPrime e) == (unPrime e :: E.EisensteinInteger)) $ take index E.primes -- | An Eisenstein integer is either zero or associated (i.e. equal up to -- multiplication by a unit) to the product of its factors raised to their @@ -140,7 +140,7 @@ (factorise (15 E.:+ 12)) testSuite :: TestTree -testSuite = testGroup "EisensteinIntegers" $ +testSuite = testGroup "EisensteinIntegers" [ testSmallAndQuick "forall z . z == signum z * abs z" signumAbsProperty , testSmallAndQuick "abs z rotates to the first sextant" absProperty , testGroup "Division" diff --git a/test-suite/Math/NumberTheory/EuclideanTests.hs b/test-suite/Math/NumberTheory/EuclideanTests.hs --- a/test-suite/Math/NumberTheory/EuclideanTests.hs +++ b/test-suite/Math/NumberTheory/EuclideanTests.hs @@ -45,7 +45,7 @@ -> Bool splitIntoCoprimesProperty1 fs' = factorback fs == factorback (unCoprimes $ splitIntoCoprimes fs) where - fs = map (id *** getPower) fs' + fs = map (second getPower) fs' factorback = abs . product . map (uncurry (^)) splitIntoCoprimesProperty2 @@ -63,7 +63,7 @@ -> Bool splitIntoCoprimesProperty3 fs' = and [ coprime x y | (x : xs) <- tails fs, y <- xs ] where - fs = map fst $ unCoprimes $ splitIntoCoprimes $ map (id *** getPower) fs' + fs = map fst $ unCoprimes $ splitIntoCoprimes $ map (second getPower) fs' -- | Check that evaluation never freezes. splitIntoCoprimesProperty4 @@ -131,8 +131,8 @@ = sort (unCoprimes (splitIntoCoprimes (xs' <> ys'))) == sort (unCoprimes (splitIntoCoprimes xs' <> splitIntoCoprimes ys')) where - xs' = map (id *** getPower) xs - ys' = map (id *** getPower) ys + xs' = map (second getPower) xs + ys' = map (second getPower) ys testSuite :: TestTree testSuite = testGroup "Euclidean" diff --git a/test-suite/Math/NumberTheory/GaussianIntegersTests.hs b/test-suite/Math/NumberTheory/GaussianIntegersTests.hs --- a/test-suite/Math/NumberTheory/GaussianIntegersTests.hs +++ b/test-suite/Math/NumberTheory/GaussianIntegersTests.hs @@ -70,7 +70,7 @@ let c : _ = sqrtsModPrime (-1) p k = integerSquareRoot (unPrime p) bs = [1 .. k] - asbs = map (\b' -> ((b' * c) `mod` (unPrime p), b')) bs + asbs = map (\b' -> ((b' * c) `mod` unPrime p, b')) bs (a, b) = head [ (a', b') | (a', b') <- asbs, a' <= k] in a :+ b @@ -117,7 +117,7 @@ numberOfPrimes :: Assertion numberOfPrimes = assertEqual "counting primes: OEIS A091100" [16,100,668,4928,38404,313752] - [4 * (length $ takeWhile ((<= 10^n) . norm . unPrime) primes) | n <- [1..6]] + [4 * length (takeWhile ((<= 10^n) . norm . unPrime) primes) | n <- [1..6]] -- | signum and abs should satisfy: z == signum z * abs z signumAbsProperty :: GaussianInteger -> Bool @@ -135,7 +135,7 @@ -- | Verify that @rem@ produces a remainder smaller than the divisor with -- regards to the Euclidean domain's function. remProperty :: GaussianInteger -> GaussianInteger -> Bool -remProperty x y = (y == 0) || (norm $ x `rem` y) < (norm y) +remProperty x y = (y == 0) || norm (x `rem` y) < norm y gcdGProperty1 :: GaussianInteger -> GaussianInteger -> Bool gcdGProperty1 z1 z2 @@ -147,7 +147,7 @@ gcdGProperty2 :: GaussianInteger -> GaussianInteger -> GaussianInteger -> Bool gcdGProperty2 z z1 z2 = z == 0 - || (gcd z1' z2') `rem` z == 0 + || gcd z1' z2' `rem` z == 0 where z1' = z * z1 z2' = z * z2 @@ -160,7 +160,7 @@ gcdGSpecialCase2 = assertEqual "gcdG" (0 :+ (-1)) $ gcd (0 :+ 3) (2 :+ 2) testSuite :: TestTree -testSuite = testGroup "GaussianIntegers" $ +testSuite = testGroup "GaussianIntegers" [ testGroup "factorise" ( [ testSmallAndQuick "factor back" factoriseProperty1 , testSmallAndQuick "powers are > 0" factoriseProperty2 @@ -168,8 +168,7 @@ , testCase "factorise 63:+36" factoriseSpecialCase1 ] ++ - map (\x -> testCase "laziness" (factoriseSpecialCase2 x)) - lazyCases) + map (testCase "laziness" . factoriseSpecialCase2) lazyCases) , testSmallAndQuick "findPrime'" findPrimeProperty1 , testSmallAndQuick "isPrime" isPrimeProperty diff --git a/test-suite/Math/NumberTheory/Moduli/CbrtTests.hs b/test-suite/Math/NumberTheory/Moduli/CbrtTests.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/NumberTheory/Moduli/CbrtTests.hs @@ -0,0 +1,94 @@ +-- | +-- Module: Math.NumberTheory.Moduli.Cbrt +-- Copyright: (c) 2020 Federico Bongiorno +-- Licence: MIT +-- Maintainer: Federico Bongiorno +-- +-- Test for Math.NumberTheory.Moduli.Cbrt +-- + +{-# LANGUAGE CPP #-} + +module Math.NumberTheory.Moduli.CbrtTests + ( testSuite + ) where + +import Math.NumberTheory.Moduli.Cbrt +import Math.NumberTheory.Quadratic.EisensteinIntegers +import Math.NumberTheory.Primes +import qualified Data.Euclidean as A +import Data.List (genericReplicate) +#if __GLASGOW_HASKELL__ < 803 +import Data.Semigroup +#endif +import Test.Tasty +import Math.NumberTheory.TestUtils + +-- Checks multiplicative property of numerators. In details, +-- @cubicSymbol1 alpha1 alpha2 beta@ checks that +-- @(cubicSymbol alpha1 beta) <> (cubicSymbol alpha2 beta) == (cubicSymbol alpha1*alpha2 beta)@ +cubicSymbol1 :: EisensteinInteger -> EisensteinInteger -> EisensteinInteger -> Bool +cubicSymbol1 alpha1 alpha2 beta = isBadDenominator beta || cubicSymbolNumerator alpha1 alpha2 beta + +cubicSymbolNumerator :: EisensteinInteger -> EisensteinInteger -> EisensteinInteger -> Bool +cubicSymbolNumerator alpha1 alpha2 beta = (symbol1 <> symbol2) == symbolProduct + where + symbol1 = cubicSymbol alpha1 beta + symbol2 = cubicSymbol alpha2 beta + symbolProduct = cubicSymbol alphaProduct beta + alphaProduct = alpha1 * alpha2 + +-- Checks multiplicative property of denominators. In details, +-- @cubicSymbol2 alpha beta1 beta2@ checks that +-- @(cubicSymbol alpha beta1) <> (cubicSymbol alpha beta2) == (cubicSymbol alpha beta1*beta2)@ +cubicSymbol2 :: EisensteinInteger -> EisensteinInteger -> EisensteinInteger -> Bool +cubicSymbol2 alpha beta1 beta2 = isBadDenominator beta1 || isBadDenominator beta2 || cubicSymbolDenominator alpha beta1 beta2 + +cubicSymbolDenominator :: EisensteinInteger -> EisensteinInteger -> EisensteinInteger -> Bool +cubicSymbolDenominator alpha beta1 beta2 = (symbol1 <> symbol2) == symbolProduct + where + symbol1 = cubicSymbol alpha beta1 + symbol2 = cubicSymbol alpha beta2 + symbolProduct = cubicSymbol alpha betaProduct + betaProduct = beta1 * beta2 + +-- Checks that `cubicSymbol` agrees with the computational definition +-- +-- when the denominator is prime. +cubicSymbol3 :: EisensteinInteger -> Prime EisensteinInteger -> Bool +cubicSymbol3 alpha prime = isBadDenominator beta || cubicSymbol alpha beta == cubicSymbolPrime alpha beta + where beta = unPrime prime + +cubicSymbolPrime :: EisensteinInteger -> EisensteinInteger -> CubicSymbol +cubicSymbolPrime alpha beta = findCubicSymbol residue beta + where + residue = foldr f 1 listOfAlphas + f x y = (x * y) `A.rem` beta + listOfAlphas = genericReplicate alphaExponent alpha + -- Exponent is defined to be 1/3*(@betaNorm@ - 1). + alphaExponent = betaNorm `div` 3 + betaNorm = norm beta + +isBadDenominator :: EisensteinInteger -> Bool +isBadDenominator x = modularNorm == 0 + where + modularNorm = norm x `mod` 3 + +-- This complication is necessary because it may happen that the residue field +-- of @beta@ has characteristic two. In this case 1=-1 and the Euclidean algorithm +-- can return both. Therefore it is not enough to pattern match for the values +-- which give a well defined @cubicSymbol@. +findCubicSymbol :: EisensteinInteger -> EisensteinInteger -> CubicSymbol +findCubicSymbol residue beta + | residue `A.rem` beta == 0 = Zero + | (residue - ω) `A.rem` beta == 0 = Omega + | (residue + 1 + ω) `A.rem` beta == 0 = OmegaSquare + | (residue - 1) `A.rem` beta == 0 = One + | otherwise = error "Math.NumberTheory.Moduli.Cbrt: invalid EisensteinInteger." + +testSuite :: TestTree +testSuite = testGroup "CubicSymbol" + [ testSmallAndQuick "multiplicative property of numerators" cubicSymbol1 + , testSmallAndQuick "multiplicative property of denominators" cubicSymbol2 + , testSmallAndQuick "cubic residue with prime denominator" cubicSymbol3 + ] diff --git a/test-suite/Math/NumberTheory/Moduli/ChineseTests.hs b/test-suite/Math/NumberTheory/Moduli/ChineseTests.hs --- a/test-suite/Math/NumberTheory/Moduli/ChineseTests.hs +++ b/test-suite/Math/NumberTheory/Moduli/ChineseTests.hs @@ -7,45 +7,23 @@ -- Tests for Math.NumberTheory.Moduli.Chinese -- -{-# LANGUAGE CPP #-} -{-# LANGUAGE ViewPatterns #-} - -{-# OPTIONS_GHC -fno-warn-deprecations #-} -{-# OPTIONS_GHC -fno-warn-type-defaults #-} - module Math.NumberTheory.Moduli.ChineseTests ( testSuite ) where import Test.Tasty -import Math.NumberTheory.Moduli hiding (invertMod) +import Math.NumberTheory.Moduli (chinese) import Math.NumberTheory.TestUtils -chineseCoprimeProperty :: Integer -> Positive Integer -> Integer -> Positive Integer -> Bool -chineseCoprimeProperty n1 (Positive m1) n2 (Positive m2) = case gcd m1 m2 of - 1 -> case chineseCoprime (n1, m1) (n2, m2) of - Nothing -> False - Just n -> (n - n1) `rem` m1 == 0 && (n - n2) `rem` m2 == 0 - _ -> case chineseCoprime (n1, m1) (n2, m2) of - Nothing -> True - Just{} -> False - chineseProperty :: Integer -> Positive Integer -> Integer -> Positive Integer -> Bool -chineseProperty n1 (Positive m1) n2 (Positive m2) = if compatible - then case chinese (n1, m1) (n2, m2) of - Nothing -> False - Just n -> (n - n1) `rem` m1 == 0 && (n - n2) `rem` m2 == 0 - else case chineseCoprime (n1, m1) (n2, m2) of - Nothing -> True - Just{} -> False +chineseProperty n1 (Positive m1) n2 (Positive m2) = not compatible || + case chinese (n1, m1) (n2, m2) of + Nothing -> not compatible + Just (n, m) -> compatible && (n - n1) `rem` m1 == 0 && (n - n2) `rem` m2 == 0 && m == lcm m1 m2 where g = gcd m1 m2 compatible = (n1 - n2) `rem` g == 0 - testSuite :: TestTree -testSuite = testGroup "Chinese" - [ testSmallAndQuick "chineseCoprime" chineseCoprimeProperty - , testSmallAndQuick "chinese" chineseProperty - ] +testSuite = testSmallAndQuick "chinese" chineseProperty 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 @@ -7,9 +7,7 @@ -- Tests for Math.NumberTheory.Moduli.Class -- -{-# LANGUAGE CPP #-} {-# LANGUAGE DataKinds #-} -{-# LANGUAGE ViewPatterns #-} {-# OPTIONS_GHC -fno-warn-type-defaults #-} diff --git a/test-suite/Math/NumberTheory/Moduli/DiscreteLogarithmTests.hs b/test-suite/Math/NumberTheory/Moduli/DiscreteLogarithmTests.hs --- a/test-suite/Math/NumberTheory/Moduli/DiscreteLogarithmTests.hs +++ b/test-suite/Math/NumberTheory/Moduli/DiscreteLogarithmTests.hs @@ -1,50 +1,59 @@ {-# LANGUAGE ScopedTypeVariables #-} -{-# LANGUAGE TypeOperators #-} {-# LANGUAGE DataKinds #-} {-# OPTIONS_GHC -fno-warn-type-defaults #-} + module Math.NumberTheory.Moduli.DiscreteLogarithmTests ( testSuite ) where -import Data.Maybe -import Numeric.Natural import Test.Tasty -import Data.Semigroup + +import Data.Maybe +import Data.Mod import Data.Proxy -import GHC.TypeNats (SomeNat(..), someNatVal) +import Data.Semigroup +import GHC.TypeNats (SomeNat(..), KnownNat, someNatVal) +import Numeric.Natural import Math.NumberTheory.ArithmeticFunctions (totient) import Math.NumberTheory.Moduli.Multiplicative import Math.NumberTheory.Moduli.Singleton +import Math.NumberTheory.Primes import Math.NumberTheory.TestUtils +nextPrimitiveRoot :: (KnownNat m, UniqueFactorisation a, Integral a) => CyclicGroup a m -> Mod m -> Maybe (PrimitiveRoot m) +nextPrimitiveRoot cg g = listToMaybe $ mapMaybe (isPrimitiveRoot cg) [g..g+100] + +nextMultElement :: KnownNat m => Mod m -> Maybe (MultMod m) +nextMultElement g = listToMaybe $ mapMaybe isMultElement [g..g+100] + -- | Ensure 'discreteLogarithm' returns in the appropriate range. discreteLogRange :: Positive Natural -> Integer -> Integer -> Bool discreteLogRange (Positive m) a b = case someNatVal m of - SomeNat (_ :: Proxy m) -> fromMaybe True $ do + SomeNat (_ :: Proxy m) -> (/= Just False) $ do cg <- cyclicGroup :: Maybe (CyclicGroup Integer m) - a' <- isPrimitiveRoot cg (fromInteger a) - b' <- isMultElement (fromInteger b) + a' <- nextPrimitiveRoot cg (fromInteger a) + b' <- nextMultElement (fromInteger b) return $ discreteLogarithm cg a' b' < totient m -- | Check that 'discreteLogarithm' inverts exponentiation. discreteLogarithmProperty :: Positive Natural -> Integer -> Integer -> Bool discreteLogarithmProperty (Positive m) a b = case someNatVal m of - SomeNat (_ :: Proxy m) -> fromMaybe True $ do + SomeNat (_ :: Proxy m) -> (/= Just False) $ do cg <- cyclicGroup :: Maybe (CyclicGroup Integer m) - a' <- isPrimitiveRoot cg (fromInteger a) - b' <- isMultElement (fromInteger b) + a' <- nextPrimitiveRoot cg (fromInteger a) + b' <- nextMultElement (fromInteger b) return $ discreteLogarithm cg a' b' `stimes` unPrimitiveRoot a' == b' -- | Check that 'discreteLogarithm' inverts exponentiation in the other direction. discreteLogarithmProperty' :: Positive Natural -> Integer -> Natural -> Bool discreteLogarithmProperty' (Positive m) a k = case someNatVal m of - SomeNat (_ :: Proxy m) -> fromMaybe True $ do + SomeNat (_ :: Proxy m) -> (/= Just False) $ do cg <- cyclicGroup :: Maybe (CyclicGroup Integer m) - a'' <- isPrimitiveRoot cg (fromInteger a) + a'' <- nextPrimitiveRoot cg (fromInteger a) let a' = unPrimitiveRoot a'' return $ discreteLogarithm cg a'' (k `stimes` a') == k `mod` totient m diff --git a/test-suite/Math/NumberTheory/Moduli/JacobiTests.hs b/test-suite/Math/NumberTheory/Moduli/JacobiTests.hs --- a/test-suite/Math/NumberTheory/Moduli/JacobiTests.hs +++ b/test-suite/Math/NumberTheory/Moduli/JacobiTests.hs @@ -8,7 +8,6 @@ -- {-# LANGUAGE CPP #-} -{-# LANGUAGE ViewPatterns #-} {-# OPTIONS_GHC -fno-warn-type-defaults #-} diff --git a/test-suite/Math/NumberTheory/Moduli/PrimitiveRootTests.hs b/test-suite/Math/NumberTheory/Moduli/PrimitiveRootTests.hs --- a/test-suite/Math/NumberTheory/Moduli/PrimitiveRootTests.hs +++ b/test-suite/Math/NumberTheory/Moduli/PrimitiveRootTests.hs @@ -60,7 +60,7 @@ allUnique = go S.empty where go _ [] = True - go acc (x : xs) = if x `S.member` acc then False else go (S.insert x acc) xs + go acc (x : xs) = not (x `S.member` acc) && go (S.insert x acc) xs isPrimitiveRoot'Property1 :: forall a. (Euclidean a, Integral a, UniqueFactorisation a) diff --git a/test-suite/Math/NumberTheory/Moduli/SqrtTests.hs b/test-suite/Math/NumberTheory/Moduli/SqrtTests.hs --- a/test-suite/Math/NumberTheory/Moduli/SqrtTests.hs +++ b/test-suite/Math/NumberTheory/Moduli/SqrtTests.hs @@ -162,7 +162,7 @@ = nubOrd ps /= sort ps || nubOrd rts == sort rts where ps = map fst pes - m = product $ map (\(p, e) -> p ^ e) pes + m = product $ map (uncurry (^)) pes rts = map (`mod` m) $ take 1000 $ sqrtsModFactorisation n pes' sqrtsModFactorisationProperty3 :: AnySign Integer -> [(Prime Integer, Power Word)] -> Bool @@ -170,7 +170,7 @@ = nubOrd ps /= sort ps || all (\rt -> rt >= 0 && rt < m) rts where ps = map fst pes - m = product $ map (\(p, e) -> p ^ e) pes + m = product $ map (uncurry (^)) pes rts = take 1000 $ sqrtsModFactorisation n pes' sqrtsModFactorisationSpecialCase1 :: Assertion diff --git a/test-suite/Math/NumberTheory/Powers/ModularTests.hs b/test-suite/Math/NumberTheory/Powers/ModularTests.hs deleted file mode 100644 --- a/test-suite/Math/NumberTheory/Powers/ModularTests.hs +++ /dev/null @@ -1,91 +0,0 @@ --- | --- Module: Math.NumberTheory.Powers.ModularTests --- Copyright: (c) 2017 Andrew Lelechenko --- Licence: MIT --- Maintainer: Andrew Lelechenko --- --- Tests for Math.NumberTheory.Powers.Modular --- - -{-# LANGUAGE CPP #-} - -{-# OPTIONS_GHC -fno-warn-type-defaults #-} - -module Math.NumberTheory.Powers.ModularTests - ( testSuite - ) where - -import Test.Tasty -import Test.Tasty.HUnit - -import Numeric.Natural - -import Math.NumberTheory.Powers.Modular -import Math.NumberTheory.TestUtils - -#include "MachDeps.h" - -powMod' :: Integer -> Natural -> Integer -> Integer -powMod' = powMod - --- | Check that 'powMod' fits between 0 and m - 1. -powModProperty1 :: NonNegative Natural -> AnySign Integer -> Positive Integer -> Bool -powModProperty1 (NonNegative e) (AnySign b) (Positive m) - = let v = powMod' b e m in 0 <= v && v < m - --- | Check that 'powMod'' is multiplicative by first argument. -powModProperty2 :: NonNegative Natural -> AnySign Integer -> AnySign Integer -> Positive Integer -> Bool -powModProperty2 (NonNegative e) (AnySign b1) (AnySign b2) (Positive m) - = (powMod' b1 e m * powMod' b2 e m) `mod` m == powMod' (b1 * b2) e m - --- | Check that 'powMod' is additive by second argument. -powModProperty3 :: NonNegative Natural -> NonNegative Natural -> AnySign Integer -> Positive Integer -> Bool -powModProperty3 (NonNegative e1) (NonNegative e2) (AnySign b) (Positive m) - = (powMod' b e1 m * powMod' b e2 m) `mod` m == powMod' b (e1 + e2) m - --- | Specialized to trigger 'powModInt'. -powModProperty_Int :: AnySign Int -> NonNegative Int -> Positive Int -> Bool -powModProperty_Int (AnySign b) (NonNegative e) (Positive m) = powModInt b e m == fromInteger (powMod' (fromIntegral b) (fromIntegral e) (fromIntegral m)) - --- | Specialized to trigger 'powModWord'. -powModProperty_Word :: AnySign Word -> NonNegative Word -> Positive Word -> Bool -powModProperty_Word (AnySign b) (NonNegative e) (Positive m) = powModWord b e m == fromInteger (powMod' (fromIntegral b) (fromIntegral e) (fromIntegral m)) - --- | Specialized to trigger 'powModInteger'. -powModProperty_Integer :: AnySign Integer -> NonNegative Integer -> Positive Integer -> Bool -powModProperty_Integer (AnySign b) (NonNegative e) (Positive m) = powMod b e m == fromInteger (powMod' b (fromIntegral e) m) - --- | Specialized to trigger 'powModNatural'. -powModProperty_Natural :: AnySign Natural -> NonNegative Natural -> Positive Natural -> Bool -powModProperty_Natural (AnySign b) (NonNegative e) (Positive m) = powMod b e m == fromInteger (powMod' (fromIntegral b) e (fromIntegral m)) - -#if WORD_SIZE_IN_BITS == 64 --- | Large modulo m such that m^2 overflows. -powModSpecialCase1_Int :: Assertion -powModSpecialCase1_Int = - assertEqual "powModInt" (powModInt 3 101 (2^60-1)) 1018105167100379328 - --- | Large modulo m such that m^2 overflows. -powModSpecialCase1_Word :: Assertion -powModSpecialCase1_Word = - assertEqual "powModWord" (powModWord 3 101 (2^60-1)) 1018105167100379328 -#endif - -testSuite :: TestTree -testSuite = testGroup "Modular" - [ testGroup "powMod" - [ testSmallAndQuick "range" powModProperty1 - , testSmallAndQuick "multiplicative by base" powModProperty2 - , testSmallAndQuick "additive by exponent" powModProperty3 - - , testSmallAndQuick "powModInt" powModProperty_Int - , testSmallAndQuick "powModWord" powModProperty_Word - , testSmallAndQuick "powModInteger" powModProperty_Integer - , testSmallAndQuick "powModNatural" powModProperty_Natural - -#if WORD_SIZE_IN_BITS == 64 - , testCase "large modulo :: Int" powModSpecialCase1_Int - , testCase "large modulo :: Word" powModSpecialCase1_Word -#endif - ] - ] diff --git a/test-suite/Math/NumberTheory/Primes/FactorisationTests.hs b/test-suite/Math/NumberTheory/Primes/FactorisationTests.hs --- a/test-suite/Math/NumberTheory/Primes/FactorisationTests.hs +++ b/test-suite/Math/NumberTheory/Primes/FactorisationTests.hs @@ -32,6 +32,7 @@ , (65537^2,[(65537,2)]) , (2147483647, [(2147483647, 1)]) , (4294967291, [(4294967291, 1)]) + , (19000000000000000001, [(19000000000000000001, 1)]) , (3 * 5^2 * 7^21, [(3,1), (5,2), (7, 21)]) , (9223372036854775783, [(9223372036854775783, 1)]) , (18446744073709551557, [(18446744073709551557, 1)]) diff --git a/test-suite/Math/NumberTheory/Primes/TestingTests.hs b/test-suite/Math/NumberTheory/Primes/TestingTests.hs --- a/test-suite/Math/NumberTheory/Primes/TestingTests.hs +++ b/test-suite/Math/NumberTheory/Primes/TestingTests.hs @@ -16,8 +16,6 @@ import Test.Tasty import Test.Tasty.HUnit -import GHC.Integer.GMP.Internals (nextPrimeInteger) - import Math.NumberTheory.Primes.Testing import Math.NumberTheory.TestUtils @@ -31,28 +29,23 @@ isPrimeProperty2 n = isPrime n == isPrime (negate n) isPrimeProperty3 :: Assertion -isPrimeProperty3 = assertBool "Carmichael pseudoprimes" $ all (not . isPrime) pseudoprimes +isPrimeProperty3 = assertBool "Carmichael pseudoprimes" $ not $ any isPrime pseudoprimes where -- OEIS A002997 pseudoprimes = [561, 1105, 1729, 2465, 2821, 6601, 8911, 10585, 15841, 29341, 41041, 46657, 52633, 62745, 63973, 75361, 101101, 115921, 126217, 162401, 172081, 188461, 252601, 278545, 294409, 314821, 334153, 340561, 399001, 410041, 449065, 488881, 512461] isPrimeProperty4 :: Assertion -isPrimeProperty4 = assertBool "strong pseudoprimes to base 2" $ all (not . isPrime) pseudoprimes +isPrimeProperty4 = assertBool "strong pseudoprimes to base 2" $ not $ any isPrime pseudoprimes where -- OEIS A001262 pseudoprimes = [2047, 3277, 4033, 4681, 8321, 15841, 29341, 42799, 49141, 52633, 65281, 74665, 80581, 85489, 88357, 90751, 104653, 130561, 196093, 220729, 233017, 252601, 253241, 256999, 271951, 280601, 314821, 357761, 390937, 458989, 476971, 486737] isPrimeProperty5 :: Assertion -isPrimeProperty5 = assertBool "strong Lucas pseudoprimes" $ all (not . isPrime) pseudoprimes +isPrimeProperty5 = assertBool "strong Lucas pseudoprimes" $ not $ any isPrime pseudoprimes where -- OEIS A217255 pseudoprimes = [5459, 5777, 10877, 16109, 18971, 22499, 24569, 25199, 40309, 58519, 75077, 97439, 100127, 113573, 115639, 130139, 155819, 158399, 161027, 162133, 176399, 176471, 189419, 192509, 197801, 224369, 230691, 231703, 243629, 253259, 268349, 288919, 313499, 324899] -isPrimeProperty6 :: NonNegative Integer -> Bool -isPrimeProperty6 (NonNegative n) = if isPrime n - then nextPrimeInteger (n - 1) == n - else isPrime (nextPrimeInteger n) - isStrongFermatPPProperty :: NonNegative Integer -> Integer -> Bool isStrongFermatPPProperty (NonNegative n) b = not (isPrime n) || isStrongFermatPP n b @@ -64,7 +57,6 @@ , testCase "Carmichael pseudoprimes" isPrimeProperty3 , testCase "strong pseudoprimes base 2" isPrimeProperty4 , testCase "strong Lucas pseudoprimes" isPrimeProperty5 - , testSmallAndQuick "matches GMP" isPrimeProperty6 ] , testGroup "isStrongFermatPP" [ testSmallAndQuick "matches isPrime" isStrongFermatPPProperty diff --git a/test-suite/Math/NumberTheory/PrimesTests.hs b/test-suite/Math/NumberTheory/PrimesTests.hs --- a/test-suite/Math/NumberTheory/PrimesTests.hs +++ b/test-suite/Math/NumberTheory/PrimesTests.hs @@ -7,6 +7,10 @@ -- Tests for Math.NumberTheory.Primes -- +{-# LANGUAGE CPP #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeApplications #-} + {-# OPTIONS_GHC -fno-warn-type-defaults #-} module Math.NumberTheory.PrimesTests @@ -15,7 +19,15 @@ import Test.Tasty -import Math.NumberTheory.Primes (primes, unPrime, nextPrime, precPrime) +import Data.Bits +import Data.Int +import Data.Proxy +#if __GLASGOW_HASKELL__ < 803 +import Data.Semigroup +#endif + +import Math.NumberTheory.Primes +import qualified Math.NumberTheory.Primes.IntSet as PS import Math.NumberTheory.TestUtils primesSumWonk :: Int -> Int @@ -27,8 +39,25 @@ primesSumProperty :: NonNegative Int -> Bool primesSumProperty (NonNegative n) = n < 2 || primesSumWonk n == primesSum n +symmetricDifferenceProperty :: [Prime Int] -> [Prime Int] -> Bool +symmetricDifferenceProperty xs ys = z1 == z2 + where + x = PS.fromList xs + y = PS.fromList ys + z1 = (x PS.\\ PS.unPrimeIntSet y) <> (y PS.\\ PS.unPrimeIntSet x) + z2 = PS.symmetricDifference x y +toPrimeIntegralTest :: forall a b. (Bits a, Integral a, Bits b, Integral b) => Proxy a -> Prime b -> Bool +toPrimeIntegralTest _ p = + toIntegralSized (unPrime p) == (fmap unPrime (toPrimeIntegral p) :: Maybe a) + testSuite :: TestTree testSuite = testGroup "Primes" [ testSmallAndQuick "primesSum" primesSumProperty + , testSmallAndQuick "symmetricDifference" symmetricDifferenceProperty + , testGroup "toPrimeIntegral" + [ testSmallAndQuick "Int -> Integer" $ toPrimeIntegralTest @Integer @Int Proxy + , testSmallAndQuick "Int -> Int8" $ toPrimeIntegralTest @Int8 @Int Proxy + , testSmallAndQuick "Integer -> Int" $ toPrimeIntegralTest @Int @Integer Proxy + ] ] diff --git a/test-suite/Math/NumberTheory/Recurrences/BilinearTests.hs b/test-suite/Math/NumberTheory/Recurrences/BilinearTests.hs --- a/test-suite/Math/NumberTheory/Recurrences/BilinearTests.hs +++ b/test-suite/Math/NumberTheory/Recurrences/BilinearTests.hs @@ -30,7 +30,7 @@ binomialProperty1 (NonNegative i) = length (binomial @Integer !! i) == i + 1 binomialProperty2 :: NonNegative Int -> Bool -binomialProperty2 (NonNegative i) = binomial @Integer !! i !! 0 == 1 +binomialProperty2 (NonNegative i) = head (binomial @Integer !! i) == 1 binomialProperty3 :: NonNegative Int -> Bool binomialProperty3 (NonNegative i) = binomial @Integer !! i !! i == 1 @@ -57,10 +57,10 @@ m = m' `mod` (n + 1) binomialRotatedProperty2 :: NonNegative Int -> Bool -binomialRotatedProperty2 (NonNegative i) = binomialRotated @Integer !! i !! 0 == 1 +binomialRotatedProperty2 (NonNegative i) = head (binomialRotated @Integer !! i) == 1 binomialRotatedProperty3 :: NonNegative Int -> Bool -binomialRotatedProperty3 (NonNegative i) = binomialRotated @Integer !! 0 !! i == 1 +binomialRotatedProperty3 (NonNegative i) = head (binomialRotated @Integer) !! i == 1 binomialRotatedProperty4 :: Positive Int -> Positive Int -> Bool binomialRotatedProperty4 (Positive i) (Positive j) @@ -90,7 +90,7 @@ stirling1Property2 :: NonNegative Int -> Bool stirling1Property2 (NonNegative i) - = stirling1 !! i !! 0 + = head (stirling1 !! i) == if i == 0 then 1 else 0 stirling1Property3 :: NonNegative Int -> Bool @@ -108,7 +108,7 @@ stirling2Property2 :: NonNegative Int -> Bool stirling2Property2 (NonNegative i) - = stirling2 !! i !! 0 + = head (stirling2 !! i) == if i == 0 then 1 else 0 stirling2Property3 :: NonNegative Int -> Bool @@ -126,7 +126,7 @@ lahProperty2 :: NonNegative Int -> Bool lahProperty2 (NonNegative i) - = lah !! i !! 0 + = head (lah !! i) == product [1 .. i+1] lahProperty3 :: NonNegative Int -> Bool @@ -142,7 +142,7 @@ eulerian1Property1 (NonNegative i) = length (eulerian1 !! i) == i eulerian1Property2 :: Positive Int -> Bool -eulerian1Property2 (Positive i) = eulerian1 !! i !! 0 == 1 +eulerian1Property2 (Positive i) = head (eulerian1 !! i) == 1 eulerian1Property3 :: Positive Int -> Bool eulerian1Property3 (Positive i) = eulerian1 !! i !! (i - 1) == 1 @@ -151,15 +151,15 @@ eulerian1Property4 (Positive i) (Positive j) = j >= i - 1 || eulerian1 !! i !! j - == (toInteger $ i - j) * eulerian1 !! (i - 1) !! (j - 1) - + (toInteger j + 1) * eulerian1 !! (i - 1) !! j + == toInteger (i - j) * eulerian1 !! (i - 1) !! (j - 1) + + (toInteger j + 1) * eulerian1 !! (i - 1) !! j eulerian2Property1 :: NonNegative Int -> Bool eulerian2Property1 (NonNegative i) = length (eulerian2 !! i) == i eulerian2Property2 :: Positive Int -> Bool eulerian2Property2 (Positive i) - = eulerian2 !! i !! 0 == 1 + = head (eulerian2 !! i) == 1 eulerian2Property3 :: Positive Int -> Bool eulerian2Property3 (Positive i) @@ -170,11 +170,11 @@ eulerian2Property4 (Positive i) (Positive j) = j >= i - 1 || eulerian2 !! i !! j - == (toInteger $ 2 * i - j - 1) * eulerian2 !! (i - 1) !! (j - 1) + == toInteger (2 * i - j - 1) * eulerian2 !! (i - 1) !! (j - 1) + (toInteger j + 1) * eulerian2 !! (i - 1) !! j bernoulliSpecialCase1 :: Assertion -bernoulliSpecialCase1 = assertEqual "B_0 = 1" (bernoulli !! 0) 1 +bernoulliSpecialCase1 = assertEqual "B_0 = 1" (head bernoulli) 1 bernoulliSpecialCase2 :: Assertion bernoulliSpecialCase2 = assertEqual "B_1 = -1/2" (bernoulli !! 1) (- 1 % 2) @@ -192,7 +192,7 @@ = bernoulli !! m == (if m == 0 then 1 else 0) - sum [ bernoulli !! k - * (binomial !! m !! k % (toInteger $ m - k + 1)) + * (binomial !! m !! k % toInteger (m - k + 1)) | k <- [0 .. m - 1] ] diff --git a/test-suite/Math/NumberTheory/Recurrences/PentagonalTests.hs b/test-suite/Math/NumberTheory/Recurrences/PentagonalTests.hs --- a/test-suite/Math/NumberTheory/Recurrences/PentagonalTests.hs +++ b/test-suite/Math/NumberTheory/Recurrences/PentagonalTests.hs @@ -38,7 +38,7 @@ -- while @2@ is the @2 * 2 - 1 == 3@-rd, and so on. pentagonalNumbersProperty1 :: AnySign Int -> Bool pentagonalNumbersProperty1 (AnySign n) - | n == 0 = pents !! 0 == 0 + | n == 0 = head pents == 0 | n > 0 = pents !! (2 * n - 1) == pent n | otherwise = pents !! (2 * abs n) == pent n where diff --git a/test-suite/Math/NumberTheory/TestUtils/Wrappers.hs b/test-suite/Math/NumberTheory/TestUtils/Wrappers.hs --- a/test-suite/Math/NumberTheory/TestUtils/Wrappers.hs +++ b/test-suite/Math/NumberTheory/TestUtils/Wrappers.hs @@ -8,8 +8,6 @@ -- {-# LANGUAGE CPP #-} -{-# LANGUAGE DeriveFoldable #-} -{-# LANGUAGE DeriveFunctor #-} {-# LANGUAGE DeriveTraversable #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-} @@ -125,6 +123,7 @@ Positive l <- arbitrary ds <- vector l return $ Huge $ foldl1 (\acc n -> acc * 2^63 + n) ds + shrink (Huge n) = Huge <$> shrink n instance Eq1 Huge where liftEq eq (Huge a) (Huge b) = a `eq` b diff --git a/test-suite/Math/NumberTheory/Zeta/DirichletTests.hs b/test-suite/Math/NumberTheory/Zeta/DirichletTests.hs --- a/test-suite/Math/NumberTheory/Zeta/DirichletTests.hs +++ b/test-suite/Math/NumberTheory/Zeta/DirichletTests.hs @@ -30,7 +30,7 @@ betasOddSpecialCase1 :: Assertion betasOddSpecialCase1 = assertEqualUpToEps "beta(1) = pi/4" epsilon - (approximateValue $ betasOdd !! 0) + (approximateValue $ head betasOdd) (pi / 4) betasOddSpecialCase2 :: Assertion @@ -63,20 +63,20 @@ betasSpecialCase1 :: Assertion betasSpecialCase1 = assertEqual "beta(0) = 1/2" - (betas' !! 0) + (head betas') (1 / 2) betasSpecialCase2 :: Assertion betasSpecialCase2 = assertEqualUpToEps "beta(2) = 0.9159655" epsilon (betas' !! 2) - (0.9159655941772190150546035149323841107) + 0.9159655941772190150546035149323841107 betasSpecialCase3 :: Assertion betasSpecialCase3 = assertEqualUpToEps "beta(4) = 0.9889445" epsilon (betas' !! 4) - (0.9889445517411053361084226332283778213) + 0.9889445517411053361084226332283778213 betasProperty1 :: Positive Int -> Bool betasProperty1 (Positive m) diff --git a/test-suite/Math/NumberTheory/Zeta/RiemannTests.hs b/test-suite/Math/NumberTheory/Zeta/RiemannTests.hs --- a/test-suite/Math/NumberTheory/Zeta/RiemannTests.hs +++ b/test-suite/Math/NumberTheory/Zeta/RiemannTests.hs @@ -27,7 +27,7 @@ zetasEvenSpecialCase1 :: Assertion zetasEvenSpecialCase1 = assertEqual "zeta(0) = -1/2" - (approximateValue $ zetasEven !! 0) + (approximateValue $ head zetasEven) (-1 / 2) zetasEvenSpecialCase2 :: Assertion diff --git a/test-suite/Test.hs b/test-suite/Test.hs --- a/test-suite/Test.hs +++ b/test-suite/Test.hs @@ -9,6 +9,7 @@ import qualified Math.NumberTheory.Moduli.ChineseTests as ModuliChinese import qualified Math.NumberTheory.Moduli.ClassTests as ModuliClass +import qualified Math.NumberTheory.Moduli.CbrtTests as ModuliCbrt import qualified Math.NumberTheory.Moduli.DiscreteLogarithmTests as ModuliDiscreteLogarithm import qualified Math.NumberTheory.Moduli.EquationsTests as ModuliEquations import qualified Math.NumberTheory.Moduli.JacobiTests as ModuliJacobi @@ -18,13 +19,13 @@ import qualified Math.NumberTheory.MoebiusInversionTests as MoebiusInversion -import qualified Math.NumberTheory.Powers.ModularTests as Modular - import qualified Math.NumberTheory.PrefactoredTests as Prefactored import qualified Math.NumberTheory.PrimesTests as Primes import qualified Math.NumberTheory.Primes.CountingTests as Counting import qualified Math.NumberTheory.Primes.FactorisationTests as Factorisation +-- import qualified Math.NumberTheory.Primes.LinearAlgebraTests as LinearAlgebra +-- import qualified Math.NumberTheory.Primes.QuadraticSieveTests as QuadraticSieve import qualified Math.NumberTheory.Primes.SequenceTests as Sequence import qualified Math.NumberTheory.Primes.SieveTests as Sieve import qualified Math.NumberTheory.Primes.TestingTests as Testing @@ -48,13 +49,14 @@ import qualified Math.NumberTheory.RootsOfUnityTests as RootsOfUnity +import qualified Math.NumberTheory.DiophantineTests as Diophantine + main :: IO () main = defaultMainWithRerun tests tests :: TestTree tests = testGroup "All" - [ Modular.testSuite - , Euclidean.testSuite + [ Euclidean.testSuite , testGroup "Recurrences" [ RecurrencesPentagonal.testSuite , RecurrencesLinear.testSuite @@ -63,6 +65,7 @@ , testGroup "Moduli" [ ModuliChinese.testSuite , ModuliClass.testSuite + , ModuliCbrt.testSuite , ModuliDiscreteLogarithm.testSuite , ModuliEquations.testSuite , ModuliJacobi.testSuite @@ -76,6 +79,8 @@ [ Primes.testSuite , Counting.testSuite , Factorisation.testSuite + -- , LinearAlgebra.testSuite + -- , QuadraticSieve.testSuite , Sequence.testSuite , Sieve.testSuite , Testing.testSuite @@ -91,7 +96,7 @@ , UniqueFactorisation.testSuite , Curves.testSuite , SmoothNumbers.testSuite - + , Diophantine.testSuite , testGroup "Zeta" [ Riemann.testSuite , Dirichlet.testSuite