arithmoi 0.11.0.1 → 0.12.0.0

raw patch · 89 files changed

+1633/−1951 lines, 89 filesdep +tasty-benchdep +vector-sizeddep −gaugedep ~basedep ~integer-gmpPVP ok

version bump matches the API change (PVP)

Dependencies added: tasty-bench, vector-sized

Dependencies removed: gauge

Dependency ranges changed: base, integer-gmp

API changes (from Hackage documentation)

Files

• Math/NumberTheory/ArithmeticFunctions/Inverse.hs +5/−3
• Math/NumberTheory/ArithmeticFunctions/Mertens.hs +6/−6
• Math/NumberTheory/ArithmeticFunctions/Moebius.hs +17/−17
• Math/NumberTheory/ArithmeticFunctions/NFreedom.hs +4/−4
• Math/NumberTheory/ArithmeticFunctions/SieveBlock.hs +4/−5
• Math/NumberTheory/ArithmeticFunctions/Standard.hs +3/−3
• Math/NumberTheory/Diophantine.hs +66/−0
• Math/NumberTheory/DirichletCharacters.hs +51/−46
• Math/NumberTheory/Euclidean.hs +0/−64
• Math/NumberTheory/Euclidean/Coprimes.hs +5/−5
• Math/NumberTheory/Moduli/Cbrt.hs +160/−0
• Math/NumberTheory/Moduli/Chinese.hs +24/−118
• Math/NumberTheory/Moduli/Class.hs +0/−6
• Math/NumberTheory/Moduli/DiscreteLogarithm.hs +0/−14
• Math/NumberTheory/Moduli/Equations.hs +2/−4
• Math/NumberTheory/Moduli/Internal.hs +11/−10
• Math/NumberTheory/Moduli/Jacobi.hs +0/−18
• Math/NumberTheory/Moduli/JacobiSymbol.hs +2/−2
• Math/NumberTheory/Moduli/Multiplicative.hs +1/−2
• Math/NumberTheory/Moduli/PrimitiveRoot.hs +0/−19
• Math/NumberTheory/Moduli/Singleton.hs +18/−15
• Math/NumberTheory/Moduli/SomeMod.hs +5/−3
• Math/NumberTheory/Moduli/Sqrt.hs +77/−67
• Math/NumberTheory/MoebiusInversion.hs +1/−1
• Math/NumberTheory/Powers.hs +0/−47
• Math/NumberTheory/Powers/Cubes.hs +0/−205
• Math/NumberTheory/Powers/Fourth.hs +0/−216
• Math/NumberTheory/Powers/General.hs +0/−89
• Math/NumberTheory/Powers/Modular.hs +10/−15
• Math/NumberTheory/Powers/Squares.hs +0/−227
• Math/NumberTheory/Powers/Squares/Internal.hs +0/−136
• Math/NumberTheory/Prefactored.hs +8/−2
• Math/NumberTheory/Primes.hs +6/−6
• Math/NumberTheory/Primes/Counting/Impl.hs +40/−46
• Math/NumberTheory/Primes/Factorisation/Montgomery.hs +9/−12
• Math/NumberTheory/Primes/Factorisation/TrialDivision.hs +2/−2
• Math/NumberTheory/Primes/IntSet.hs +338/−0
• Math/NumberTheory/Primes/Sieve/Eratosthenes.hs +91/−99
• Math/NumberTheory/Primes/Sieve/Indexing.hs +3/−3
• Math/NumberTheory/Primes/Small.hs +1/−1
• Math/NumberTheory/Primes/Testing/Certified.hs +1/−1
• Math/NumberTheory/Primes/Testing/Probabilistic.hs +16/−10
• Math/NumberTheory/Primes/Types.hs +101/−6
• Math/NumberTheory/Quadratic/EisensteinIntegers.hs +32/−26
• Math/NumberTheory/Quadratic/GaussianIntegers.hs +28/−22
• Math/NumberTheory/Recurrences/Bilinear.hs +21/−27
• Math/NumberTheory/Recurrences/Linear.hs +4/−7
• Math/NumberTheory/Recurrences/Pentagonal.hs +7/−35
• Math/NumberTheory/RootsOfUnity.hs +1/−1
• Math/NumberTheory/SmoothNumbers.hs +1/−2
• Math/NumberTheory/Utils.hs +38/−13
• Math/NumberTheory/Utils/DirichletSeries.hs +0/−1
• Math/NumberTheory/Utils/FromIntegral.hs +99/−2
• Math/NumberTheory/Utils/Hyperbola.hs +9/−16
• Math/NumberTheory/Zeta/Dirichlet.hs +1/−0
• Math/NumberTheory/Zeta/Hurwitz.hs +3/−3
• Math/NumberTheory/Zeta/Utils.hs +1/−1
• arithmoi.cabal +24/−24
• benchmark/Math/NumberTheory/MertensBench.hs +0/−3
• benchmark/Math/NumberTheory/PrimitiveRootsBench.hs +1/−1
• benchmark/Math/NumberTheory/RecurrencesBench.hs +1/−1
• benchmark/Math/NumberTheory/SequenceBench.hs +2/−4
• benchmark/Math/NumberTheory/SieveBlockBench.hs +1/−1
• benchmark/Math/NumberTheory/ZetaBench.hs +2/−2
• changelog.md +21/−0
• test-suite/Math/NumberTheory/ArithmeticFunctionsTests.hs +1/−1
• test-suite/Math/NumberTheory/CurvesTests.hs +1/−3
• test-suite/Math/NumberTheory/DiophantineTests.hs +40/−0
• test-suite/Math/NumberTheory/DirichletCharactersTests.hs +2/−8
• test-suite/Math/NumberTheory/EisensteinIntegersTests.hs +5/−5
• test-suite/Math/NumberTheory/EuclideanTests.hs +4/−4
• test-suite/Math/NumberTheory/GaussianIntegersTests.hs +6/−7
• test-suite/Math/NumberTheory/Moduli/CbrtTests.hs +94/−0
• test-suite/Math/NumberTheory/Moduli/ChineseTests.hs +6/−28
• test-suite/Math/NumberTheory/Moduli/ClassTests.hs +0/−2
• test-suite/Math/NumberTheory/Moduli/DiscreteLogarithmTests.hs +22/−13
• test-suite/Math/NumberTheory/Moduli/JacobiTests.hs +0/−1
• test-suite/Math/NumberTheory/Moduli/PrimitiveRootTests.hs +1/−1
• test-suite/Math/NumberTheory/Moduli/SqrtTests.hs +2/−2
• test-suite/Math/NumberTheory/Powers/ModularTests.hs +0/−91
• test-suite/Math/NumberTheory/Primes/FactorisationTests.hs +1/−0
• test-suite/Math/NumberTheory/Primes/TestingTests.hs +3/−11
• test-suite/Math/NumberTheory/PrimesTests.hs +30/−1
• test-suite/Math/NumberTheory/Recurrences/BilinearTests.hs +13/−13
• test-suite/Math/NumberTheory/Recurrences/PentagonalTests.hs +1/−1
• test-suite/Math/NumberTheory/TestUtils/Wrappers.hs +1/−2
• test-suite/Math/NumberTheory/Zeta/DirichletTests.hs +4/−4
• test-suite/Math/NumberTheory/Zeta/RiemannTests.hs +1/−1
• test-suite/Test.hs +10/−5

@@ -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.

@@ -7,7 +7,6 @@ -- Values of <https://en.wikipedia.org/wiki/Mertens_function Mertens function>. -- -{-# 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

@@ -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

@@ -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

@@ -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#)

@@ -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 #-} 

@@ -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')

@@ -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 <https://en.wikipedia.org/wiki/Jacobi_symbol Jacobi symbol> 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

@@ -1,64 +0,0 @@--- |--- Module:      Math.NumberTheory.Euclidean--- Copyright:   (c) 2018 Alexandre Rodrigues Baldé--- Licence:     MIT--- Maintainer:  Alexandre Rodrigues Baldé <alexandrer_b@outlook.com>--- 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

@@ -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 

@@ -0,0 +1,160 @@+-- |+-- Module:      Math.NumberTheory.Moduli.Cbrt+-- Copyright:   (c) 2020 Federico Bongiorno+-- Licence:     MIT+-- Maintainer:  Federico Bongiorno <federicobongiorno97@gmail.com>+--+-- <https://en.wikipedia.org/wiki/Cubic_reciprocity#Cubic_residue_character Cubic symbol>+-- 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+-- <https://en.wikipedia.org/wiki/Cubic_reciprocity#Cubic_residue_character cubic residue character>+-- 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+-- <https://en.wikipedia.org/wiki/Cubic_reciprocity#Cubic_residue_character cubic symbol>+-- to an Eisenstein Integer.+symbolToNum :: CubicSymbol -> EisensteinInteger+symbolToNum = \case+  Zero        -> 0+  Omega       -> ω+  OmegaSquare -> -1 - ω+  One         -> 1++-- The algorithm `cubicSymbol` is adapted from+-- <https://cs.au.dk/~gudmund/Documents/cubicres.pdf here>.+-- 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.++-- | <https://en.wikipedia.org/wiki/Cubic_reciprocity#Cubic_residue_character Cubic symbol>+-- 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."

@@ -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

@@ -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

@@ -1,14 +0,0 @@--- |--- Module:       Math.NumberTheory.Moduli.DiscreteLogarithm--- Copyright:    (c) 2018 Bhavik Mehta--- License:      MIT--- Maintainer:   Andrew Lelechenko <andrew.lelechenko@gmail.com>--- Description: Deprecated-----module Math.NumberTheory.Moduli.DiscreteLogarithm-  {-# DEPRECATED "Use Math.NumberTheory.Moduli.Multiplicative instead" #-}-  ( discreteLogarithm-  ) where--import Math.NumberTheory.Moduli.Multiplicative

@@ -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-

@@ -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

@@ -1,18 +0,0 @@--- |--- Module:      Math.NumberTheory.Moduli.Jacobi--- Copyright:   (c) 2011 Daniel Fischer, 2017-2018 Andrew Lelechenko--- Licence:     MIT--- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>--- Description: Deprecated------ <https://en.wikipedia.org/wiki/Jacobi_symbol Jacobi symbol>--- 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

@@ -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,

@@ -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)

@@ -1,19 +0,0 @@--- |--- Module:      Math.NumberTheory.Moduli.PrimitiveRoot--- Copyright:   (c) 2017 Andrew Lelechenko, 2018 Bhavik Mehta--- Licence:     MIT--- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>--- 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

@@ -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)]

@@ -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

@@ -8,9 +8,11 @@ -- <https://en.wikipedia.org/wiki/Jacobi_symbol Jacobi symbol>. -- -{-# 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]

@@ -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) 

@@ -1,47 +0,0 @@--- |--- Module:      Math.NumberTheory.Powers--- Copyright:   (c) 2011 Daniel Fischer--- Licence:     MIT--- Maintainer:  Daniel Fischer <daniel.is.fischer@googlemail.com>--- 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

@@ -1,205 +0,0 @@--- |--- Module:      Math.NumberTheory.Powers.Cubes--- Copyright:   (c) 2011 Daniel Fischer--- Licence:     MIT--- Maintainer:  Daniel Fischer <daniel.is.fischer@googlemail.com>--- 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

@@ -1,216 +0,0 @@--- |--- Module:      Math.NumberTheory.Powers.Squares--- Copyright:   (c) 2011 Daniel Fischer--- Licence:     MIT--- Maintainer:  Daniel Fischer <daniel.is.fischer@googlemail.com>--- 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)

@@ -1,89 +0,0 @@--- |--- Module:      Math.NumberTheory.Powers.General--- Copyright:   (c) 2011 Daniel Fischer--- Licence:     MIT--- Maintainer:  Daniel Fischer <daniel.is.fischer@googlemail.com>--- 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)

@@ -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. --

@@ -1,227 +0,0 @@--- |--- Module:      Math.NumberTheory.Powers.Squares--- Copyright:   (c) 2011 Daniel Fischer--- Licence:     MIT--- Maintainer:  Daniel Fischer <daniel.is.fischer@googlemail.com>--- 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

@@ -1,136 +0,0 @@--- |--- Module:      Math.NumberTheory.Powers.Squares.Internal--- Copyright:   (c) 2016 Andrew Lelechenko--- Licence:     MIT--- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>--- 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. <inria-00072854>--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 #-}

@@ -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 _

@@ -5,7 +5,6 @@ -- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com> -- -{-# 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'

@@ -7,7 +7,6 @@ -- Number of primes not exceeding @n@, @&#960;(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)))

@@ -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

@@ -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)

@@ -0,0 +1,338 @@+-- |+-- Module:      Math.NumberTheory.Primes.IntSet+-- Copyright:   (c) 2020 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+--+-- 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 #-}

@@ -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)

@@ -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]

@@ -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#

@@ -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

@@ -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

@@ -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 #-}

@@ -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.

@@ -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 -- <http://www.ams.org/journals/mcom/1972-26-120/S0025-5718-1972-0314745-6/S0025-5718-1972-0314745-6.pdf Hermite-Serret algorithm>. --+-- >>> 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

@@ -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 <http://hackage.haskell.org/package/combinat combinat> 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 <http://hackage.haskell.org/package/combinat combinat> 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 <http://hackage.haskell.org/package/combinat combinat> 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]     #-}

@@ -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

@@ -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 <https://oeis.org/A000041 partition numbers>. --

@@ -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 = (<>)

@@ -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

@@ -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

@@ -9,7 +9,6 @@  {-# LANGUAGE ScopedTypeVariables   #-} {-# LANGUAGE TypeApplications      #-}-{-# LANGUAGE ViewPatterns          #-}  module Math.NumberTheory.Utils.DirichletSeries   ( DirichletSeries

@@ -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' #-}

@@ -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'

@@ -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

@@ -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

@@ -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

@@ -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 <andrew.lelechenko@gmail.com> 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

@@ -1,6 +1,3 @@-{-# LANGUAGE CPP        #-}-{-# LANGUAGE LambdaCase #-}- {-# OPTIONS_GHC -fno-warn-type-defaults #-}  module Math.NumberTheory.MertensBench

@@ -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

@@ -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

@@ -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

@@ -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   } 

@@ -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)   ]

@@ -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 

@@ -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]))

@@ -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

@@ -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+  ]

@@ -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) =

@@ -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"

@@ -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"

@@ -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

@@ -0,0 +1,94 @@+-- |+-- Module:      Math.NumberTheory.Moduli.Cbrt+-- Copyright:   (c) 2020 Federico Bongiorno+-- Licence:     MIT+-- Maintainer:  Federico Bongiorno <federicobongiorno97@gmail.com>+--+-- 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+-- <https://en.wikipedia.org/wiki/Cubic_reciprocity#Definition here>+-- 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+  ]

@@ -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

@@ -7,9 +7,7 @@ -- Tests for Math.NumberTheory.Moduli.Class -- -{-# LANGUAGE CPP             #-} {-# LANGUAGE DataKinds       #-}-{-# LANGUAGE ViewPatterns    #-}  {-# OPTIONS_GHC -fno-warn-type-defaults #-} 

@@ -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 

@@ -8,7 +8,6 @@ --  {-# LANGUAGE CPP                        #-}-{-# LANGUAGE ViewPatterns               #-}  {-# OPTIONS_GHC -fno-warn-type-defaults #-} 

@@ -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)

@@ -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

@@ -1,91 +0,0 @@--- |--- Module:      Math.NumberTheory.Powers.ModularTests--- Copyright:   (c) 2017 Andrew Lelechenko--- Licence:     MIT--- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>------ 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-    ]-  ]

@@ -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)])

@@ -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

@@ -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+    ]   ]

@@ -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]          ] 

@@ -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

@@ -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

@@ -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)

@@ -27,7 +27,7 @@ zetasEvenSpecialCase1 :: Assertion zetasEvenSpecialCase1   = assertEqual "zeta(0) = -1/2"-    (approximateValue $ zetasEven !! 0)+    (approximateValue $ head zetasEven)     (-1 / 2)  zetasEvenSpecialCase2 :: Assertion

@@ -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