arithmoi 0.6.0.1 → 0.7.0.0

raw patch · 59 files changed

+2453/−436 lines, 59 filesdep +gaugedep +vectordep −criteriondep ~QuickCheckdep ~basedep ~containersPVP ok

version bump matches the API change (PVP)

Dependencies added: gauge, vector

Dependencies removed: criterion

Dependency ranges changed: QuickCheck, base, containers, ghc-prim, tasty

API changes (from Hackage documentation)

Files

• Changes +44/−2
• GHC/TypeNats/Compat.hs +4/−3
• LICENSE +1/−1
• Math/NumberTheory/ArithmeticFunctions/Class.hs +2/−0
• Math/NumberTheory/ArithmeticFunctions/Mertens.hs +78/−0
• Math/NumberTheory/ArithmeticFunctions/Moebius.hs +176/−0
• Math/NumberTheory/ArithmeticFunctions/SieveBlock.hs +132/−0
• Math/NumberTheory/ArithmeticFunctions/SieveBlock/Unboxed.hs +132/−0
• Math/NumberTheory/ArithmeticFunctions/Standard.hs +12/−35
• Math/NumberTheory/GCD.hs +41/−1
• Math/NumberTheory/GaussianIntegers.hs +29/−33
• Math/NumberTheory/Moduli/Class.hs +24/−13
• Math/NumberTheory/Moduli/Jacobi.hs +2/−0
• Math/NumberTheory/Moduli/PrimitiveRoot.hs +148/−0
• Math/NumberTheory/MoebiusInversion.hs +8/−8
• Math/NumberTheory/MoebiusInversion/Int.hs +8/−8
• Math/NumberTheory/Powers.hs +4/−0
• Math/NumberTheory/Powers/Integer.hs +0/−44
• Math/NumberTheory/Powers/Modular.hs +97/−0
• Math/NumberTheory/Prefactored.hs +130/−0
• Math/NumberTheory/Primes/Factorisation/Montgomery.hs +107/−92
• Math/NumberTheory/Primes/Heap.hs +1/−1
• Math/NumberTheory/Primes/Sieve/Misc.hs +2/−0
• Math/NumberTheory/Primes/Testing.hs +5/−0
• Math/NumberTheory/Primes/Types.hs +53/−0
• Math/NumberTheory/SmoothNumbers.hs +130/−0
• Math/NumberTheory/UniqueFactorisation.hs +15/−45
• Math/NumberTheory/Utils.hs +17/−8
• Math/NumberTheory/Utils/FromIntegral.hs +55/−0
• Math/NumberTheory/Utils/Hyperbola.hs +92/−0
• arithmoi.cabal +34/−11
• benchmark/Bench.hs +7/−1
• benchmark/Math/NumberTheory/ArithmeticFunctionsBench.hs +3/−2
• benchmark/Math/NumberTheory/GCDBench.hs +21/−0
• benchmark/Math/NumberTheory/MertensBench.hs +21/−0
• benchmark/Math/NumberTheory/PowersBench.hs +4/−1
• benchmark/Math/NumberTheory/PrimesBench.hs +2/−2
• benchmark/Math/NumberTheory/RecurrenciesBench.hs +2/−3
• benchmark/Math/NumberTheory/SieveBlockBench.hs +95/−0
• test-suite/Math/NumberTheory/ArithmeticFunctions/MertensTests.hs +76/−0
• test-suite/Math/NumberTheory/ArithmeticFunctions/SieveBlockTests.hs +109/−0
• test-suite/Math/NumberTheory/ArithmeticFunctionsTests.hs +8/−12
• test-suite/Math/NumberTheory/GCDTests.hs +53/−0
• test-suite/Math/NumberTheory/GaussianIntegersTests.hs +20/−5
• test-suite/Math/NumberTheory/Moduli/JacobiTests.hs +2/−0
• test-suite/Math/NumberTheory/Moduli/PrimitiveRootTests.hs +118/−0
• test-suite/Math/NumberTheory/Moduli/SqrtTests.hs +24/−19
• test-suite/Math/NumberTheory/MoebiusInversion/IntTests.hs +1/−1
• test-suite/Math/NumberTheory/MoebiusInversionTests.hs +1/−1
• test-suite/Math/NumberTheory/Powers/ModularTests.hs +98/−0
• test-suite/Math/NumberTheory/PrefactoredTests.hs +101/−0
• test-suite/Math/NumberTheory/Primes/FactorisationTests.hs +16/−1
• test-suite/Math/NumberTheory/Primes/HeapTests.hs +0/−67
• test-suite/Math/NumberTheory/Primes/SieveTests.hs +1/−0
• test-suite/Math/NumberTheory/PrimesTests.hs +1/−0
• test-suite/Math/NumberTheory/TestUtils.hs +39/−4
• test-suite/Math/NumberTheory/TestUtils/Wrappers.hs +33/−9
• test-suite/Math/NumberTheory/UniqueFactorisationTests.hs +1/−1
• test-suite/Test.hs +13/−2

@@ -1,3 +1,45 @@+0.7.0.0+    This release supports GHC 7.8, 7.10, 8.0, 8.2 and 8.4.++    Breaking changes:++        Remove 'Math.NumberTheory.Powers.Integer', deprecated in 0.5.0.0.++        Deprecate 'Math.NumberTheory.Primes.Heap'.+        Use 'Math.NumberTheory.Primes.Sieve' instead.++        Deprecate 'FactorSieve', 'TotientSieve', 'CarmichaelSieve' and+        accompanying functions. Use new general approach for bulk evaluation+        of arithmetic functions instead (#77).++        Now 'moebius' returns not a number, but a value of 'Moebius' type (#90).++    New functions:++        A general framework for bulk evaluation of arithmetic functions (#77):++        > runFunctionOverBlock carmichaelA 1 10+        [1,1,2,2,4,2,6,2,6,4]++        Implement a sublinear algorithm for Mertens function (#90):++        > map (mertens . (10 ^)) [0..9]+        [1,-1,1,2,-23,-48,212,1037,1928,-222]++        Add basic support for cyclic groups and primitive roots (#86).++        Implement an efficient modular exponentiation (#86).++        Write routines for lazy generation of smooth numbers (#91).++        > smoothOverInRange (fromJust (fromList [3,5,7])) 1000 2000+        [1029,1125,1215,1225,1323,1575,1701,1715,1875]++    Improvements:++        Now factorisation of large integers and Gaussian integers produces+        factors as lazy as possible (#72, #76).+ 0.6.0.1:     Switch to smallcheck 1.1.3. @@ -96,7 +138,7 @@     Add new cabal flag check-bounds, which replaces all unsafe array functions with safe ones.      Add basic functions on Gaussian integers.-    Add Moebius mu-function.+    Add Möbius mu-function.      Forbid non-positive moduli in Math.NumberTheory.Moduli. @@ -123,7 +165,7 @@ 0.4.0.1:     Fixed Haddock bug 0.4.0.0:-    Added generalised Moebius inversion, to be continued+    Added generalised Möbius inversion, to be continued 0.3.0.0:     Added modular square roots and Chinese remainder theorem 0.2.0.6:

@@ -1,15 +1,16 @@ {-# LANGUAGE CPP         #-}  {-# OPTIONS_HADDOCK hide #-}- #if MIN_VERSION_base(4,10,0)- module GHC.TypeNats.Compat   ( module GHC.TypeNats   ) where +#if MIN_VERSION_base(4,11,0)+import GHC.TypeNats hiding (Mod)+#else import GHC.TypeNats-+#endif #else  module GHC.TypeNats.Compat

@@ -1,4 +1,4 @@-Copyright (c) 2011 Daniel Fischer, 2016-2017 Andrew Lelechenko, Carter Schonwald+Copyright (c) 2011 Daniel Fischer, 2016-2017 Andrew Lelechenko, Carter Schonwald, Google Inc.  Permission is hereby granted, free of charge, to any person obtaining a copy of this software and  associated documentation files (the "Software"), to deal in the Software without restriction,

@@ -21,7 +21,9 @@   ) where  import Control.Applicative+#if __GLASGOW_HASKELL__ < 803 import Data.Semigroup+#endif  #if MIN_VERSION_base(4,8,0) #else

@@ -0,0 +1,78 @@+-- |+-- Module:      Math.NumberTheory.ArithmeticFunctions.Mertens+-- Copyright:   (c) 2018 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Values of <https://en.wikipedia.org/wiki/Mertens_function Mertens function>.+--++{-# LANGUAGE CPP        #-}+{-# LANGUAGE LambdaCase #-}++module Math.NumberTheory.ArithmeticFunctions.Mertens+  ( mertens+  ) where++import qualified Data.Vector.Unboxed as U+#if __GLASGOW_HASKELL__ < 709+import Data.Word+#endif++import Math.NumberTheory.Powers.Cubes+import Math.NumberTheory.Powers.Squares+import Math.NumberTheory.ArithmeticFunctions.Moebius++-- | Compute individual values of Mertens function in O(n^(2/3)) time and space.+--+-- > > map (mertens . (10 ^)) [0..9]+-- > [1,-1,1,2,-23,-48,212,1037,1928,-222]+--+-- The implementation follows Theorem 3.1 from <https://arxiv.org/pdf/1610.08551.pdf Computations of the Mertens function and improved bounds on the Mertens conjecture> by G. Hurst, excluding segmentation of sieves.+mertens :: Word -> Int+mertens 0 = 0+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)++    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]+      where+        nu = integerSquareRoot y+        kappa = y `quot` (nu + 1)++    -- cacheSize ~ u+    cacheSize :: Word+    cacheSize = u `max` (x `quot` u) `max` integerSquareRoot x++    -- 1-based index+    mus :: U.Vector Moebius+    mus = sieveBlockMoebius 1 cacheSize++    lookupMus :: Word -> Moebius+    lookupMus i = U.unsafeIndex mus (fromIntegral (i - 1))++    -- 0-based index+    mes :: U.Vector Int+    mes = U.scanl' go 0 mus+      where+        go acc = \case+          MoebiusN -> acc - 1+          MoebiusZ -> acc+          MoebiusP -> acc + 1++-- | Compute sum (map (\x -> runMoebius (mu x) * f x))+sumMultMoebius :: (Word -> Moebius) -> (Word -> Int) -> [Word] -> Int+sumMultMoebius mu f = foldl go 0+  where+    go acc i = case mu i of+      MoebiusN -> acc - f i+      MoebiusZ -> acc+      MoebiusP -> acc + f i

@@ -0,0 +1,176 @@+-- |+-- Module:      Math.NumberTheory.ArithmeticFunctions.Moebius+-- Copyright:   (c) 2018 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Values of <https://en.wikipedia.org/wiki/Möbius_function Möbius function>.+--++{-# LANGUAGE BangPatterns          #-}+{-# LANGUAGE CPP                   #-}+{-# LANGUAGE MagicHash             #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE TypeFamilies          #-}++module Math.NumberTheory.ArithmeticFunctions.Moebius+  ( Moebius(..)+  , runMoebius+  , sieveBlockMoebius+  ) where++import Control.Monad (forM_, liftM)+import Control.Monad.ST (runST)+import Data.Bits+import Data.Int+import Data.Word+#if __GLASGOW_HASKELL__ < 803+import Data.Semigroup+#endif+import qualified Data.Vector.Generic         as G+import qualified Data.Vector.Generic.Mutable as M+import qualified Data.Vector.Primitive as P+import qualified Data.Vector.Unboxed         as U+import qualified Data.Vector.Unboxed.Mutable as MU+import GHC.Exts+import GHC.Integer.GMP.Internals+import Unsafe.Coerce++import Math.NumberTheory.Primes (primes)+import Math.NumberTheory.Powers.Squares (integerSquareRoot)+import Math.NumberTheory.Utils.FromIntegral (wordToInt)++import Math.NumberTheory.Logarithms++-- | Represents three possible values of <https://en.wikipedia.org/wiki/Möbius_function Möbius function>.+data Moebius+  = MoebiusN -- ^ −1+  | MoebiusZ -- ^  0+  | MoebiusP -- ^  1+  deriving (Eq, Ord, Show)++-- | Convert to any numeric type.+runMoebius :: Num a => Moebius -> a+runMoebius m = fromInteger (S# (dataToTag# m -# 1#))++fromMoebius :: Moebius -> Int8+fromMoebius m = fromIntegral $ I# (dataToTag# m)+{-# INLINE fromMoebius #-}++toMoebius :: Int8 -> Moebius+toMoebius i = let !(I# i#) = fromIntegral i in tagToEnum# i#+{-# INLINE toMoebius #-}++newtype instance U.MVector s Moebius = MV_Moebius (P.MVector s Int8)+newtype instance U.Vector    Moebius = V_Moebius  (P.Vector    Int8)++instance U.Unbox Moebius++instance M.MVector U.MVector Moebius 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_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+  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+  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++instance G.Vector U.Vector Moebius where+  {-# INLINE basicUnsafeFreeze #-}+  {-# INLINE basicUnsafeThaw #-}+  {-# INLINE basicLength #-}+  {-# 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+  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+  basicUnsafeCopy (MV_Moebius mv) (V_Moebius v) = G.basicUnsafeCopy mv v+  elemseq _ = seq++instance Semigroup Moebius where+  MoebiusZ <> _ = MoebiusZ+  _ <> MoebiusZ = MoebiusZ+  MoebiusP <> a = a+  a <> MoebiusP = a+  _ <> _ = MoebiusP++instance Monoid Moebius where+  mempty  = MoebiusP+  mappend = (<>)++-- | Evaluate the Möbius function over a block.+-- Value of @f@ at 0, if zero falls into block, is undefined.+--+-- Based on the sieving algorithm from p. 3 of <https://arxiv.org/pdf/1610.08551.pdf Computations of the Mertens function and improved bounds on the Mertens conjecture> by G. Hurst. It is approximately 5x faster than 'Math.NumberTheory.ArithmeticFunctions.SieveBlock.sieveBlockUnboxed'.+--+-- > > sieveBlockMoebius 1 10+-- > [MoebiusP, MoebiusN, MoebiusN, MoebiusZ, MoebiusN, MoebiusP, MoebiusN, MoebiusZ, MoebiusZ, MoebiusP]+sieveBlockMoebius+  :: Word+  -> Word+  -> U.Vector Moebius+sieveBlockMoebius _ 0 = U.empty+sieveBlockMoebius lowIndex' len'+  = (unsafeCoerce :: U.Vector Word8 -> U.Vector Moebius) $ runST $ do+    as <- MU.replicate len (0x80 :: Word8)+    forM_ ps $ \p -> do+      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+        MU.unsafeWrite as ix 0+    forM_ [0 .. len - 1] $ \ix -> do+      MU.unsafeModify as (mapper ix) ix+    U.unsafeFreeze as++  where+    lowIndex :: Int+    lowIndex = wordToInt lowIndex'++    len :: Int+    len = wordToInt len'++    highIndex :: Int+    highIndex = lowIndex + len - 1++    -- Bit fiddling in 'mapper' is correct only+    -- if all sufficiently small (<= 191) primes has been sieved out.+    ps :: [Int]+    ps = takeWhile (<= (191 `max` integerSquareRoot highIndex)) $ map fromInteger primes++    mapper :: Int -> Word8 -> Word8+    mapper ix val+      | val .&. 0x80 == 0x00+      = 1+      | fromIntegral (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+      | otherwise+      = ((val .&. 1) `xor` 1) `shiftL` 1

@@ -0,0 +1,132 @@+-- |+-- Module:      Math.NumberTheory.ArithmeticFunctions.SieveBlock+-- Copyright:   (c) 2017 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Bulk evaluation of arithmetic functions over continuous intervals+-- without factorisation.+--++{-# LANGUAGE BangPatterns        #-}+{-# LANGUAGE CPP                 #-}+{-# LANGUAGE MagicHash           #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE UnboxedTuples       #-}++module Math.NumberTheory.ArithmeticFunctions.SieveBlock+  ( runFunctionOverBlock+  , SieveBlockConfig(..)+  , multiplicativeSieveBlockConfig+  , additiveSieveBlockConfig+  , sieveBlock+  , sieveBlockUnboxed+  , sieveBlockMoebius+  ) where++import Control.Monad (forM_)+import Control.Monad.ST (runST)+import Data.Coerce+import qualified Data.Vector as V+import qualified Data.Vector.Mutable as MV+import GHC.Exts+#if __GLASGOW_HASKELL__ < 709+import Data.Monoid+#endif++import Math.NumberTheory.ArithmeticFunctions.Class+import Math.NumberTheory.ArithmeticFunctions.Moebius (sieveBlockMoebius)+import Math.NumberTheory.ArithmeticFunctions.SieveBlock.Unboxed+import Math.NumberTheory.Logarithms (integerLogBase')+import Math.NumberTheory.Primes (primes)+import Math.NumberTheory.Primes.Types+import Math.NumberTheory.Powers.Squares (integerSquareRoot)+import Math.NumberTheory.Utils (splitOff#)+import Math.NumberTheory.Utils.FromIntegral (wordToInt, intToWord)++-- | 'runFunctionOverBlock' @f@ @x@ @l@ evaluates an arithmetic function+-- for integers between @x@ and @x+l-1@ and returns a vector of length @l@.+-- It completely avoids factorisation, so it is asymptotically faster than+-- pointwise evaluation of @f@.+--+-- Value of @f@ at 0, if zero falls into block, is undefined.+--+-- Beware that for underlying non-commutative monoids the results may potentially+-- differ from pointwise application via 'runFunction'.+--+-- This is a thin wrapper over 'sieveBlock', read more details there.+--+-- > > runFunctionOverBlock carmichaelA 1 10+-- > [1,1,2,2,4,2,6,2,6,4]+runFunctionOverBlock+  :: ArithmeticFunction Word a+  -> Word+  -> Word+  -> V.Vector a+runFunctionOverBlock (ArithmeticFunction f g) = (V.map g .) . sieveBlock SieveBlockConfig+  { sbcEmpty                = mempty+  , sbcAppend               = mappend+  , sbcFunctionOnPrimePower = coerce f+  }++-- | Evaluate a function over a block in accordance to provided configuration.+-- Value of @f@ at 0, if zero falls into block, is undefined.+--+-- Based on Algorithm M of <https://arxiv.org/pdf/1305.1639.pdf Parity of the number of primes in a given interval and algorithms of the sublinear summation> by A. V. Lelechenko. See Lemma 2 on p. 5 on its algorithmic complexity. For the majority of use-cases its time complexity is O(x^(1+ε)).+--+-- 'sieveBlock' is similar to 'sieveBlockUnboxed' up to flavour of 'Data.Vector',+-- but is typically 7x-10x slower and consumes 3x memory.+-- Use unboxed version whenever possible.+--+-- For example, following code lists smallest prime factors:+--+-- > > sieveBlock (SieveBlockConfig maxBound min (\p _ -> p) 2 10)+-- > [2,3,2,5,2,7,2,3,2,11]+sieveBlock+  :: SieveBlockConfig a+  -> Word+  -> Word+  -> V.Vector a+sieveBlock _ _ 0 = V.empty+sieveBlock (SieveBlockConfig empty f append) lowIndex' len' = runST $ do++    let lowIndex :: Int+        lowIndex = wordToInt lowIndex'++        len :: Int+        len = wordToInt len'++    as <- V.unsafeThaw $ V.enumFromN lowIndex' len+    bs <- MV.replicate len empty++    let highIndex :: Int+        highIndex = lowIndex + len - 1++        ps :: [Int]+        ps = takeWhile (<= integerSquareRoot highIndex) $ map fromInteger primes++    forM_ ps $ \p -> do++      let p# :: Word#+          !p'@(W# p#) = intToWord p++          fs = V.generate+            (integerLogBase' (toInteger p) (toInteger highIndex))+            (\k -> f p' (intToWord k + 1))++          offset :: Int+          offset = negate lowIndex `mod` p++      forM_ [offset, offset + p .. len - 1] $ \ix -> do+        W# a# <- MV.unsafeRead as ix+        let !(# pow#, a'# #) = splitOff# p# (a# `quotWord#` p#)+        MV.unsafeWrite as ix (W# a'#)+        MV.unsafeModify bs (\y -> y `append` V.unsafeIndex fs (I# pow#)) ix++    forM_ [0 .. len - 1] $ \k -> do+      a <- MV.unsafeRead as k+      MV.unsafeModify bs (\b -> if a /= 1 then b `append` f a 1 else b) k++    V.unsafeFreeze bs

@@ -0,0 +1,132 @@+-- |+-- Module:      Math.NumberTheory.ArithmeticFunctions.SieveBlock.Unboxed+-- Copyright:   (c) 2017 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Bulk evaluation of arithmetic functions without factorisation+-- of arguments.+--++{-# LANGUAGE BangPatterns        #-}+{-# LANGUAGE MagicHash           #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE UnboxedTuples       #-}++module Math.NumberTheory.ArithmeticFunctions.SieveBlock.Unboxed+  ( SieveBlockConfig(..)+  , multiplicativeSieveBlockConfig+  , additiveSieveBlockConfig+  , sieveBlockUnboxed+  ) where++import Control.Monad (forM_)+import Control.Monad.ST (runST)+import qualified Data.Vector.Unboxed as V+import qualified Data.Vector.Unboxed.Mutable as MV+import GHC.Exts++import Math.NumberTheory.ArithmeticFunctions.Moebius (Moebius)+import Math.NumberTheory.Logarithms (integerLogBase')+import Math.NumberTheory.Primes (primes)+import Math.NumberTheory.Powers.Squares (integerSquareRoot)+import Math.NumberTheory.Utils (splitOff#)+import Math.NumberTheory.Utils.FromIntegral (wordToInt, intToWord)++-- | A record, which specifies a function to evaluate over a block.+--+-- For example, here is a configuration for the totient function:+--+-- > SieveBlockConfig+-- >   { sbcEmpty                = 1+-- >   , sbcFunctionOnPrimePower = (\p a -> (p - 1) * p ^ (a - 1)+-- >   , sbcAppend               = (*)+-- >   }+data SieveBlockConfig a = SieveBlockConfig+  { sbcEmpty                :: a+    -- ^ value of a function on 1+  , sbcFunctionOnPrimePower :: Word -> Word -> a+    -- ^ how to evaluate a function on prime powers+  , sbcAppend               :: a -> a -> a+    -- ^ how to combine values of a function on coprime arguments+  }++-- | Create a config for a multiplicative function from its definition on prime powers.+multiplicativeSieveBlockConfig :: Num a => (Word -> Word -> a) -> SieveBlockConfig a+multiplicativeSieveBlockConfig f = SieveBlockConfig+  { sbcEmpty                = 1+  , sbcFunctionOnPrimePower = f+  , sbcAppend               = (*)+  }++-- | Create a config for an additive function from its definition on prime powers.+additiveSieveBlockConfig :: Num a => (Word -> Word -> a) -> SieveBlockConfig a+additiveSieveBlockConfig f = SieveBlockConfig+  { sbcEmpty                = 0+  , sbcFunctionOnPrimePower = f+  , sbcAppend               = (+)+  }++-- | Evaluate a function over a block in accordance to provided configuration.+-- Value of @f@ at 0, if zero falls into block, is undefined.+--+-- Based on Algorithm M of <https://arxiv.org/pdf/1305.1639.pdf Parity of the number of primes in a given interval and algorithms of the sublinear summation> by A. V. Lelechenko. See Lemma 2 on p. 5 on its algorithmic complexity. For the majority of use-cases its time complexity is O(x^(1+ε)).+--+-- For example, here is an analogue of divisor function 'tau':+--+-- > > sieveBlockUnboxed (SieveBlockConfig 1 (*) (\_ a -> a + 1) 1 10)+-- > [1,2,2,3,2,4,2,4,3,4]+sieveBlockUnboxed+  :: V.Unbox a+  => SieveBlockConfig a+  -> Word+  -> Word+  -> V.Vector a+sieveBlockUnboxed _ _ 0 = V.empty+sieveBlockUnboxed (SieveBlockConfig empty f append) lowIndex' len' = runST $ do++    let lowIndex :: Int+        lowIndex = wordToInt lowIndex'++        len :: Int+        len = wordToInt len'++    as <- V.unsafeThaw $ V.enumFromN lowIndex' len+    bs <- MV.replicate len empty++    let highIndex :: Int+        highIndex = lowIndex + len - 1++        ps :: [Int]+        ps = takeWhile (<= integerSquareRoot highIndex) $ map fromInteger primes++    forM_ ps $ \p -> do++      let p# :: Word#+          !p'@(W# p#) = intToWord p++          fs = V.generate+            (integerLogBase' (toInteger p) (toInteger highIndex))+            (\k -> f p' (intToWord k + 1))++          offset :: Int+          offset = negate lowIndex `mod` p++      forM_ [offset, offset + p .. len - 1] $ \ix -> do+        W# a# <- MV.unsafeRead as ix+        let !(# pow#, a'# #) = splitOff# p# (a# `quotWord#` p#)+        MV.unsafeWrite as ix (W# a'#)+        MV.unsafeModify bs (\y -> y `append` V.unsafeIndex fs (I# pow#)) ix++    forM_ [0 .. len - 1] $ \k -> do+      a <- MV.unsafeRead as k+      MV.unsafeModify bs (\b -> if a /= 1 then b `append` f a 1 else b) k++    V.unsafeFreeze bs++{-# SPECIALIZE sieveBlockUnboxed :: SieveBlockConfig Int  -> Word -> Word -> V.Vector Int  #-}+{-# SPECIALIZE sieveBlockUnboxed :: SieveBlockConfig Word -> Word -> Word -> V.Vector Word #-}+{-# SPECIALIZE sieveBlockUnboxed :: SieveBlockConfig Bool -> Word -> Word -> V.Vector Bool #-}+{-# SPECIALIZE sieveBlockUnboxed :: SieveBlockConfig Moebius -> Word -> Word -> V.Vector Moebius #-}

@@ -26,7 +26,7 @@   , sigma, sigmaA   , totient, totientA   , jordan, jordanA-  , moebius, moebiusA+  , moebius, moebiusA, Moebius(..), runMoebius   , liouville, liouvilleA     -- * Additive functions   , additive@@ -45,7 +45,9 @@ import Data.Semigroup  import Math.NumberTheory.ArithmeticFunctions.Class+import Math.NumberTheory.ArithmeticFunctions.Moebius import Math.NumberTheory.UniqueFactorisation+import Math.NumberTheory.Utils.FromIntegral  import Numeric.Natural @@ -55,9 +57,6 @@ import Data.Word #endif -wordToInt :: Word -> Int-wordToInt = fromIntegral- -- | Create a multiplicative function from the function on prime's powers. See examples below. multiplicative :: Num a => (Prime n -> Word -> a) -> ArithmeticFunction n a multiplicative f = ArithmeticFunction ((Product .) . f) getProduct@@ -130,38 +129,38 @@ sigmaHelper pa k = (pa ^ wordToInt (k + 1) - 1) `quot` (pa - 1) {-# INLINE sigmaHelper #-} -totient :: (UniqueFactorisation n, Integral n) => n -> n+totient :: (UniqueFactorisation n, Num n) => n -> n totient = runFunction totientA  -- | Calculates the totient of a positive number @n@, i.e. --   the number of @k@ with @1 <= k <= n@ and @'gcd' n k == 1@, --   in other words, the order of the group of units in @&#8484;/(n)@.-totientA :: forall n. (UniqueFactorisation n, Integral n) => ArithmeticFunction n n+totientA :: forall n. (UniqueFactorisation n, Num n) => ArithmeticFunction n n totientA = multiplicative $ \((unPrime :: Prime n -> n) -> p) -> jordanHelper p -jordan :: (UniqueFactorisation n, Integral n) => Word -> n -> n+jordan :: (UniqueFactorisation n, Num n) => Word -> n -> n jordan = runFunction . jordanA  -- | Calculates the k-th Jordan function of an argument. -- -- > jordanA 1 = totientA-jordanA :: forall n. (UniqueFactorisation n, Integral n) => Word -> ArithmeticFunction n n+jordanA :: forall n. (UniqueFactorisation n, Num n) => Word -> ArithmeticFunction n n jordanA 0 = multiplicative $ \_ _ -> 0 jordanA 1 = totientA jordanA a = multiplicative $ \((unPrime :: Prime n -> n) -> p) -> jordanHelper (p ^ wordToInt a) -jordanHelper :: Integral n => n -> Word -> n+jordanHelper :: Num n => n -> Word -> n jordanHelper pa 1 = pa - 1 jordanHelper pa 2 = (pa - 1) * pa jordanHelper pa k = (pa - 1) * pa ^ wordToInt (k - 1) {-# INLINE jordanHelper #-} -moebius :: (UniqueFactorisation n, Num a) => n -> a+moebius :: UniqueFactorisation n => n -> Moebius moebius = runFunction moebiusA --- | Calculates the Moebius function of an argument.-moebiusA :: Num a => ArithmeticFunction n a-moebiusA = ArithmeticFunction (const f) runMoebius+-- | Calculates the Möbius function of an argument.+moebiusA :: ArithmeticFunction n Moebius+moebiusA = ArithmeticFunction (const f) id   where     f 1 = MoebiusN     f 0 = MoebiusP@@ -220,28 +219,6 @@ -- | The exponent of von Mangoldt function. Use @log expMangoldtA@ to recover von Mangoldt function itself. expMangoldtA :: forall n. (UniqueFactorisation n, Num n) => ArithmeticFunction n n expMangoldtA = ArithmeticFunction (\((unPrime :: Prime n -> n) -> p) _ -> MangoldtOne p) runMangoldt--data Moebius-  = MoebiusZ-  | MoebiusP-  | MoebiusN--runMoebius :: Num a => Moebius -> a-runMoebius m = case m of-  MoebiusZ ->  0-  MoebiusP ->  1-  MoebiusN -> -1--instance Semigroup Moebius where-  MoebiusZ <> _ = MoebiusZ-  _ <> MoebiusZ = MoebiusZ-  MoebiusP <> a = a-  a <> MoebiusP = a-  _ <> _ = MoebiusP--instance Monoid Moebius where-  mempty = MoebiusP-  mappend = (<>)  data Mangoldt a   = MangoldtZero

@@ -20,11 +20,17 @@ -- -- When using this module, always compile with optimisations turned on to -- benefit from GHC's primops and the rewrite rules.-{-# LANGUAGE CPP, BangPatterns, MagicHash #-}++{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE CPP          #-}+{-# LANGUAGE LambdaCase   #-}+{-# LANGUAGE MagicHash    #-}+ module Math.NumberTheory.GCD     ( binaryGCD     , extendedGCD     , coprime+    , splitIntoCoprimes     ) where  import Data.Bits@@ -245,3 +251,37 @@  cw32 :: Word32 -> Word32 -> Bool cw32 (W32# x#) (W32# y#) = coprimeWord# x# y#++-- | The input list is assumed to be a factorisation of some number+-- into a list of powers of (possibly, composite) non-zero factors. The output+-- list is a factorisation of the same number such that all factors+-- are coprime. Such transformation is crucial to continue factorisation+-- (lazily, in parallel or concurrent fashion) without+-- having to merge multiplicities of primes, which occurs more than in one+-- composite factor.+--+-- > > splitIntoCoprimes [(140, 1), (165, 1)]+-- > [(5,2),(28,1),(33,1)]+-- > > splitIntoCoprimes [(360, 1), (210, 1)]+-- > [(2,4),(3,3),(5,2),(7,1)]+splitIntoCoprimes :: (Integral a, Eq b, Num b) => [(a, b)] -> [(a, b)]+splitIntoCoprimes xs = if any ((== 0) . fst) ys then [(0, 1)] else go ys+  where+    ys = filter (/= (0, 0)) xs++    go = \case+      [] -> []+      ((1,  _) : rest) -> go rest+      ((x, xm) : rest) -> case popSuchThat (\(r, _) -> gcd x r /= 1) rest of+        Nothing            -> (x, xm) : go rest+        Just ((y, ym), zs) -> let g = gcd x y in+          go ((g, xm + ym) : (x `quot` g, xm) : (y `quot` g, ym) : zs)++popSuchThat :: (a -> Bool) -> [a] -> Maybe (a, [a])+popSuchThat predicate = go+  where+    go = \case+      [] -> Nothing+      (x : xs)+        | predicate x -> Just (x, xs)+        | otherwise   -> fmap (fmap (x :)) (go xs)

@@ -1,6 +1,6 @@ -- | -- Module:      Math.NumberTheory.GaussianIntegers--- Copyright:   (c) 2016 Chris Fredrickson+-- Copyright:   (c) 2016 Chris Fredrickson, Google Inc. -- Licence:     MIT -- Maintainer:  Chris Fredrickson <chris.p.fredrickson@gmail.com> -- Stability:   Provisional@@ -190,42 +190,38 @@  -- |Compute the prime factorisation of a Gaussian integer. This is unique up to units (+/- 1, +/- i). factorise :: GaussianInteger -> [(GaussianInteger, Int)]-factorise g-    | g == 0    = error "0 has no prime factorisation"-    | g == 1    = []-    | otherwise =-        let helper :: [(Integer, Int)] -> GaussianInteger -> [(GaussianInteger, Int)] -> [(GaussianInteger, Int)]-            helper [] g' fs = (if g' == 1 then [] else [(g', 1)]) ++ fs    -- include the unit, if it isn't 1-            helper ((!p, !e) : pt) g' fs-                | p `mod` 4 == 3 =-                    -- prime factors congruent to 3 mod 4 are simple.-                    let pow = div e 2-                        gp = fromInteger p-                    in helper pt (g' `divG` (gp .^ pow)) ((gp, pow) : fs)-                | otherwise      =-                    -- general case: for every prime factor of the magnitude-                    -- squared, find a Gaussian prime whose magnitude squared-                    -- is that prime. Then find out how many times the original-                    -- number is divisible by that Gaussian prime and its-                    -- conjugate. The order that the prime factors are-                    -- processed doesn't really matter, but it is reversed so-                    -- that the Gaussian primes will be in order of increasing-                    -- magnitude.-                    let gp = findPrime' p-                        (!gNext, !facs) = trialDivide g' [gp, abs $ conjugate gp] []-                    in helper pt gNext (facs ++ fs)-        in helper (reverse . Factorisation.factorise $ norm g) g []+factorise g = helper (Factorisation.factorise $ norm g) g+    where+    helper [] g' = if g' == 1 then [] else [(g', 1)] -- include the unit, if it isn't 1+    helper ((!p, !e) : pt) g' =+        -- For a given prime factor p of the magnitude squared...+        let (!g'', !facs) = if p `mod` 4 == 3+            then+                -- if the p is congruent to 3 (mod 4), then g' is divisible by+                -- p^(e/2).+                let pow = div e 2+                    gp = fromInteger p+                in (g' `divG` (gp .^ pow), [(gp, pow)])+            else+                -- otherwise: find a Gaussian prime gp for which `norm gp ==+                -- p`. Then do trial divisions to find out how many times g' is+                -- divisible by gp or its conjugate.+                let gp = findPrime' p+                in trialDivide g' [gp, abs $ conjugate gp]+        in facs ++ helper pt g''  -- Divide a Gaussian integer by a set of (relatively prime) Gaussian integers, -- as many times as possible, and return the final quotient as well as a count -- of how many times each factor divided the original.-trialDivide :: GaussianInteger -> [GaussianInteger] -> [(GaussianInteger, Int)] -> (GaussianInteger, [(GaussianInteger, Int)])-trialDivide g [] fs = (g, fs)-trialDivide g (pf : pft) fs-    | g `modG` pf == 0 =-        let (cnt, g') = countEvenDivisions g pf-        in trialDivide g' pft ((pf, cnt) : fs)-    | otherwise    = trialDivide g pft fs+trialDivide :: GaussianInteger -> [GaussianInteger] -> (GaussianInteger, [(GaussianInteger, Int)])+trialDivide = helper []+    where+    helper fs g [] = (g, fs)+    helper fs g (pf : pft)+        | g `modG` pf == 0 =+            let (cnt, g') = countEvenDivisions g pf+            in helper ((pf, cnt) : fs) g' pft+        | otherwise        = helper fs g pft  -- Divide a Gaussian integer by a possible factor, and return how many times -- the factor divided it evenly, as well as the result of dividing the original

@@ -9,14 +9,15 @@ -- Safe modular arithmetic with modulo on type level. -- -{-# LANGUAGE CPP                 #-}-{-# LANGUAGE DataKinds           #-}-{-# LANGUAGE GADTs               #-}-{-# LANGUAGE KindSignatures      #-}-{-# LANGUAGE LambdaCase          #-}-{-# LANGUAGE RankNTypes          #-}-{-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE StandaloneDeriving  #-}+{-# LANGUAGE CPP                        #-}+{-# LANGUAGE DataKinds                  #-}+{-# LANGUAGE GADTs                      #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE KindSignatures             #-}+{-# LANGUAGE LambdaCase                 #-}+{-# LANGUAGE RankNTypes                 #-}+{-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE StandaloneDeriving         #-}  module Math.NumberTheory.Moduli.Class   ( -- * Known modulo@@ -40,12 +41,18 @@ import Data.Proxy import Data.Ratio import Data.Type.Equality+import GHC.Integer.GMP.Internals+import GHC.TypeNats.Compat+ #if __GLASGOW_HASKELL__ < 709 import Data.Word+import Numeric.Natural (Natural)++powModNatural :: Natural -> Natural -> Natural -> Natural+powModNatural b e m = fromInteger $ powModInteger (toInteger b) (toInteger e) (toInteger m)+#else+import GHC.Natural (Natural, powModNatural) #endif-import GHC.Integer.GMP.Internals-import GHC.TypeNats.Compat-import Numeric.Natural  -- | Wrapper for residues modulo @m@. --@@ -60,11 +67,15 @@ -- -- Note that modulo cannot be negative. newtype Mod (m :: Nat) = Mod Natural-  deriving (Eq, Ord)+  deriving (Eq, Ord, Enum)  instance KnownNat m => Show (Mod m) where   show m = "(" ++ show (getVal m) ++ " `modulo` " ++ show (getMod m) ++ ")" +instance KnownNat m => Bounded (Mod m) where+  minBound = Mod 0+  maxBound = let mx = Mod (getNatMod mx - 1) in mx+ instance KnownNat m => Num (Mod m) where   mx@(Mod x) + Mod y =     Mod $ if xy >= m then xy - m else xy@@ -148,7 +159,7 @@ powMod :: (KnownNat m, Integral a) => Mod m -> a -> Mod m powMod mx a   | a < 0     = error $ "^{Mod}: negative exponent"-  | otherwise = Mod $ fromInteger $ powModInteger (getVal mx) (toInteger a) (getMod mx)+  | otherwise = Mod $ powModNatural (getNatVal mx) (fromIntegral a) (getNatMod mx) {-# INLINABLE [1] powMod #-}  {-# SPECIALISE [1] powMod ::

@@ -20,7 +20,9 @@  import Data.Array.Unboxed import Data.Bits+#if __GLASGOW_HASKELL__ < 803 import Data.Semigroup+#endif #if __GLASGOW_HASKELL__ < 709 import Data.Word #endif

@@ -0,0 +1,148 @@+-- |+-- Module:      Math.NumberTheory.Moduli.PrimitiveRoot+-- Copyright:   (c) 2017 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Primitive roots and cyclic groups.+--++{-# LANGUAGE CPP                  #-}+{-# LANGUAGE FlexibleContexts     #-}+{-# LANGUAGE LambdaCase           #-}+{-# LANGUAGE StandaloneDeriving   #-}+{-# LANGUAGE TupleSections        #-}+{-# LANGUAGE UndecidableInstances #-}++module Math.NumberTheory.Moduli.PrimitiveRoot+  ( isPrimitiveRoot+    -- * Cyclic groups+  , CyclicGroup(..)+  , cyclicGroupFromModulo+  , cyclicGroupToModulo+  , isPrimitiveRoot'+  ) where++#if !(MIN_VERSION_base(4,8,0))+import Control.Applicative+import Data.Word+#endif++import Math.NumberTheory.ArithmeticFunctions (totient)+import Math.NumberTheory.Moduli (Mod, getNatMod, getNatVal, KnownNat)+import Math.NumberTheory.Powers.General (highestPower)+import Math.NumberTheory.Powers.Modular+import Math.NumberTheory.Prefactored+import Math.NumberTheory.UniqueFactorisation+import Math.NumberTheory.Utils.FromIntegral++-- | A multiplicative group of residues is called cyclic,+-- if there is a primitive root @g@,+-- whose powers generates all elements.+-- Any cyclic multiplicative group of residues+-- falls into one of the following cases.+data CyclicGroup a+  = CG2 -- ^ Residues modulo 2.+  | CG4 -- ^ Residues modulo 4.+  | CGOddPrimePower       (Prime a) Word+  -- ^ Residues modulo @p@^@k@ for __odd__ prime @p@.+  | CGDoubleOddPrimePower (Prime a) Word+  -- ^ Residues modulo 2@p@^@k@ for __odd__ prime @p@.++deriving instance Show (Prime a) => Show (CyclicGroup a)++-- | Check whether a multiplicative group of residues,+-- characterized by its modulo, is cyclic and, if yes, return its form.+--+-- > > cyclicGroupFromModulo 4+-- > Just CG4+-- > > cyclicGroupFromModulo (2 * 13 ^ 3)+-- > Just (CGDoubleOddPrimePower (PrimeNat 13) 3)+-- > > cyclicGroupFromModulo (4 * 13)+-- > Nothing+cyclicGroupFromModulo+  :: (Ord a, Integral a, UniqueFactorisation a)+  => a+  -> Maybe (CyclicGroup a)+cyclicGroupFromModulo = \case+  2 -> Just CG2+  4 -> Just CG4+  n+    | n <= 1    -> Nothing+    | odd n     -> uncurry CGOddPrimePower       <$> isPrimePower n+    | odd halfN -> uncurry CGDoubleOddPrimePower <$> isPrimePower halfN+    | otherwise -> Nothing+    where+      halfN = n `quot` 2++isPrimePower+  :: (Integral a, UniqueFactorisation a)+  => a+  -> Maybe (Prime a, Word)+isPrimePower n = (, intToWord k) <$> isPrime m+  where+    (m, k) = highestPower n++-- | Extract modulo and its factorisation from+-- a cyclic multiplicative group of residues.+--+-- > > cyclicGroupToModulo CG4+-- > Prefactored {prefValue = 4, prefFactors = [(2, 2)]}+-- > > cyclicGroupToModulo (CGDoubleOddPrimePower (PrimeNat 13) 3)+-- > Prefactored {prefValue = 4394, prefFactors = [(2, 1), (13, 3)]}+cyclicGroupToModulo+  :: (Integral a, UniqueFactorisation a)+  => CyclicGroup a+  -> Prefactored a+cyclicGroupToModulo = fromFactors . \case+  CG2                       -> [(2, 1)]+  CG4                       -> [(2, 2)]+  CGOddPrimePower p k       -> [(unPrime p, k)]+  CGDoubleOddPrimePower p k -> [(2, 1), (unPrime p, k)]++-- | 'isPrimitiveRoot'' @cg@ @a@ checks whether @a@ is+-- a <https://en.wikipedia.org/wiki/Primitive_root_modulo_n primitive root>+-- of a given cyclic multiplicative group of residues @cg@.+--+-- > > let Just cg = cyclicGroupFromModulo 13+-- > > isPrimitiveRoot cg 1+-- > False+-- > > isPrimitiveRoot cg 2+-- > True+isPrimitiveRoot'+  :: (Integral a, UniqueFactorisation a)+  => CyclicGroup a+  -> a+  -> Bool+isPrimitiveRoot' cg r = r /= 0 && gcd r (prefValue m) == 1 && all (/= 1) exps+  where+    -- https://en.wikipedia.org/wiki/Primitive_root_modulo_n#Finding_primitive_roots+    m    = cyclicGroupToModulo cg+    phi  = totient m+    pows = map (\pk -> prefValue phi `quot` unPrime (fst pk)) (factorise phi)+    exps = map (\p -> powMod r p (prefValue m)) pows++-- | Check whether a given modular residue is+-- a <https://en.wikipedia.org/wiki/Primitive_root_modulo_n primitive root>.+--+-- > > isPrimitiveRoot (1 :: Mod 13)+-- > False+-- > > isPrimitiveRoot (2 :: Mod 13)+-- > True+--+-- Here is how to list all primitive roots:+--+-- > > filter isPrimitiveRoot [minBound .. maxBound] :: [Mod 13]+-- > [(2 `modulo` 13), (6 `modulo` 13), (7 `modulo` 13), (11 `modulo` 13)]+--+-- This function is a convenient wrapper around 'isPrimitiveRoot''. The latter+-- provides better control and performance, if you need them.+isPrimitiveRoot+  :: KnownNat n+  => Mod n+  -> Bool+isPrimitiveRoot r = case cyclicGroupFromModulo (getNatMod r) of+  Nothing -> False+  Just cg -> isPrimitiveRoot' cg (getNatVal r)

@@ -6,7 +6,7 @@ -- Stability:   Provisional -- Portability: Non-portable (GHC extensions) ----- Generalised Moebius inversion+-- Generalised Möbius inversion -- {-# LANGUAGE BangPatterns, FlexibleContexts #-} module Math.NumberTheory.MoebiusInversion@@ -22,7 +22,7 @@ import Math.NumberTheory.Unsafe  -- | @totientSum n@ is, for @n > 0@, the sum of @[totient k | k <- [1 .. n]]@,---   computed via generalised Moebius inversion.+--   computed via generalised Möbius inversion. --   See <http://mathworld.wolfram.com/TotientSummatoryFunction.html> for the --   formula used for @totientSum@. totientSum :: Int -> Integer@@ -32,10 +32,10 @@   where     triangle k = (k*(k+1)) `quot` 2 --- | @generalInversion g n@ evaluates the generalised Moebius inversion of @g@+-- | @generalInversion g n@ evaluates the generalised Möbius inversion of @g@ --   at the argument @n@. -----   The generalised Moebius inversion implemented here allows an efficient+--   The generalised Möbius inversion implemented here allows an efficient --   calculation of isolated values of the function @f : N -> Z@ if the function --   @g@ defined by --@@ -43,15 +43,15 @@ -- > g n = sum [f (n `quot` k) | k <- [1 .. n]] -- > -----   can be cheaply computed. By the generalised Moebius inversion formula, then+--   can be cheaply computed. By the generalised Möbius inversion formula, then -- -- > -- > f n = sum [moebius k * g (n `quot` k) | k <- [1 .. n]] -- > -- --   which allows the computation in /O/(n) steps, if the values of the---   Moebius function are known. A slightly different formula, used here,---   does not need the values of the Moebius function and allows the+--   Möbius function are known. A slightly different formula, used here,+--   does not need the values of the Möbius function and allows the --   computation in /O/(n^0.75) steps, using /O/(n^0.5) memory. -- --   An example of a pair of such functions where the inversion allows a@@ -78,7 +78,7 @@ --   method is only appropriate to compute isolated values of @f@. generalInversion :: (Int -> Integer) -> Int -> Integer generalInversion fun n-    | n < 1     = error "Moebius inversion only defined on positive domain"+    | n < 1     = error "Möbius inversion only defined on positive domain"     | n == 1    = fun 1     | n == 2    = fun 2 - fun 1     | n == 3    = fun 3 - 2*fun 1

@@ -6,7 +6,7 @@ -- Stability:   Provisional -- Portability: Non-portable (GHC extensions) ----- Generalised Moebius inversion for 'Int' valued functions.+-- Generalised Möbius inversion for 'Int' valued functions. -- {-# LANGUAGE BangPatterns, FlexibleContexts #-} {-# OPTIONS_GHC -fspec-constr-count=8 #-}@@ -23,7 +23,7 @@ import Math.NumberTheory.Unsafe  -- | @totientSum n@ is, for @n > 0@, the sum of @[totient k | k <- [1 .. n]]@,---   computed via generalised Moebius inversion.+--   computed via generalised Möbius inversion. --   See <http://mathworld.wolfram.com/TotientSummatoryFunction.html> for the --   formula used for @totientSum@. totientSum :: Int -> Int@@ -33,10 +33,10 @@   where     triangle k = (k*(k+1)) `quot` 2 --- | @generalInversion g n@ evaluates the generalised Moebius inversion of @g@+-- | @generalInversion g n@ evaluates the generalised Möbius inversion of @g@ --   at the argument @n@. -----   The generalised Moebius inversion implemented here allows an efficient+--   The generalised Möbius inversion implemented here allows an efficient --   calculation of isolated values of the function @f : N -> Z@ if the function --   @g@ defined by --@@ -44,15 +44,15 @@ -- > g n = sum [f (n `quot` k) | k <- [1 .. n]] -- > -----   can be cheaply computed. By the generalised Moebius inversion formula, then+--   can be cheaply computed. By the generalised Möbius inversion formula, then -- -- > -- > f n = sum [moebius k * g (n `quot` k) | k <- [1 .. n]] -- > -- --   which allows the computation in /O/(n) steps, if the values of the---   Moebius function are known. A slightly different formula, used here,---   does not need the values of the Moebius function and allows the+--   Möbius function are known. A slightly different formula, used here,+--   does not need the values of the Möbius function and allows the --   computation in /O/(n^0.75) steps, using /O/(n^0.5) memory. -- --   An example of a pair of such functions where the inversion allows a@@ -79,7 +79,7 @@ --   method is only appropriate to compute isolated values of @f@. generalInversion :: (Int -> Int) -> Int -> Int generalInversion fun n-    | n < 1     = error "Moebius inversion only defined on positive domain"+    | n < 1     = error "Möbius inversion only defined on positive domain"     | n == 1    = fun 1     | n == 2    = fun 2 - fun 1     | n == 3    = fun 3 - 2*fun 1

@@ -32,9 +32,13 @@   , 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,44 +0,0 @@--- |--- Module:      Math.NumberTheory.Powers.Integer--- Copyright:   (c) 2011-2014 Daniel Fischer--- Licence:     MIT--- Maintainer:  Daniel Fischer <daniel.is.fischer@googlemail.com>--- Stability:   Provisional--- Portability: Non-portable (GHC extensions)------ Potentially faster power function for 'Integer' base and 'Int'--- or 'Word' exponent.----{-# LANGUAGE CPP #-}-module Math.NumberTheory.Powers.Integer-    {-# DEPRECATED "It is no faster than (^)" #-}-    ( integerPower-    , integerWordPower-    ) where--#if MIN_VERSION_base(4,8,0)-#else-import Data.Word-#endif---- | Power of an 'Integer' by the left-to-right repeated squaring algorithm.---   This needs two multiplications in each step while the right-to-left---   algorithm needs only one multiplication for 0-bits, but here the---   two factors always have approximately the same size, which on average---   gains a bit when the result is large.------   For small results, it is unlikely to be any faster than '(^)', quite---   possibly slower (though the difference shouldn't be large), and for---   exponents with few bits set, the same holds. But for exponents with---   many bits set, the speedup can be significant.------   /Warning:/ No check for the negativity of the exponent is performed,---   a negative exponent is interpreted as a large positive exponent.-integerPower :: Integer -> Int -> Integer-integerPower = (^)-{-# DEPRECATED integerPower "Use (^) instead" #-}---- | Same as 'integerPower', but for exponents of type 'Word'.-integerWordPower :: Integer -> Word -> Integer-integerWordPower = (^)-{-# DEPRECATED integerWordPower "Use (^) instead" #-}

@@ -0,0 +1,97 @@+-- |+-- Module:      Math.NumberTheory.Powers.Modular+-- Copyright:   (c) 2017 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Modular powers (a. k. a. modular exponentiation).+--++{-# LANGUAGE CPP       #-}+{-# LANGUAGE MagicHash #-}++module Math.NumberTheory.Powers.Modular+  ( powMod+#if __GLASGOW_HASKELL__ > 709+  , powModWord+  , powModInt+#endif+  ) where++import qualified GHC.Integer.GMP.Internals as GMP (powModInteger)++#if __GLASGOW_HASKELL__ > 709+import GHC.Exts (Word(..))+import GHC.Natural (powModNatural)+import qualified GHC.Integer.GMP.Internals as GMP (powModWord)+import Math.NumberTheory.Utils.FromIntegral+#endif++-- | @powMod@ @b@ @e@ @m@ computes (@b^e@) \`mod\` @m@ in effective way.+-- An exponent @e@ must be non-negative, a modulo @m@ must be positive.+-- Otherwise the behaviour of @powMod@ is undefined.+--+-- > > powMod 2 3 5+-- > 3+-- > > powMod 3 12345678901234567890 1001+-- > 1+--+-- See also 'Math.NumberTheory.Moduli.Class.powMod' and 'Math.NumberTheory.Moduli.Class.powSomeMod'.+--+-- For finite numeric types ('Int', 'Word', etc.)+-- modulo @m@ should be such that @m^2@ does not overflow,+-- otherwise the behaviour is undefined. If you+-- 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 :: Int64)+-- > 1115647832265427613 -- incorrect due to overflow+-- > > powModInt 3 101 (2^60-1 :: Int)+-- > 1018105167100379328 -- correct+powMod :: (Integral a, Integral b) => a -> b -> a -> a+powMod x y m+  | m <= 0    = error "powModInt: non-positive modulo"+  | y <  0    = error "powModInt: negative exponent"+  | otherwise = f (x `rem` m) y 1 `mod` m+  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)++{-# INLINE [1] powMod #-}+{-# RULES+"powMod/Integer" powMod = powModInteger+  #-}++-- Work around https://ghc.haskell.org/trac/ghc/ticket/14085+powModInteger :: Integer -> Integer -> Integer -> Integer+powModInteger b e m = GMP.powModInteger (b `mod` m) e m++#if __GLASGOW_HASKELL__ > 709+{-# RULES+"powMod/Natural" powMod = powModNatural+"powMod/Word"    powMod = powModWord+"powMod/Int"     powMod = powModInt+  #-}++-- | Specialised version of 'powMod', able to handle large moduli correctly.+--+-- > > 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)++-- | Specialised version of 'powMod', able to handle large moduli correctly.+--+-- > > powModInt 3 101 (2^60-1)+-- > 1018105167100379328+powModInt :: Int -> Int -> Int -> Int+powModInt x y m+  | m <= 0 = error "powModInt: non-positive modulo"+  | y <  0 = error "powModInt: negative exponent"+  | otherwise = wordToInt $ powModWord (intToWord (x `mod` m)) (intToWord y) (intToWord m)+#endif

@@ -0,0 +1,130 @@+-- |+-- Module:      Math.NumberTheory.Prefactored+-- Copyright:   (c) 2017 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Type for numbers, accompanied by their factorisation.+--++{-# LANGUAGE CPP          #-}+{-# LANGUAGE TypeFamilies #-}++module Math.NumberTheory.Prefactored+  ( Prefactored(prefValue, prefFactors)+  , fromValue+  , fromFactors+  ) where++import Control.Arrow++import Math.NumberTheory.GCD (splitIntoCoprimes)+import Math.NumberTheory.UniqueFactorisation++#if MIN_VERSION_base(4,8,0)+#else+import Data.Word+#endif++-- | A container for a number and its pairwise coprime (but not neccessarily prime)+-- factorisation.+-- It is designed to preserve information about factors under multiplication.+-- One can use this representation to speed up prime factorisation+-- and computation of arithmetic functions.+--+-- For instance, let @p@ and @q@ be big primes:+--+-- > > let p, q :: Integer+-- > >     p = 1000000000000000000000000000057+-- > >     q = 2000000000000000000000000000071+--+-- It would be  difficult to compute prime factorisation of their product+-- as is:+-- 'factorise' would take ages. Things become different if we simply+-- change types of @p@ and @q@ to prefactored ones:+--+-- > > let p, q :: Prefactored Integer+-- > >     p = 1000000000000000000000000000057+-- > >     q = 2000000000000000000000000000071+--+-- Now prime factorisation is done instantly:+--+-- > > factorise (p * q)+-- > [(PrimeNat 1000000000000000000000000000057, 1), (PrimeNat 2000000000000000000000000000071, 1)]+-- > > factorise (p^2 * q^3)+-- > [(PrimeNat 1000000000000000000000000000057, 2), (PrimeNat 2000000000000000000000000000071, 3)]+--+-- Moreover, we can instantly compute 'totient' and its iterations.+-- It works fine, because output of 'totient' is also prefactored.+--+-- > > prefValue $ totient (p^2 * q^3)+-- > 8000000000000000000000000001752000000000000000000000000151322000000000000000000000006445392000000000000000000000135513014000000000000000000001126361040+-- > > prefValue $ totient $ totient (p^2 * q^3)+-- > 2133305798262843681544648472180210822742702690942899511234131900112583590230336435053688694839034890779375223070157301188739881477320529552945446912000+--+-- Let us look under the hood:+--+-- > > prefFactors $ totient (p^2 * q^3)+-- > [(2, 4), (41666666666666666666666666669, 1), (3, 3), (111111111111111111111111111115, 1),+-- > (1000000000000000000000000000057, 1), (2000000000000000000000000000071, 2)]+-- > > prefFactors $ totient $ totient (p^2 * q^3)+-- > [(2, 22), (39521, 1), (5, 3), (199937, 1), (3, 8), (6046667, 1), (227098769, 1),+-- > (85331809838489, 1), (361696272343, 1), (22222222222222222222222222223, 1),+-- > (41666666666666666666666666669, 1), (2000000000000000000000000000071, 1)]+--+-- Pairwise coprimality of factors is crucial, because it allows+-- us to process them independently, possibly even+-- in parallel or concurrent fashion.+--+-- Following invariant is guaranteed to hold:+--+-- > abs (prefValue x) = abs (product (map (uncurry (^)) (prefFactors x)))+data Prefactored a = Prefactored+  { prefValue   :: a+    -- ^ Number itself.+  , prefFactors :: [(a, Word)]+    -- ^ List of pairwise coprime (but not neccesarily prime) factors,+    -- accompanied by their multiplicities.+  } deriving (Show)++-- | Create 'Prefactored' from a given number.+--+-- > > fromValue 123+-- > Prefactored {prefValue = 123, prefFactors = [(123, 1)]}+fromValue :: a -> Prefactored a+fromValue a = Prefactored a [(a, 1)]++-- | Create 'Prefactored' from a given list of pairwise coprime+-- (but not neccesarily prime) factors with multiplicities.+-- If you cannot ensure coprimality, use 'splitIntoCoprimes'.+--+-- > > fromFactors (splitIntoCoprimes [(140, 1), (165, 1)])+-- > Prefactored {prefValue = 23100, prefFactors = [(5, 2), (28, 1), (33, 1)]}+-- > > fromFactors (splitIntoCoprimes [(140, 2), (165, 3)])+-- > Prefactored {prefValue = 88045650000, prefFactors = [(5, 5), (28, 2), (33, 3)]}+fromFactors :: Integral a => [(a, Word)] -> Prefactored a+fromFactors as = Prefactored (product (map (uncurry (^)) as)) (splitIntoCoprimes as)++instance (Integral a, UniqueFactorisation a) => Num (Prefactored a) where+  Prefactored v1 _ + Prefactored v2 _+    = fromValue (v1 + v2)+  Prefactored v1 _ - Prefactored v2 _+    = fromValue (v1 - v2)+  Prefactored v1 f1 * Prefactored v2 f2+    = Prefactored (v1 * v2) (splitIntoCoprimes (f1 ++ f2))+  negate (Prefactored v f) = Prefactored (negate v) f+  abs (Prefactored v f)    = Prefactored (abs v) f+  signum (Prefactored v _) = Prefactored (signum v) []+  fromInteger n = fromValue (fromInteger n)++type instance Prime (Prefactored a) = Prime a++instance UniqueFactorisation a => UniqueFactorisation (Prefactored a) where+  unPrime p = fromValue (unPrime p)+  factorise (Prefactored _ f)+    = concatMap (\(x, xm) -> map (second (* xm)) (factorise x)) f+  isPrime (Prefactored _ f) = case f of+    [(n, 1)] -> isPrime n+    _        -> Nothing

@@ -21,11 +21,12 @@ -- Given enough time, the algorithm should be able to factor numbers of 100-120 digits, but it -- is best suited for numbers of up to 50-60 digits. -{-# LANGUAGE BangPatterns   #-}-{-# LANGUAGE CPP            #-}-{-# LANGUAGE DataKinds      #-}-{-# LANGUAGE KindSignatures #-}-{-# LANGUAGE LambdaCase     #-}+{-# LANGUAGE BangPatterns        #-}+{-# LANGUAGE CPP                 #-}+{-# LANGUAGE DataKinds           #-}+{-# LANGUAGE KindSignatures      #-}+{-# LANGUAGE LambdaCase          #-}+{-# LANGUAGE ScopedTypeVariables #-}  {-# OPTIONS_GHC -fno-warn-type-defaults #-} {-# OPTIONS_HADDOCK hide #-}@@ -48,10 +49,9 @@   , findParms   ) where -#include "MachDeps.h"-+import Control.Arrow import System.Random-import Control.Monad.State.Strict+import Control.Monad.State.Lazy #if __GLASGOW_HASKELL__ < 709 import Control.Applicative import Data.Word@@ -61,10 +61,14 @@ import qualified Data.IntMap as IM import Data.List (foldl') import Data.Maybe+#if __GLASGOW_HASKELL__ < 803+import Data.Semigroup+#endif  import GHC.TypeNats.Compat  import Math.NumberTheory.Curves.Montgomery+import Math.NumberTheory.GCD (splitIntoCoprimes) import Math.NumberTheory.Moduli.Class import Math.NumberTheory.Powers.General     (highestPower, largePFPower) import Math.NumberTheory.Powers.Squares     (integerSquareRoot')@@ -138,7 +142,7 @@ stdGenFactorisation primeBound sg digits n     = curveFactorisation primeBound bailliePSW (\m -> randomR (6,m-2)) sg digits n --- | @'curveFactorisation'@ is the driver for the factorisation. Its performance (and success)+-- | '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 --   below @b@, any divisor found less than @b*b@ must be prime, thus giving @Just (b*b)@ as the --   first argument allows skipping the comparatively expensive primality test for those.@@ -150,78 +154,105 @@ --   make a huge difference. So, if the default takes too long, try another one; or you can improve your --   chances for a quick result by running several instances in parallel. -----   @'curveFactorisation'@ requires that small prime factors have been stripped before. Also, it is---   unlikely to succeed if @n@ has more than one (really) large prime factor.-curveFactorisation :: Maybe Integer                 -- ^ Lower bound for composite divisors-                   -> (Integer -> Bool)             -- ^ A primality test-                   -> (Integer -> g -> (Integer,g)) -- ^ A PRNG-                   -> g                             -- ^ Initial PRNG state-                   -> Maybe Int                     -- ^ Estimated number of digits of the smallest prime factor-                   -> Integer                       -- ^ The number to factorise-                   -> [(Integer,Int)]               -- ^ List of prime factors and exponents+--   'curveFactorisation' @n@ requires that small (< 100000) prime factors of @n@+--   have been stripped before. Otherwise it is likely to cycle forever. When in doubt,+--   use 'defaultStdGenFactorisation'.+--+--   'curveFactorisation' is unlikely to succeed if @n@ has more than one (really) large prime factor.+--+curveFactorisation+  :: forall g.+     Maybe Integer                  -- ^ Lower bound for composite divisors+  -> (Integer -> Bool)              -- ^ A primality test+  -> (Integer -> g -> (Integer, g)) -- ^ A PRNG+  -> g                              -- ^ Initial PRNG state+  -> Maybe Int                      -- ^ Estimated number of digits of the smallest prime factor+  -> Integer                        -- ^ The number to factorise+  -> [(Integer, Int)]               -- ^ List of prime factors and exponents curveFactorisation primeBound primeTest prng seed mbdigs n-    | ptest n   = [(n,1)]+    | n == 1    = []+    | ptest n   = [(n, 1)]     | otherwise = evalState (fact n digits) seed       where+        digits :: Int         digits = fromMaybe 8 mbdigs-        mult 1 xs = xs-        mult j xs = [(p,j*k) | (p,k) <- xs]-        dbl (u,v) = (mult 2 u, mult 2 v)-        ptest = case primeBound of-                  Just bd -> \k -> k <= bd || primeTest k-                  Nothing -> primeTest-        rndR k = state (\gen -> prng k gen)-        perfPw = case primeBound of-                   Nothing -> highestPower-                   Just bd -> largePFPower (integerSquareRoot' bd)-        fact m digs = do let (b1,b2,ct) = findParms digs-                         (pfs,cfs) <- repFact m b1 b2 ct-                         if null cfs-                           then return pfs-                           else do-                               nfs <- forM cfs $ \(k,j) ->-                                   mult j <$> fact k (if null pfs then digs+5 else digs)-                               return (mergeAll $ pfs:nfs)-        repFact m b1 b2 count = case perfPw m of-                                  (_,1) -> workFact m b1 b2 count-                                  (b,e)-                                    | ptest b -> return ([(b,e)],[])-                                    | otherwise -> do-                                      (as,bs) <- workFact b b1 b2 count-                                      return $ (mult e as, mult e bs)-        workFact m b1 b2 count-            | count == 0 = return ([],[(m,1)])-            | otherwise = do-                s <- rndR m-                case s `modulo` fromInteger m of-                  InfMod{} -> error "impossible case"-                  SomeMod sm -> case montgomeryFactorisation b1 b2 sm of-                    Nothing -> workFact m b1 b2 (count-1)-                    Just d  -> do-                      let !cof = m `quot` d-                      case gcd cof d of-                        1 -> do-                            (dp,dc) <- if ptest d-                                         then return ([(d,1)],[])-                                         else repFact d b1 b2 (count-1)-                            (cp,cc) <- if ptest cof-                                         then return ([(cof,1)],[])-                                         else repFact cof b1 b2 (count-1)-                            return (merge dp cp, dc ++ cc)-                        g -> do-                            let d' = d `quot` g-                                c' = cof `quot` g-                            (dp,dc) <- if ptest d'-                                         then return ([(d',1)],[])-                                         else repFact d' b1 b2 (count-1)-                            (cp,cc) <- if ptest c'-                                         then return ([(c',1)],[])-                                         else repFact c' b1 b2 (count-1)-                            (gp,gc) <- if ptest g-                                         then return ([(g,2)],[])-                                         else dbl <$> repFact g b1 b2 (count-1)-                            return  (mergeAll [dp,cp,gp], dc ++ cc ++ gc) +        ptest :: Integer -> Bool+        ptest = maybe primeTest (\bd k -> k <= bd || primeTest k) primeBound++        rndR :: Integer -> State g Integer+        rndR k = state (prng k)++        perfPw :: Integer -> (Integer, Int)+        perfPw = maybe highestPower (largePFPower . integerSquareRoot') primeBound++        fact :: Integer -> Int -> State g [(Integer, Int)]+        fact 1 _ = return mempty+        fact m digs = do+          let (b1, b2, ct) = findParms digs+          -- All factors (both @pfs@ and @cfs@), are pairwise coprime. This is+          -- because 'repFact' returns either a single factor, or output of 'workFact'.+          -- In its turn, 'workFact' returns either a single factor,+          -- or concats 'repFact's over coprime integers. Induction completes the proof.+          Factors pfs cfs <- repFact m b1 b2 ct+          case cfs of+            [] -> return pfs+            _  -> do+              nfs <- forM cfs $ \(k, j) ->+                  map (second (* j)) <$> fact k (if null pfs then digs + 5 else digs)+              return $ mconcat (pfs : nfs)++        repFact :: Integer -> Word -> Word -> Word -> State g Factors+        repFact 1 _ _ _ = return mempty+        repFact m b1 b2 count =+          case perfPw m of+            (_, 1) -> workFact m b1 b2 count+            (b, e)+              | ptest b   -> return $ singlePrimeFactor b e+              | otherwise -> modifyPowers (* e) <$> workFact b b1 b2 count++        workFact :: Integer -> Word -> Word -> Word -> State g Factors+        workFact 1 _ _ _ = return mempty+        workFact m _ _ 0 = return $ singleCompositeFactor m 1+        workFact m b1 b2 count = do+          s <- rndR m+          case s `modulo` fromInteger m of+            InfMod{} -> error "impossible case"+            SomeMod sm -> case montgomeryFactorisation b1 b2 sm of+              Nothing -> workFact m b1 b2 (count - 1)+              Just d  -> do+                let cs = splitIntoCoprimes [(d, 1), (m `quot` d, 1)]+                -- 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 $+                  \(x, xm) -> if ptest x+                              then pure $ singlePrimeFactor x xm+                              else repFact x b1 b2 (count - 1)++data Factors = Factors+  { _primeFactors     :: [(Integer, Int)]+  , _compositeFactors :: [(Integer, Int)]+  }++singlePrimeFactor :: Integer -> Int -> Factors+singlePrimeFactor a b = Factors [(a, b)] []++singleCompositeFactor :: Integer -> Int -> Factors+singleCompositeFactor a b = Factors [] [(a, b)]++instance Semigroup Factors where+  Factors pfs1 cfs1 <> Factors pfs2 cfs2+    = Factors (pfs1 <> pfs2) (cfs1 <> cfs2)++instance Monoid Factors where+  mempty = Factors [] []+  mappend = (<>)++modifyPowers :: (Int -> Int) -> Factors -> Factors+modifyPowers f (Factors pfs cfs)+  = Factors (map (second f) pfs) (map (second f) cfs)+ ---------------------------------------------------------------------------------------------------- --                                         The workhorse                                          -- ----------------------------------------------------------------------------------------------------@@ -331,22 +362,6 @@                         (k,r) | r == 1 -> ([(p,k)], Nothing)                               | otherwise -> (p,k) <: go r ps     go m [] = ([(m,1)], Nothing)---- helpers: merge sorted lists-merge :: (Ord a, Num b) => [(a, b)] -> [(a, b)] -> [(a, b)]-merge xs [] = xs-merge [] ys = ys-merge xxs@(x@(p, k) : xs) yys@(y@(q, m) : ys)-  = case p `compare` q of-    LT -> x          : merge xs yys-    EQ -> (p, k + m) : merge xs  ys-    GT -> y          : merge xxs ys--mergeAll :: (Ord a, Num b) => [[(a, b)]] -> [(a, b)]-mergeAll = \case-  []              -> []-  [xs]            -> xs-  (xs : ys : zss) -> merge (merge xs ys) (mergeAll zss)  -- | For a given estimated decimal length of the smallest prime factor -- ("tier") return parameters B1, B2 and the number of curves to try

@@ -17,7 +17,7 @@ {-# LANGUAGE BangPatterns, CPP, MonoLocalBinds #-} {-# OPTIONS_GHC -funbox-strict-fields #-} {-# OPTIONS_GHC -fno-float-in -fno-spec-constr -fno-full-laziness #-}-module Math.NumberTheory.Primes.Heap (primes, sieveFrom) where+module Math.NumberTheory.Primes.Heap {-# DEPRECATED "Use Math.NumberTheory.Primes.Sieve instead" #-} (primes, sieveFrom) where  import Data.Array.Unboxed import Data.Array.ST

@@ -46,6 +46,8 @@ import Math.NumberTheory.Unsafe import Math.NumberTheory.Utils +{-# DEPRECATED FactorSieve, TotientSieve, CarmichaelSieve, factorSieve, sieveFactor, fsBound, fsPrimeTest, totientSieve, sieveTotient, carmichaelSieve, sieveCarmichael "Use new interface for sieves, provided by \"Math.NumberTheory.ArithmeticFunctions.SieveBlock\"" #-}+ {- IMPORTANT NOTICE: Not all sieves use the same layout! 

@@ -7,6 +7,9 @@ -- Portability: Non-portable (GHC extensions) -- -- Primality tests.++{-# OPTIONS_GHC -fno-warn-deprecations #-}+ module Math.NumberTheory.Primes.Testing     ( -- * Standard tests       isPrime@@ -27,6 +30,8 @@ import Math.NumberTheory.Primes.Testing.Certified import Math.NumberTheory.Primes.Factorisation.TrialDivision import Math.NumberTheory.Primes.Sieve.Misc++{-# DEPRECATED fsIsPrime "Use new interface for sieves, provided by Math.NumberTheory.ArithmeticFunctions.SieveBlock" #-}  -- | Test primality using a 'FactorSieve'. If @n@ is out of bounds --   of the sieve, fall back to 'isPrime'.

@@ -0,0 +1,53 @@+-- |+-- Module:      Math.NumberTheory.Primes.Types+-- Copyright:   (c) 2017 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- This is an internal module, defining types for primes.+-- Should not be exposed to users.+--++{-# LANGUAGE CPP          #-}+{-# LANGUAGE TypeFamilies #-}++module Math.NumberTheory.Primes.Types+  ( Prime+  , Prm(..)+  , PrimeNat(..)+  ) where++#if MIN_VERSION_base(4,8,0)+#else+import Data.Word+#endif++import Numeric.Natural++newtype Prm = Prm { unPrm :: Word }+  deriving (Eq, Ord)++instance Show Prm where+  showsPrec d (Prm p) r = (if d > 10 then "(" ++ s ++ ")" else s) ++ r+    where+      s = "Prm " ++ show p++newtype PrimeNat = PrimeNat { unPrimeNat :: Natural }+  deriving (Eq, Ord)++instance Show PrimeNat where+  showsPrec d (PrimeNat p) r = (if d > 10 then "(" ++ s ++ ")" else s) ++ r+    where+      s = "PrimeNat " ++ show p++-- | Type of primes of a given unique factorisation domain.+--+-- @abs (unPrime n) == unPrime n@ must hold for all @n@ of type @Prime t@+type family Prime (f :: *) :: *++type instance Prime Int     = Prm+type instance Prime Word    = Prm+type instance Prime Integer = PrimeNat+type instance Prime Natural = PrimeNat

@@ -0,0 +1,130 @@+-- |+-- Module:      Math.NumberTheory.SmoothNumbers+-- Copyright:   (c) 2018 Frederick Schneider+-- Licence:     MIT+-- Maintainer:  Frederick Schneider <frederick.schneider2011@gmail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- A <https://en.wikipedia.org/wiki/Smooth_number smooth number>+-- is an integer, which can be represented as a product of powers of elements+-- from a given set (smooth basis). E. g., 48 = 3 * 4 * 4 is smooth+-- over a set {3, 4}, and 24 is not.+--++module Math.NumberTheory.SmoothNumbers+  ( -- * Create a smooth basis+    SmoothBasis+  , fromSet+  , fromList+  , fromSmoothUpperBound+    -- * Generate smooth numbers+  , smoothOver+  , smoothOverInRange+  , smoothOverInRangeBF+  ) where++import Data.List (sort, nub)+import qualified Data.Set as S+import Math.NumberTheory.Primes.Sieve (primes)++-- | An abstract representation of a smooth basis.+-- It consists of a set of coprime numbers ≥2.+newtype SmoothBasis a = SmoothBasis { unSmoothBasis :: [a] } deriving Show++-- | Build a 'SmoothBasis' from a set of coprime numbers ≥2.+--+-- > > fromSet (Set.fromList [2, 3])+-- > Just (SmoothBasis [2, 3])+-- > > fromSet (Set.fromList [2, 4]) -- should be coprime+-- > Nothing+-- > > fromSet (Set.fromList [1, 3]) -- should be >= 2+-- > Nothing+fromSet :: Integral a => S.Set a -> Maybe (SmoothBasis a)+fromSet s = if isValid l then Just (SmoothBasis l) else Nothing where l = S.elems s++-- | Build a 'SmoothBasis' from a list of coprime numbers ≥2.+--+-- > > fromList [2, 3]+-- > Just (SmoothBasis [2, 3])+-- > > fromList [2, 2]+-- > Just (SmoothBasis [2])+-- > > fromList [2, 4] -- should be coprime+-- > Nothing+-- > > fromList [1, 3] -- should be >= 2+-- > Nothing+fromList :: Integral a => [a] -> Maybe (SmoothBasis a)+fromList l = if isValid l' then Just (SmoothBasis l') else Nothing where l' = nub $ sort l++-- | Build a 'SmoothBasis' from a list of primes below given bound.+--+-- > > fromSmoothUpperBound 10+-- > Just (SmoothBasis [2, 3, 5, 7])+-- > > fromSmoothUpperBound 1+-- > Nothing+fromSmoothUpperBound :: Integral a => a -> Maybe (SmoothBasis a)+fromSmoothUpperBound n = if (n < 2)+                         then Nothing+                         else Just $ SmoothBasis $ map fromInteger $ takeWhile (<= nI) primes+                         where nI = toInteger n++-- | Generate an infinite ascending list of+-- <https://en.wikipedia.org/wiki/Smooth_number smooth numbers>+-- over a given smooth basis.+--+-- > > take 10 (smoothOver (fromJust (fromList [2, 5])))+-- > [1, 2, 4, 5, 8, 10, 16, 20, 25, 32]+smoothOver :: Integral a => SmoothBasis a -> [a]+smoothOver pl = foldr (\p l -> mergeListLists $ iterate (map (p*)) l) [1] (unSmoothBasis pl)+                where+                      {-# INLINE mergeListLists #-}+                      mergeListLists      = foldr go1 []+                        where go1 (h:t) b = h:(go2 t b)+                              go1 _     b = b++                              go2 a@(ah:at) b@(bh:bt)+                                | bh < ah   = bh : (go2 a bt)+                                | otherwise = ah : (go2 at b) -- no possibility of duplicates+                              go2 a b = if null a then b else a++-- | Generate an ascending list of+-- <https://en.wikipedia.org/wiki/Smooth_number smooth numbers>+-- over a given smooth basis in a given range.+--+-- It may appear inefficient+-- for short, but distant ranges;+-- consider using 'smoothOverInRangeBF' in such cases.+--+-- > > smoothOverInRange (fromJust (fromList [2, 5])) 100 200+-- > [100, 125, 128, 160, 200]+smoothOverInRange   :: Integral a => SmoothBasis a -> a -> a -> [a]+smoothOverInRange s lo hi = takeWhile (<= hi) $ dropWhile (< lo) (smoothOver s)++-- | Generate an ascending list of+-- <https://en.wikipedia.org/wiki/Smooth_number smooth numbers>+-- over a given smooth basis in a given range.+--+-- It is inefficient+-- for large or starting near 0 ranges;+-- consider using 'smoothOverInRange' in such cases.+--+-- Suffix BF stands for the brute force algorithm, involving a lot of divisions.+--+-- > > smoothOverInRangeBF (fromJust (fromList [2, 5])) 100 200+-- > [100, 125, 128, 160, 200]+smoothOverInRangeBF :: Integral a => SmoothBasis a -> a -> a -> [a]+smoothOverInRangeBF prs lo hi = filter (mf prs') [lo..hi]+                                where mf []     n    = (n == 1) -- mf means manually factor+                                      mf pl@(p:ps) n = if (mod n p == 0)+                                                       then mf pl (div n p)+                                                       else mf ps n+                                      prs'           = unSmoothBasis prs++-- | isValid assumes that the list is sorted and unique and then checks if the list is suitable to be a SmoothBasis.+isValid :: (Integral a) => [a] -> Bool+isValid pl = if (length pl == 0) then False else v' pl+             where v' []        = True+                   v' (x:xs)    = if (x < 2 || (not $ rpl x xs)) then False else v' xs+                   rpl _ []     = True  -- rpl means relatively prime to the rest of the list+                   rpl n (x:xs) = if (gcd n x > 1) then False else rpl n xs+

@@ -26,29 +26,13 @@ #endif  import Math.NumberTheory.Primes.Factorisation as F (factorise')-import Math.NumberTheory.GaussianIntegers as G+import Math.NumberTheory.Primes.Testing.Probabilistic as T (isPrime)+import Math.NumberTheory.Primes.Types (Prime, Prm(..), PrimeNat(..))+import qualified Math.NumberTheory.GaussianIntegers as G+import Math.NumberTheory.Utils.FromIntegral  import Numeric.Natural -newtype SmallPrime = SmallPrime { _unSmallPrime :: Word }-  deriving (Eq, Ord, Show)--newtype BigPrime = BigPrime { _unBigPrime :: Natural }-  deriving (Eq, Ord, Show)---- | Type of primes of a given unique factorisation domain.--- When the domain has exactly one unit, @Prime t = t@,--- but when units are multiple more restricted types--- (or at least newtypes) should be specified.------ @abs (unPrime n) == unPrime n@ must hold for all @n@ of type @Prime t@-type family Prime (f :: *) :: *--type instance Prime Int     = SmallPrime-type instance Prime Word    = SmallPrime-type instance Prime Integer = BigPrime-type instance Prime Natural = BigPrime- type instance Prime G.GaussianInteger = GaussianPrime  -- | The following invariant must hold for @n /= 0@:@@ -60,6 +44,13 @@   unPrime   :: Prime a -> a   factorise :: a -> [(Prime a, Word)] +  isPrime   :: a -> Maybe (Prime a)+  isPrime n = case factorise n of+    [(p, 1)] -> Just p+    _        -> Nothing++  {-# MINIMAL unPrime, factorise #-}+ instance UniqueFactorisation Int where   unPrime   = coerce wordToInt   factorise m' = if m <= 1@@ -67,13 +58,14 @@                 else map (coerce integerToWord *** intToWord) . F.factorise' . intToInteger $ m                   where                     m = abs m'+  isPrime n = if T.isPrime (toInteger n) then Just (coerce $ intToWord $ abs n) else Nothing  instance UniqueFactorisation Word where   unPrime     = coerce   factorise m = if m <= 1                   then []                   else map (coerce integerToWord *** intToWord) . F.factorise' . wordToInteger $ m-+  isPrime n = if T.isPrime (toInteger n) then Just (coerce n) else Nothing  instance UniqueFactorisation Integer where   unPrime      = coerce naturalToInteger@@ -82,12 +74,14 @@                 else map (coerce integerToNatural *** intToWord) . F.factorise' $ m                   where                     m = abs m'+  isPrime n = if T.isPrime n then Just (coerce $ integerToNatural $ abs n) else Nothing  instance UniqueFactorisation Natural where   unPrime   = coerce   factorise m = if m <= 1                   then []                   else map (coerce integerToNatural *** intToWord) . F.factorise' . naturalToInteger $ m+  isPrime n = if T.isPrime (toInteger n) then Just (coerce n) else Nothing  newtype GaussianPrime = GaussianPrime { _unGaussianPrime :: G.GaussianInteger }   deriving (Eq, Show)@@ -97,27 +91,3 @@    factorise 0 = []   factorise g = map (coerce *** intToWord) $ filter (\(h, _) -> abs h /= 1) $ G.factorise g---------------- Utils--wordToInt :: Word -> Int-wordToInt = fromIntegral--wordToInteger :: Word -> Integer-wordToInteger = fromIntegral--intToWord :: Int -> Word-intToWord = fromIntegral--intToInteger :: Int -> Integer-intToInteger = fromIntegral--naturalToInteger :: Natural -> Integer-naturalToInteger = fromIntegral--integerToNatural :: Integer -> Natural-integerToNatural = fromIntegral--integerToWord :: Integer -> Word-integerToWord = fromIntegral

@@ -20,6 +20,7 @@     , bitCountWord#     , uncheckedShiftR     , splitOff+    , splitOff#     ) where  #include "MachDeps.h"@@ -191,14 +192,22 @@                 case (v2 `plusWord#` uncheckedShiftRL# v2 4#) `and#` mf## of                   v3 -> word2Int# (uncheckedShiftRL# (v3 `timesWord#` m1##) sd#) --- {-# SPECIALISE splitOff :: Integer -> Integer -> (Int, Integer),---                            Int -> Int -> (Int, Int),---                            Word -> Word -> (Int, Word)---   #-}-{-# INLINABLE splitOff #-}-splitOff :: Integral a => a -> a -> (Int, a)+splitOff :: Integer -> Integer -> (Int, Integer)+splitOff _ 0 = (0, 0) -- prevent infinite loop splitOff p n = go 0 n   where     go !k m = case m `quotRem` p of-                (q,r) | r == 0 -> go (k+1) q-                      | otherwise -> (k,m)+      (q, 0) -> go (k + 1) q+      _      -> (k, m)+{-# INLINABLE splitOff #-}++-- | It is difficult to convince GHC to unbox output of 'splitOff' and 'splitOff.go',+-- so we fallback to a specialized unboxed version to minimize allocations.+splitOff# :: Word# -> Word# -> (# Int#, Word# #)+splitOff# _ 0## = (# 0#, 0## #)+splitOff# p n = go 0# n+  where+    go k m = case m `quotRemWord#` p of+      (# q, 0## #) -> go (k +# 1#) q+      _            -> (# k, m #)+{-# INLINABLE splitOff# #-}

@@ -0,0 +1,55 @@+-- |+-- Module:      Math.NumberTheory.Utils.FromIntegral+-- Copyright:   (c) 2017 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Monomorphic `fromIntegral`.+--++{-# LANGUAGE CPP                 #-}++module Math.NumberTheory.Utils.FromIntegral+  ( wordToInt+  , wordToInteger+  , intToWord+  , intToInteger+  , naturalToInteger+  , integerToNatural+  , integerToWord+  ) where++import Numeric.Natural+#if __GLASGOW_HASKELL__ < 709+import Data.Word+#endif++wordToInt :: Word -> Int+wordToInt = fromIntegral+{-# INLINE wordToInt #-}++wordToInteger :: Word -> Integer+wordToInteger = fromIntegral+{-# INLINE wordToInteger #-}++intToWord :: Int -> Word+intToWord = fromIntegral+{-# INLINE intToWord #-}++intToInteger :: Int -> Integer+intToInteger = fromIntegral+{-# INLINE intToInteger #-}++naturalToInteger :: Natural -> Integer+naturalToInteger = fromIntegral+{-# INLINE naturalToInteger #-}++integerToNatural :: Integer -> Natural+integerToNatural = fromIntegral+{-# INLINE integerToNatural #-}++integerToWord :: Integer -> Word+integerToWord = fromIntegral+{-# INLINE integerToWord #-}

@@ -0,0 +1,92 @@+-- |+-- Module:      Math.NumberTheory.Utils.Hyperbola+-- Copyright:   (c) 2018 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Highest points under hyperbola.+--++{-# LANGUAGE CPP        #-}+{-# LANGUAGE LambdaCase #-}++module Math.NumberTheory.Utils.Hyperbola+  ( pointsUnderHyperbola+  ) where++import Data.Bits++import Math.NumberTheory.Powers.Cubes++-- | Straightforward computation of+-- [ n `quot` x | x <- [hi, hi - 1 .. lo] ].+-- Unfortunately, such list generator performs poor,+-- so we fall back to manual recursion.+pointsUnderHyperbola0 :: Int -> Int -> Int -> [Int]+pointsUnderHyperbola0 n lo hi+  | n < 0     = error "pointsUnderHyperbola0: first argument must be non-negative"+  | lo <= 0   = error "pointsUnderHyperbola0: second argument must be positive"+  | otherwise = go hi+    where+      go x+        | x < lo    = []+        | otherwise = n `quot` x : go (x - 1)++data Bresenham = Bresenham+  {  bresX       :: !Int+  ,  bresBeta    :: !Int+  , _bresGamma   :: !Int+  , _bresDelta1  :: !Int+  , _bresEpsilon :: !Int+  }++initBresenham :: Int -> Int -> Bresenham+initBresenham n x = Bresenham x beta gamma delta1 epsilon+  where+    beta    = n `quot` x+    epsilon = n `rem` x+    delta1  = n `quot` (x - 1) - beta+    gamma   = beta - (x - 1) * delta1++-- | 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)))+  where+    x       = x' - 1+    eps     = epsilon' + gamma'+{-# INLINE stepBack #-}++-- | Division-free computation of+-- [ n `quot` x | x <- [hi, hi - 1 .. lo] ].+-- In other words, we compute y-coordinates of highest integral points+-- under hyperbola @x * y = n@ between @x = lo@ and @x = hi@ in reverse order.+--+-- The implementation follows section 5 of <https://arxiv.org/pdf/1206.3369.pdf A successive approximation algorithm for computing the divisor summatory function>+-- by R. Sladkey.+-- It is 2x faster than a trivial implementation for 'Int'.+pointsUnderHyperbola :: Int -> Int -> Int -> [Int]+pointsUnderHyperbola n lo hi+  | n < 0        = error "pointsUnderHyperbola: first argument must be non-negative"+  | lo <= 0      = error "pointsUnderHyperbola: second argument must be positive"+  | hi <  lo     = []+  | hi == lo     = [n `quot` lo]+  | otherwise    = go (initBresenham n hi)+  where+    mid = (integerCubeRoot (2 * n) + 1) `max` lo+    go h+      | bresX h < mid = pointsUnderHyperbola0 n lo ((mid - 1) `min` hi)+      | otherwise = bresBeta h : go (stepBack h)

@@ -1,8 +1,8 @@ name                : arithmoi-version             : 0.6.0.1+version             : 0.7.0.0 cabal-version       : >= 1.10 author              : Daniel Fischer-copyright           : (c) 2011 Daniel Fischer, 2016-2017 Andrew Lelechenko, Carter Schonwald+copyright           : (c) 2011 Daniel Fischer, 2016-2018 Andrew Lelechenko, Carter Schonwald license             : MIT license-file        : LICENSE maintainer          : Carter Schonwald  carter at wellposed dot com,@@ -42,6 +42,7 @@                         , mtl >= 2.0 && < 2.3                         , exact-pi >= 0.4.1.1                         , integer-logarithms >= 1.0+                        , vector >= 0.12     if impl(ghc < 7.10)       build-depends     : nats >= 1 && <1.2     if impl(ghc < 8.0)@@ -49,12 +50,16 @@      exposed-modules     : Math.NumberTheory.ArithmeticFunctions                           Math.NumberTheory.ArithmeticFunctions.Class+                          Math.NumberTheory.ArithmeticFunctions.Mertens+                          Math.NumberTheory.ArithmeticFunctions.Moebius+                          Math.NumberTheory.ArithmeticFunctions.SieveBlock                           Math.NumberTheory.ArithmeticFunctions.Standard                           Math.NumberTheory.Curves.Montgomery                           Math.NumberTheory.Moduli                           Math.NumberTheory.Moduli.Chinese                           Math.NumberTheory.Moduli.Class                           Math.NumberTheory.Moduli.Jacobi+                          Math.NumberTheory.Moduli.PrimitiveRoot                           Math.NumberTheory.Moduli.Sqrt                           Math.NumberTheory.MoebiusInversion                           Math.NumberTheory.MoebiusInversion.Int@@ -69,7 +74,7 @@                           Math.NumberTheory.Powers.Cubes                           Math.NumberTheory.Powers.Fourth                           Math.NumberTheory.Powers.General-                          Math.NumberTheory.Powers.Integer+                          Math.NumberTheory.Powers.Modular                           Math.NumberTheory.Primes                           Math.NumberTheory.Primes.Sieve                           Math.NumberTheory.Primes.Factorisation@@ -80,8 +85,14 @@                           Math.NumberTheory.Primes.Heap                           Math.NumberTheory.UniqueFactorisation                           Math.NumberTheory.Zeta+                          Math.NumberTheory.Prefactored                           GHC.TypeNats.Compat-    other-modules       : Math.NumberTheory.Utils+                          Math.NumberTheory.SmoothNumbers++    other-modules       : Math.NumberTheory.ArithmeticFunctions.SieveBlock.Unboxed+                          Math.NumberTheory.Utils+                          Math.NumberTheory.Utils.FromIntegral+                          Math.NumberTheory.Utils.Hyperbola                           Math.NumberTheory.Unsafe                           Math.NumberTheory.Primes.Counting.Impl                           Math.NumberTheory.Primes.Counting.Approximate@@ -93,7 +104,8 @@                           Math.NumberTheory.Primes.Testing.Probabilistic                           Math.NumberTheory.Primes.Testing.Certified                           Math.NumberTheory.Primes.Testing.Certificates.Internal-    other-extensions          : BangPatterns+                          Math.NumberTheory.Primes.Types+    other-extensions    : BangPatterns      ghc-options         : -O2 -Wall     ghc-prof-options    : -O2 -auto@@ -108,16 +120,22 @@ benchmark criterion   build-depends:    base                     , arithmoi-                    , criterion+                    , gauge                     , containers                     , random                     , integer-logarithms+                    , vector   if impl(ghc < 7.10)     build-depends     : nats >= 1 && <1.2+  if impl(ghc < 8.0)+    build-depends     : semigroups >= 0.8   other-modules:    Math.NumberTheory.ArithmeticFunctionsBench+                  , Math.NumberTheory.GCDBench+                  , Math.NumberTheory.MertensBench                   , Math.NumberTheory.PowersBench                   , Math.NumberTheory.PrimesBench                   , Math.NumberTheory.RecurrenciesBench+                  , Math.NumberTheory.SieveBlockBench   hs-source-dirs:   benchmark   main-is:          Bench.hs   type:             exitcode-stdio-1.0@@ -132,13 +150,14 @@   build-depends:        arithmoi                       , base >= 4.6 && < 5                       , containers >= 0.5 && < 0.6-                      , tasty >= 0.10 && < 0.12+                      , tasty >= 0.10 && < 1.1                       , tasty-smallcheck >= 0.8 && < 0.9-                      , tasty-quickcheck >= 0.9 && < 0.10-                      , tasty-hunit >= 0.9 && < 0.10-                      , QuickCheck >= 2.10 && < 2.11+                      , tasty-quickcheck >= 0.9 && < 0.11+                      , tasty-hunit >= 0.9 && < 0.11+                      , QuickCheck >= 2.10 && < 2.12                       , transformers >= 0.5                       , integer-gmp < 1.1+                      , vector   if impl(ghc < 7.10)     build-depends     : smallcheck >= 1.1 && < 1.1.3,                         nats >= 1 && <1.2@@ -148,6 +167,8 @@     build-depends     : semigroups >= 0.8    other-modules :   Math.NumberTheory.ArithmeticFunctionsTests+                  , Math.NumberTheory.ArithmeticFunctions.MertensTests+                  , Math.NumberTheory.ArithmeticFunctions.SieveBlockTests                   , Math.NumberTheory.CurvesTests                   , Math.NumberTheory.GaussianIntegersTests                   , Math.NumberTheory.GCDTests@@ -157,17 +178,19 @@                   , Math.NumberTheory.Moduli.ChineseTests                   , Math.NumberTheory.Moduli.ClassTests                   , Math.NumberTheory.Moduli.JacobiTests+                  , Math.NumberTheory.Moduli.PrimitiveRootTests                   , Math.NumberTheory.Moduli.SqrtTests                   , Math.NumberTheory.Powers.CubesTests                   , Math.NumberTheory.MoebiusInversionTests                   , Math.NumberTheory.MoebiusInversion.IntTests                   , Math.NumberTheory.Powers.FourthTests                   , Math.NumberTheory.Powers.GeneralTests+                  , Math.NumberTheory.Powers.ModularTests                   , Math.NumberTheory.Powers.SquaresTests+                  , Math.NumberTheory.PrefactoredTests                   , Math.NumberTheory.PrimesTests                   , Math.NumberTheory.Primes.CountingTests                   , Math.NumberTheory.Primes.FactorisationTests-                  , Math.NumberTheory.Primes.HeapTests                   , Math.NumberTheory.Primes.SieveTests                   , Math.NumberTheory.Primes.TestingTests                   , Math.NumberTheory.TestUtils

@@ -1,15 +1,21 @@ module Main where -import Criterion.Main+import Gauge.Main  import Math.NumberTheory.ArithmeticFunctionsBench as ArithmeticFunctions+import Math.NumberTheory.GCDBench as GCD+import Math.NumberTheory.MertensBench as Mertens import Math.NumberTheory.PowersBench as Powers import Math.NumberTheory.PrimesBench as Primes import Math.NumberTheory.RecurrenciesBench as Recurrencies+import Math.NumberTheory.SieveBlockBench as SieveBlock  main = defaultMain   [ ArithmeticFunctions.benchSuite+  , GCD.benchSuite+  , Mertens.benchSuite   , Powers.benchSuite   , Primes.benchSuite   , Recurrencies.benchSuite+  , SieveBlock.benchSuite   ]

@@ -2,7 +2,7 @@   ( benchSuite   ) where -import Criterion.Main+import Gauge.Main import Data.Set (Set)  import Math.NumberTheory.ArithmeticFunctions as A@@ -13,12 +13,13 @@ compareSetFunctions :: String -> (Integer -> Set Integer) -> Benchmark compareSetFunctions name new = bench name $ nf (map new) [1..100000] +benchSuite :: Benchmark benchSuite = bgroup "ArithmeticFunctions"   [ compareSetFunctions "divisors" A.divisors   , bench "divisors/int" $ nf (map A.divisorsSmall) [1 :: Int .. 100000]   , compareFunctions "totient" A.totient   , compareFunctions "carmichael" A.carmichael-  , compareFunctions "moebius" A.moebius+  , compareFunctions "moebius" (A.runMoebius . A.moebius)   , compareFunctions "tau" A.tau   , compareFunctions "sigma 1" (A.sigma 1)   , compareFunctions "sigma 2" (A.sigma 2)

@@ -0,0 +1,21 @@+module Math.NumberTheory.GCDBench+  ( benchSuite+  ) where++import Gauge.Main++import Math.NumberTheory.GCD as A+import Prelude as P++benchSuite = bgroup "GCD"+  [ subSuite "large coprimes" 1073741823 100003+  , subSuite "powers of 2" (2^12) (2^19)+  , subSuite "power of 23" (23^3) (23^7)+  ]+  where subSuite :: String -> Int -> Int -> Benchmark+        subSuite name m n = bgroup name+          [ bench "Prelude.gcd" $ nf (P.gcd m) n+          , bench "binaryGCD" $ nf (A.binaryGCD m) n+          , bench "Prelude.coprime" $ nf (\n -> 1 == P.gcd m n) n+          , bench "coprime" $ nf (A.coprime m) n+          ]

@@ -0,0 +1,21 @@+{-# LANGUAGE CPP        #-}+{-# LANGUAGE LambdaCase #-}++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.MertensBench+  ( benchSuite+  ) where++import Gauge.Main+#if __GLASGOW_HASKELL__ < 709+import Data.Word+#endif++import Math.NumberTheory.ArithmeticFunctions.Mertens++mertensBench :: Word -> Benchmark+mertensBench n = bench (show n) (nf mertens n)++benchSuite :: Benchmark+benchSuite = bgroup "Mertens" $ map mertensBench $ take 4 $ iterate (* 10) 10000000

@@ -1,8 +1,10 @@+{-# OPTIONS_GHC -fno-warn-type-defaults #-}+ module Math.NumberTheory.PowersBench   ( benchSuite   ) where -import Criterion.Main+import Gauge.Main import System.Random  import Math.NumberTheory.Logarithms (integerLog2)@@ -25,4 +27,5 @@   where     n = genInteger 0 bits +benchSuite :: Benchmark benchSuite = bgroup "Powers" $ map compareRoots [2300, 2400 .. 2600]

@@ -1,10 +1,10 @@-{-# OPTIONS_GHC -fno-warn-warnings-deprecations #-}+{-# OPTIONS_GHC -fno-warn-type-defaults #-}  module Math.NumberTheory.PrimesBench   ( benchSuite   ) where -import Criterion.Main+import Gauge.Main import System.Random  import Math.NumberTheory.Logarithms (integerLog2)

@@ -4,9 +4,7 @@   ( benchSuite   ) where -import Criterion.Main-import Numeric.Natural-import System.Random+import Gauge.Main  import Math.NumberTheory.Recurrencies.Bilinear @@ -21,6 +19,7 @@     benchAt i j = bench ("!! " ++ show i ++ " !! " ++ show j)                 $ nf (\(x, y) -> triangle !! x !! y :: Integer) (i, j) +benchSuite :: Benchmark benchSuite = bgroup "Bilinear"   [ benchTriangle "binomial"  binomial 1000   , benchTriangle "stirling1" stirling1 100

@@ -0,0 +1,95 @@+{-# LANGUAGE CPP        #-}+{-# LANGUAGE LambdaCase #-}++{-# OPTIONS_GHC -fno-warn-deprecations  #-}+{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.SieveBlockBench+  ( benchSuite+  ) where++import Gauge.Main+#if __GLASGOW_HASKELL__ < 803+import Data.Semigroup+#endif+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as U+#if __GLASGOW_HASKELL__ < 709+import Data.Word+#endif++import Math.NumberTheory.ArithmeticFunctions.Moebius+import Math.NumberTheory.ArithmeticFunctions.SieveBlock+import Math.NumberTheory.Primes.Factorisation (totientSieve, sieveTotient, carmichaelSieve, sieveCarmichael)++blockLen :: Word+blockLen = 10^6++totientHelper :: Word -> Word -> Word+totientHelper p 1 = p - 1+totientHelper p 2 = (p - 1) * p+totientHelper p k = (p - 1) * p ^ (k - 1)++totientBlockConfig :: SieveBlockConfig Word+totientBlockConfig = SieveBlockConfig+  { sbcEmpty                = 1+  , sbcAppend               = (*)+  , sbcFunctionOnPrimePower = totientHelper+  }++sumOldTotientSieve :: Word -> Word+sumOldTotientSieve len' = sum $ map (fromInteger . sieveTotient sieve) [1 .. len]+  where+    len = toInteger len'+    sieve = totientSieve len++carmichaelHelper :: Word -> Word -> Word+carmichaelHelper 2 1 = 1+carmichaelHelper 2 2 = 2+carmichaelHelper 2 k = 2 ^ (k - 2)+carmichaelHelper p 1 = p - 1+carmichaelHelper p 2 = (p - 1) * p+carmichaelHelper p k = (p - 1) * p ^ (k - 1)++carmichaelBlockConfig :: SieveBlockConfig Word+carmichaelBlockConfig = SieveBlockConfig+  { sbcEmpty                = 1+  -- There is a specialized 'gcd' for Word, but not 'lcm'.+  , sbcAppend               = (\x y -> (x `quot` (gcd x y)) * y)+  , sbcFunctionOnPrimePower = carmichaelHelper+  }++sumOldCarmichaelSieve :: Word -> Word+sumOldCarmichaelSieve len' = sum $ map (fromInteger . sieveCarmichael sieve) [1 .. len]+  where+    len = toInteger len'+    sieve = carmichaelSieve len++moebiusConfig :: SieveBlockConfig Moebius+moebiusConfig = SieveBlockConfig+  { sbcEmpty                = MoebiusP+  , sbcAppend               = (<>)+  , sbcFunctionOnPrimePower = const $ \case+      0 -> MoebiusP+      1 -> MoebiusN+      _ -> MoebiusZ+  }++benchSuite :: Benchmark+benchSuite = bgroup "SieveBlock"+  [ bgroup "totient"+    [ bench "old"     $ nf sumOldTotientSieve blockLen+    , bench "boxed"   $ nf (V.sum . sieveBlock        totientBlockConfig 1) blockLen+    , bench "unboxed" $ nf (U.sum . sieveBlockUnboxed totientBlockConfig 1) blockLen+    ]+  , bgroup "carmichael"+    [ bench "old"     $ nf sumOldCarmichaelSieve blockLen+    , bench "boxed"   $ nf (V.sum . sieveBlock        carmichaelBlockConfig 1) blockLen+    , bench "unboxed" $ nf (U.sum . sieveBlockUnboxed carmichaelBlockConfig 1) blockLen+    ]+  , bgroup "moebius"+    [ bench "boxed"   $ nf (V.sum . V.map runMoebius . sieveBlock        moebiusConfig 1 :: Word -> Int) blockLen+    , bench "unboxed" $ nf (U.sum . U.map runMoebius . sieveBlockUnboxed moebiusConfig 1 :: Word -> Int) blockLen+    , bench "special" $ nf (U.sum . U.map runMoebius . sieveBlockMoebius 1 :: Word -> Int) blockLen+    ]+  ]

@@ -0,0 +1,76 @@+-- |+-- Module:      Math.NumberTheory.ArithmeticFunctions.MertensTests+-- Copyright:   (c) 2018 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.ArithmeticFunctions.Mertens+--++{-# LANGUAGE CPP        #-}+{-# LANGUAGE LambdaCase #-}++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.ArithmeticFunctions.MertensTests+  ( testSuite+  ) where++import Test.Tasty++#if __GLASGOW_HASKELL__ < 803+import Data.Semigroup+#endif+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as U+#if __GLASGOW_HASKELL__ < 709+import Data.Word+#endif++import Math.NumberTheory.ArithmeticFunctions+import Math.NumberTheory.ArithmeticFunctions.Mertens+import Math.NumberTheory.ArithmeticFunctions.SieveBlock+import Math.NumberTheory.TestUtils++moebiusConfig :: SieveBlockConfig Moebius+moebiusConfig = SieveBlockConfig+  { sbcEmpty                = MoebiusP+  , sbcAppend               = (<>)+  , sbcFunctionOnPrimePower = const $ \case+      0 -> MoebiusP+      1 -> MoebiusN+      _ -> MoebiusZ+  }++mertensDiffPointwise :: Word -> Word -> Int+mertensDiffPointwise lo len = sum $ map (runMoebius . moebius) [lo + 1 .. lo + len]++mertensDiffBlockSpecial :: Word -> Word -> Int+mertensDiffBlockSpecial lo len = U.sum $ U.map runMoebius+  $ sieveBlockMoebius (lo + 1) len++mertensDiffBlockUnboxed :: Word -> Word -> Int+mertensDiffBlockUnboxed lo len = U.sum $ U.map runMoebius+  $ sieveBlockUnboxed moebiusConfig (lo + 1) len++mertensDiffBlockBoxed :: Word -> Word -> Int+mertensDiffBlockBoxed lo len = V.sum $ V.map runMoebius+  $ sieveBlock moebiusConfig (lo + 1) len++mertensDiff :: Word -> Word -> Int+mertensDiff lo len = mertens (lo + len) - mertens lo++propertyCompare :: (Word -> Word -> Int) -> Word -> Word -> Bool+propertyCompare func lo' len' = mertensDiff lo len == func lo len+  where+    lo  = lo'  `rem` 10000000+    len = len' `rem` 1000++testSuite :: TestTree+testSuite = testGroup "Mertens"+  [ testSmallAndQuick "pointwise"     $ propertyCompare mertensDiffPointwise+  , testSmallAndQuick "block special" $ propertyCompare mertensDiffBlockSpecial+  , testSmallAndQuick "block unboxed" $ propertyCompare mertensDiffBlockUnboxed+  , testSmallAndQuick "block boxed"   $ propertyCompare mertensDiffBlockBoxed+  ]

@@ -0,0 +1,109 @@+-- |+-- Module:      Math.NumberTheory.ArithmeticFunctions.SieveBlockTests+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.ArithmeticFunctions.SieveBlock+--++{-# LANGUAGE CPP        #-}+{-# LANGUAGE LambdaCase #-}++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.ArithmeticFunctions.SieveBlockTests+  ( testSuite+  ) where++import Test.Tasty+import Test.Tasty.HUnit++#if __GLASGOW_HASKELL__ < 803+import Data.Semigroup+#endif+import qualified Data.Vector as V+import qualified Data.Vector.Unboxed as U+#if __GLASGOW_HASKELL__ < 709+import Data.Word+#endif++import Math.NumberTheory.ArithmeticFunctions+import Math.NumberTheory.ArithmeticFunctions.SieveBlock++pointwiseTest :: (Eq a, Show a) => ArithmeticFunction Word a -> Word -> Word -> IO ()+pointwiseTest f lowIndex len = assertEqual "pointwise"+    (runFunctionOverBlock f lowIndex len)+    (V.generate (fromIntegral len) (runFunction f . (+ lowIndex) . fromIntegral))++unboxedTest :: (Eq a, U.Unbox a, Show a) => SieveBlockConfig a -> IO ()+unboxedTest config = assertEqual "unboxed"+    (sieveBlock config 1 1000)+    (U.convert $ sieveBlockUnboxed config 1 1000)++moebiusTest :: Word -> Word -> Bool+moebiusTest m n+  = m == 0+  || sieveBlockUnboxed moebiusConfig m n+  == sieveBlockMoebius m n++moebiusSpecialCases :: [TestTree]+moebiusSpecialCases = map (uncurry pairToTest)+  [ (1, 1)+  , (1, 2)+  , (208, 298)+  , (1, 12835)+  , (10956, 4430)+  , (65, 16171)+  , (120906, 19456)+  , (33800000, 27002)+  , (17266222643, 5051)+  , (1000158, 48758)+  , (1307265, 3725)+  , (2600000, 14686)+  , (4516141422507 - 100000, 100001)+  , (1133551497049257 - 100000, 100001)+  -- too long for regular runs+  -- , (1157562178759482171 - 100000, 100001)+  ]+  where+    pairToTest :: Word -> Word -> TestTree+    pairToTest m n = testCase (show m ++ "," ++ show n) $ assertBool "should be equal" $ moebiusTest m n++multiplicativeConfig :: (Word -> Word -> Word) -> SieveBlockConfig Word+multiplicativeConfig f = SieveBlockConfig+  { sbcEmpty                = 1+  , sbcAppend               = (*)+  , sbcFunctionOnPrimePower = f+  }++moebiusConfig :: SieveBlockConfig Moebius+moebiusConfig = SieveBlockConfig+  { sbcEmpty = MoebiusP+  , sbcAppend = (<>)+  , sbcFunctionOnPrimePower = const $ \case+      0 -> MoebiusP+      1 -> MoebiusN+      _ -> MoebiusZ+  }++testSuite :: TestTree+testSuite = testGroup "SieveBlock"+  [ testGroup "pointwise"+    [ testCase "divisors"   $ pointwiseTest divisorsA   1 1000+    , testCase "tau"        $ pointwiseTest tauA        1 1000+    , testCase "totient"    $ pointwiseTest totientA    1 1000+    , testCase "moebius"    $ pointwiseTest moebiusA    1 1000+    , testCase "smallOmega" $ pointwiseTest smallOmegaA 1 1000+    , testCase "bigOmega"   $ pointwiseTest bigOmegaA   1 1000+    , testCase "carmichael" $ pointwiseTest carmichaelA 1 1000+    ]+  , testGroup "unboxed"+    [ testCase "id"      $ unboxedTest $ multiplicativeConfig (^)+    , testCase "tau"     $ unboxedTest $ multiplicativeConfig (const id)+    , testCase "moebius" $ unboxedTest moebiusConfig+    , testCase "totient" $ unboxedTest $ multiplicativeConfig (\p a -> (p - 1) * p ^ (a - 1))+    ]+  , testGroup "special moebius" moebiusSpecialCases+  ]

@@ -10,6 +10,7 @@  {-# LANGUAGE CPP       #-} +{-# OPTIONS_GHC -fno-warn-deprecations  #-} {-# OPTIONS_GHC -fno-warn-type-defaults #-}  module Math.NumberTheory.ArithmeticFunctionsTests@@ -137,21 +138,17 @@   , 1728, 1584, 2208, 1536   ] --- | moebius values are [-1, 0, 1]-moebiusProperty1 :: Natural -> Bool-moebiusProperty1 n = runFunction moebiusA n `elem` [-1, 0, 1]- -- | moebius does not require full factorisation moebiusLazy :: Assertion-moebiusLazy = assertEqual "moebius" 0 (runFunction moebiusA (2^2 * (2^100000-1) :: Natural))+moebiusLazy = assertEqual "moebius" MoebiusZ (runFunction moebiusA (2^2 * (2^100000-1) :: Natural))  -- | moebius matches baseline from OEIS. moebiusOeis :: Assertion moebiusOeis = oeisAssertion "A008683" moebiusA-  [ 1, -1, -1, 0, -1, 1, -1, 0, 0, 1, -1, 0, -1, 1, 1, 0, -1, 0, -1, 0, 1, 1, -1-  , 0, 0, 1, 0, 0, -1, -1, -1, 0, 1, 1, 1, 0, -1, 1, 1, 0, -1, -1, -1, 0, 0, 1-  , -1, 0, 0, 0, 1, 0, -1, 0, 1, 0, 1, 1, -1, 0, -1, 1, 0, 0, 1, -1, -1, 0, 1-  , -1, -1, 0, -1, 1, 0, 0, 1+  [ MoebiusP, MoebiusN, MoebiusN, MoebiusZ, MoebiusN, MoebiusP, MoebiusN, MoebiusZ, MoebiusZ, MoebiusP, MoebiusN, MoebiusZ, MoebiusN, MoebiusP, MoebiusP, MoebiusZ, MoebiusN, MoebiusZ, MoebiusN, MoebiusZ, MoebiusP, MoebiusP, MoebiusN+  , MoebiusZ, MoebiusZ, MoebiusP, MoebiusZ, MoebiusZ, MoebiusN, MoebiusN, MoebiusN, MoebiusZ, MoebiusP, MoebiusP, MoebiusP, MoebiusZ, MoebiusN, MoebiusP, MoebiusP, MoebiusZ, MoebiusN, MoebiusN, MoebiusN, MoebiusZ, MoebiusZ, MoebiusP+  , MoebiusN, MoebiusZ, MoebiusZ, MoebiusZ, MoebiusP, MoebiusZ, MoebiusN, MoebiusZ, MoebiusP, MoebiusZ, MoebiusP, MoebiusP, MoebiusN, MoebiusZ, MoebiusN, MoebiusP, MoebiusZ, MoebiusZ, MoebiusP, MoebiusN, MoebiusN, MoebiusZ, MoebiusP+  , MoebiusN, MoebiusN, MoebiusZ, MoebiusN, MoebiusP, MoebiusZ, MoebiusZ, MoebiusP   ]  -- | liouville values are [-1, 1]@@ -160,7 +157,7 @@  -- | moebius is zero or equal to liouville liouvilleProperty2 :: Natural -> Bool-liouvilleProperty2 n = m == 0 || l == m+liouvilleProperty2 n = m == MoebiusZ || l == runMoebius m   where     l = runFunction liouvilleA n     m = runFunction moebiusA   n@@ -260,8 +257,7 @@     , testCase          "OEIS jordan_2"      jordan2Oeis     ]   , testGroup "Moebius"-    [ testSmallAndQuick "moebius values" moebiusProperty1-    , testCase          "OEIS"           moebiusOeis+    [ testCase          "OEIS"           moebiusOeis     , testCase          "Lazy"           moebiusLazy     ]   , testGroup "Liouville"

@@ -8,6 +8,7 @@ -- Tests for Math.NumberTheory.GCD -- +{-# LANGUAGE CPP                 #-} {-# LANGUAGE ScopedTypeVariables #-}  {-# OPTIONS_GHC -fno-warn-type-defaults #-}@@ -17,8 +18,16 @@   ) where  import Test.Tasty+import Test.Tasty.HUnit +import Control.Arrow import Data.Bits+import Data.List (tails)+#if MIN_VERSION_base(4,8,0)+#else+import Data.Word+#endif+import Numeric.Natural  import Math.NumberTheory.GCD import Math.NumberTheory.TestUtils@@ -42,9 +51,53 @@ coprimeProperty :: (Integral a, Bits a) => AnySign a -> AnySign a -> Bool coprimeProperty (AnySign a) (AnySign b) = coprime a b == (gcd a b == 1) +splitIntoCoprimesProperty1 :: [(Positive Natural, Power Word)] -> Bool+splitIntoCoprimesProperty1 fs' = factorback fs == factorback (splitIntoCoprimes fs)+  where+    fs = map (getPositive *** getPower) fs'+    factorback = product . map (uncurry (^))++splitIntoCoprimesProperty2 :: [(Positive Natural, Power Word)] -> Bool+splitIntoCoprimesProperty2 fs' = multiplicities fs <= multiplicities (splitIntoCoprimes fs)+  where+    fs = map (getPositive *** getPower) fs'+    multiplicities = sum . map snd . filter ((/= 1) . fst)++splitIntoCoprimesProperty3 :: [(Positive Natural, Power Word)] -> Bool+splitIntoCoprimesProperty3 fs' = and [ coprime x y | (x : xs) <- tails fs, y <- xs ]+  where+    fs = map fst $ splitIntoCoprimes $ map (getPositive *** getPower) fs'++-- | Check that evaluation never freezes.+splitIntoCoprimesProperty4 :: [(Integer, Word)] -> Bool+splitIntoCoprimesProperty4 fs' = fs == fs+  where+    fs = splitIntoCoprimes fs'++-- | This is an undefined behaviour, but at least it should not+-- throw exceptions or loop forever.+splitIntoCoprimesSpecialCase1 :: Assertion+splitIntoCoprimesSpecialCase1 =+  assertBool "should not fail" $ splitIntoCoprimesProperty4 [(0, 0), (0, 0)]++-- | This is an undefined behaviour, but at least it should not+-- throw exceptions or loop forever.+splitIntoCoprimesSpecialCase2 :: Assertion+splitIntoCoprimesSpecialCase2 =+  assertBool "should not fail" $ splitIntoCoprimesProperty4 [(0, 1), (-2, 0)]+ testSuite :: TestTree testSuite = testGroup "GCD"   [ testSameIntegralProperty "binaryGCD"   binaryGCDProperty   , testSameIntegralProperty "extendedGCD" extendedGCDProperty   , testSameIntegralProperty "coprime"     coprimeProperty+  , testGroup "splitIntoCoprimes"+    [ testSmallAndQuick "preserves product of factors"        splitIntoCoprimesProperty1+    , testSmallAndQuick "number of factors is non-decreasing" splitIntoCoprimesProperty2+    , testSmallAndQuick "output factors are coprime"          splitIntoCoprimesProperty3++    , testCase          "does not freeze 1"                   splitIntoCoprimesSpecialCase1+    , testCase          "does not freeze 2"                   splitIntoCoprimesSpecialCase2+    , testSmallAndQuick "does not freeze random"              splitIntoCoprimesProperty4+    ]   ]

@@ -2,7 +2,7 @@  -- | -- Module:      Math.NumberTheory.GaussianIntegersTests--- Copyright:   (c) 2016 Chris Fredrickson+-- Copyright:   (c) 2016 Chris Fredrickson, Google Inc. -- Licence:     MIT -- Maintainer:  Chris Fredrickson <chris.p.fredrickson@gmail.com> -- Stability:   Provisional@@ -14,15 +14,25 @@   ( testSuite   ) where +import Control.Monad (zipWithM_) import Test.Tasty import Test.Tasty.HUnit  import Math.NumberTheory.GaussianIntegers import Math.NumberTheory.TestUtils +lazyCases :: [(GaussianInteger, [(GaussianInteger, Int)])]+lazyCases =+  [ ( 14145130733+    * 10000000000000000000000000000000000000121+    * 100000000000000000000000000000000000000000000000447+    , [(21037 :+ 117058, 1), (117058 :+ 21037, 1)]+    )+  ]+ -- | Number is zero or is equal to the product of its factors.-factoriseProperty :: Integer -> Integer -> Bool-factoriseProperty x y+factoriseProperty1 :: Integer -> Integer -> Bool+factoriseProperty1 x y   =  x == 0 && y == 0   || g == g'   where@@ -30,6 +40,9 @@     factors = factorise g     g' = product $ map (uncurry (.^)) factors +factoriseProperty2 :: (GaussianInteger, [(GaussianInteger, Int)]) -> Assertion+factoriseProperty2 (n, fs) = zipWithM_ (assertEqual (show n)) fs (factorise n)+ -- | Number is prime iff it is non-zero --   and has exactly one (non-unit) factor. isPrimeProperty :: Integer -> Integer -> Bool@@ -69,11 +82,13 @@ gcdGSpecialCase1 = assertEqual "gcdG" 1 $ gcdG (12 :+ 23) (23 :+ 34)  testSuite :: TestTree-testSuite = testGroup "GaussianIntegers"-  [ testSmallAndQuick "factorise"         factoriseProperty+testSuite = testGroup "GaussianIntegers" $+  [ testSmallAndQuick "factorise"         factoriseProperty1   , testSmallAndQuick "isPrime"           isPrimeProperty   , testSmallAndQuick "primes"            primesGeneratesPrimesProperty   , testSmallAndQuick "signumAbsProperty" signumAbsProperty   , testSmallAndQuick "absProperty"       absProperty   , testCase          "gcdG (12 :+ 23) (23 :+ 34)" gcdGSpecialCase1   ]+  +++  map (\x -> testCase ("laziness " ++ show (fst x)) (factoriseProperty2 x)) lazyCases

@@ -20,7 +20,9 @@ import Test.Tasty  import Data.Bits+#if __GLASGOW_HASKELL__ < 803 import Data.Semigroup+#endif  import Math.NumberTheory.Moduli hiding (invertMod) import Math.NumberTheory.TestUtils

@@ -0,0 +1,118 @@+-- |+-- Module:      Math.NumberTheory.Moduli.PrimitiveRootTests+-- Copyright:   (c) 2017 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.Moduli.PrimitiveRoot+--++{-# LANGUAGE CPP                 #-}+{-# LANGUAGE ScopedTypeVariables #-}++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.Moduli.PrimitiveRootTests+  ( testSuite+  ) where++import Test.Tasty++import qualified Data.Set as S+import Data.List (genericTake)+import Data.Maybe+import Numeric.Natural++#if !(MIN_VERSION_base(4,8,0))+import Data.Word+#endif++import Math.NumberTheory.ArithmeticFunctions (totient)+import Math.NumberTheory.Moduli.Class (Mod, SomeMod(..), modulo)+import Math.NumberTheory.Moduli.PrimitiveRoot+import Math.NumberTheory.Prefactored (prefValue)+import Math.NumberTheory.TestUtils+import Math.NumberTheory.UniqueFactorisation++cyclicGroupProperty1 :: (Integral a, UniqueFactorisation a, Show a) => AnySign a -> Bool+cyclicGroupProperty1 (AnySign n) = case cyclicGroupFromModulo n of+  Nothing -> True+  Just cg -> prefValue (cyclicGroupToModulo cg) == n++-- | Multiplicative groups modulo primes are always cyclic.+cyclicGroupProperty2 :: (Integral a, UniqueFactorisation a) => Positive a -> Bool+cyclicGroupProperty2 (Positive n) = case isPrime n of+  Nothing -> True+  Just _  -> isJust (cyclicGroupFromModulo n)++-- | Multiplicative groups modulo double primes are always cyclic.+cyclicGroupProperty3 :: (Integral a, UniqueFactorisation a) => Positive a -> Bool+cyclicGroupProperty3 (Positive n) = case isPrime n of+  Nothing -> True+  Just _  -> 2 * n < n {- overflow check -}+          || isJust (cyclicGroupFromModulo n)++allUnique :: Ord a => [a] -> Bool+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++isPrimitiveRoot'Property1+  :: (Eq a, Integral a, UniqueFactorisation a)+  => AnySign a -> CyclicGroup a -> Bool+isPrimitiveRoot'Property1 (AnySign n) cg+  = gcd n (prefValue (cyclicGroupToModulo cg)) == 1+  || not (isPrimitiveRoot' cg n)++isPrimitiveRootProperty1 :: AnySign Integer -> Positive Natural -> Bool+isPrimitiveRootProperty1 (AnySign n) (Positive m)+  = case n `modulo` m of+    SomeMod n' -> gcd n (toInteger m) == 1+               || not (isPrimitiveRoot n')+    InfMod{}   -> False++isPrimitiveRootProperty2 :: Positive Natural -> Bool+isPrimitiveRootProperty2 (Positive m)+  = isNothing (cyclicGroupFromModulo m)+  || case 0 `modulo` m of+    SomeMod (_ :: Mod t) -> any isPrimitiveRoot [(minBound :: Mod t) .. maxBound]+    InfMod{}             -> False++isPrimitiveRootProperty3 :: AnySign Integer -> Positive Natural -> Bool+isPrimitiveRootProperty3 (AnySign n) (Positive m)+  = case n `modulo` m of+    SomeMod n' -> not (isPrimitiveRoot n')+               || allUnique (genericTake (totient m - 1) (iterate (* n') 1))+    InfMod{}   -> False++isPrimitiveRootProperty4 :: AnySign Integer -> Positive Natural -> Bool+isPrimitiveRootProperty4 (AnySign n) (Positive m)+  = isJust (cyclicGroupFromModulo m)+  || case n `modulo` m of+    SomeMod n' -> not (isPrimitiveRoot n')+    InfMod{}   -> False++testSuite :: TestTree+testSuite = testGroup "Primitive root"+  [ testGroup "CyclicGroup"+    [ testIntegralProperty "cyclicGroupToModulo . cyclicGroupFromModulo" cyclicGroupProperty1+    , testIntegralProperty "cyclic group mod p" cyclicGroupProperty2+    , testIntegralProperty "cyclic group mod 2p" cyclicGroupProperty3+    ]+  , testGroup "isPrimitiveRoot'"+    [ testGroup "primitive root is coprime with modulo"+      [ testSmallAndQuick "Integer" (isPrimitiveRoot'Property1 :: AnySign Integer -> CyclicGroup Integer -> Bool)+      , testSmallAndQuick "Natural" (isPrimitiveRoot'Property1 :: AnySign Natural -> CyclicGroup Natural -> Bool)+      , testSmallAndQuick "Int"     (isPrimitiveRoot'Property1 :: AnySign Int     -> CyclicGroup Int     -> Bool)+      , testSmallAndQuick "Word"    (isPrimitiveRoot'Property1 :: AnySign Word    -> CyclicGroup Word    -> Bool)+      ]+    ]+  , testGroup "isPrimitiveRoot"+    [ testSmallAndQuick "primitive root is coprime with modulo" isPrimitiveRootProperty1+    , testSmallAndQuick "cyclic group has a primitive root"     isPrimitiveRootProperty2+    , testSmallAndQuick "primitive root generates cyclic group" isPrimitiveRootProperty3+    , testSmallAndQuick "no primitive root in non-cyclic group" isPrimitiveRootProperty4+    ]+  ]

@@ -21,35 +21,40 @@ import Test.Tasty.HUnit  import Data.List (nub)+import Data.Maybe (fromJust)  import Math.NumberTheory.Moduli hiding (invertMod)+import Math.NumberTheory.UniqueFactorisation (unPrime, isPrime) import Math.NumberTheory.TestUtils -unwrapPP :: (Prime, Power Int) -> (Integer, Int)-unwrapPP (Prime p, Power e) = (p, e)+unwrapP :: PrimeWrapper Integer -> Integer+unwrapP (PrimeWrapper p) = unPrime p +unwrapPP :: (PrimeWrapper Integer, Power Int) -> (Integer, Int)+unwrapPP (p, Power e) = (unwrapP p, e)+ -- | Check that 'sqrtMod' is defined iff a quadratic residue exists. --   Also check that the result is a solution of input modular equation.-sqrtModPProperty :: AnySign Integer -> Prime -> Bool-sqrtModPProperty (AnySign n) (Prime p) = case sqrtModP n p of+sqrtModPProperty :: AnySign Integer -> PrimeWrapper Integer -> Bool+sqrtModPProperty (AnySign n) (unwrapP -> p) = case sqrtModP n p of   Nothing -> jacobi n p == MinusOne   Just rt -> (p == 2 || jacobi n p /= MinusOne) && rt ^ 2 `mod` p == n `mod` p -sqrtModPListProperty :: AnySign Integer -> Prime -> Bool-sqrtModPListProperty (AnySign n) (Prime p) = all (\rt -> rt ^ 2 `mod` p == n `mod` p) (sqrtModPList n p)+sqrtModPListProperty :: AnySign Integer -> PrimeWrapper Integer -> Bool+sqrtModPListProperty (AnySign n) (unwrapP -> p) = all (\rt -> rt ^ 2 `mod` p == n `mod` p) (sqrtModPList n p) -sqrtModP'Property :: Positive Integer -> Prime -> Bool-sqrtModP'Property (Positive n) (Prime p) = (p /= 2 && jacobi n p /= One) || rt ^ 2 `mod` p == n `mod` p+sqrtModP'Property :: Positive Integer -> PrimeWrapper Integer -> Bool+sqrtModP'Property (Positive n) (unwrapP -> p) = (p /= 2 && jacobi n p /= One) || rt ^ 2 `mod` p == n `mod` p   where     rt = sqrtModP' n p -tonelliShanksProperty1 :: Positive Integer -> Prime -> Bool-tonelliShanksProperty1 (Positive n) (Prime p) = p `mod` 4 /= 1 || jacobi n p /= One || rt ^ 2 `mod` p == n `mod` p+tonelliShanksProperty1 :: Positive Integer -> PrimeWrapper Integer -> Bool+tonelliShanksProperty1 (Positive n) (unwrapP -> p) = p `mod` 4 /= 1 || jacobi n p /= One || rt ^ 2 `mod` p == n `mod` p   where     rt = tonelliShanks n p -tonelliShanksProperty2 :: Prime -> Bool-tonelliShanksProperty2 (Prime p) = p `mod` 4 /= 1 || rt ^ 2 `mod` p == n `mod` p+tonelliShanksProperty2 :: PrimeWrapper Integer -> Bool+tonelliShanksProperty2 (unwrapP -> p) = p `mod` 4 /= 1 || rt ^ 2 `mod` p == n `mod` p   where     n  = head $ filter (\s -> jacobi s p == One) [2..p-1]     rt = tonelliShanks n p@@ -61,13 +66,13 @@     ps = [17, 73, 241, 1009, 2689, 8089, 33049, 53881, 87481, 483289, 515761, 1083289, 3818929, 9257329, 22000801, 48473881, 175244281, 427733329, 898716289, 8114538721, 9176747449, 23616331489]     rts = map (\p -> tonelliShanks 2 p) ps -sqrtModPPProperty :: AnySign Integer -> (Prime, Power Int) -> Bool-sqrtModPPProperty (AnySign n) (Prime p, Power e) = gcd n p > 1 || case sqrtModPP n (p, e) of+sqrtModPPProperty :: AnySign Integer -> (PrimeWrapper Integer, Power Int) -> Bool+sqrtModPPProperty (AnySign n) (unwrapP -> p, Power e) = gcd n p > 1 || case sqrtModPP n (p, e) of   Nothing -> True   Just rt -> rt ^ 2 `mod` (p ^ e) == n `mod` (p ^ e)  sqrtModPPBase2Property :: AnySign Integer -> Power Int -> Bool-sqrtModPPBase2Property n e = sqrtModPPProperty n (Prime 2, e)+sqrtModPPBase2Property n e = sqrtModPPProperty n (PrimeWrapper $ fromJust $ isPrime (2 :: Integer), e)  sqrtModPPSpecialCase1 :: Assertion sqrtModPPSpecialCase1 =@@ -77,16 +82,16 @@ sqrtModPPSpecialCase2 =   assertEqual "sqrtModPP 16 3 2 = 4" (Just 4) (sqrtModPP 16 (3, 2)) -sqrtModPPListProperty :: AnySign Integer -> (Prime, Power Int) -> Bool-sqrtModPPListProperty (AnySign n) (Prime p, Power e) = gcd n p > 1+sqrtModPPListProperty :: AnySign Integer -> (PrimeWrapper Integer, Power Int) -> Bool+sqrtModPPListProperty (AnySign n) (unwrapP -> p, Power e) = gcd n p > 1   || all (\rt -> rt ^ 2 `mod` (p ^ e) == n `mod` (p ^ e)) (sqrtModPPList n (p, e)) -sqrtModFProperty :: AnySign Integer -> [(Prime, Power Int)] -> Bool+sqrtModFProperty :: AnySign Integer -> [(PrimeWrapper Integer, Power Int)] -> Bool sqrtModFProperty (AnySign n) (map unwrapPP -> pes) = case sqrtModF n pes of   Nothing -> True   Just rt -> all (\(p, e) -> rt ^ 2 `mod` (p ^ e) == n `mod` (p ^ e)) pes -sqrtModFListProperty :: AnySign Integer -> [(Prime, Power Int)] -> Bool+sqrtModFListProperty :: AnySign Integer -> [(PrimeWrapper Integer, Power Int)] -> Bool sqrtModFListProperty (AnySign n) (map unwrapPP -> pes)   = nub ps /= ps || all     (\rt -> all (\(p, e) -> rt ^ 2 `mod` (p ^ e) == n `mod` (p ^ e)) pes)

@@ -37,7 +37,7 @@ generalInversionProperty :: (Int -> Int) -> Positive Int -> Bool generalInversionProperty g (Positive n)   =  g n == sum [f (n `quot` k) | k <- [1 .. n]]-  && f n == sum [fromInteger (moebius (toInteger k)) * g (n `quot` k) | k <- [1 .. n]]+  && f n == sum [runMoebius (moebius k) * g (n `quot` k) | k <- [1 .. n]]   where     f = generalInversion g 

@@ -37,7 +37,7 @@ generalInversionProperty :: (Int -> Integer) -> Positive Int -> Bool generalInversionProperty g (Positive n)   =  g n == sum [f (n `quot` k) | k <- [1 .. n]]-  && f n == sum [moebius (toInteger k) * g (n `quot` k) | k <- [1 .. n]]+  && f n == sum [runMoebius (moebius k) * g (n `quot` k) | k <- [1 .. n]]   where     f = generalInversion g 

@@ -0,0 +1,98 @@+-- |+-- Module:      Math.NumberTheory.Powers.ModularTests+-- Copyright:   (c) 2017 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.Powers.Modular+--++{-# LANGUAGE CPP #-}++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.Powers.ModularTests+  ( testSuite+  ) where++import Test.Tasty+#if MIN_VERSION_base(4,8,0)+import Test.Tasty.HUnit+#endif++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++#if __GLASGOW_HASKELL__ > 709+-- | 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))+#endif++-- | 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' (fromIntegral b) (fromIntegral e) (fromIntegral 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) (fromIntegral e) (fromIntegral m))++#if WORD_SIZE_IN_BITS == 64 && __GLASGOW_HASKELL__ > 709+-- | 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++#if __GLASGOW_HASKELL__ > 709+    , testSmallAndQuick "powModInt"              powModProperty_Int+    , testSmallAndQuick "powModWord"             powModProperty_Word+#endif+    , testSmallAndQuick "powModInteger"          powModProperty_Integer+    , testSmallAndQuick "powModNatural"          powModProperty_Natural++#if WORD_SIZE_IN_BITS == 64 && __GLASGOW_HASKELL__ > 709+    , testCase          "large modulo :: Int"    powModSpecialCase1_Int+    , testCase          "large modulo :: Word"   powModSpecialCase1_Word+#endif+    ]+  ]

@@ -0,0 +1,101 @@+-- |+-- Module:      Math.NumberTheory.PrefactoredTests+-- Copyright:   (c) 2017 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.Prefactored+--++{-# LANGUAGE CPP #-}++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.PrefactoredTests+  ( testSuite+  ) where++import Test.Tasty++import Control.Arrow (second)+import Data.Bits (Bits)+import Data.List (tails)+#if MIN_VERSION_base(4,8,0)+#else+import Data.Word+#endif+import Numeric.Natural++import Math.NumberTheory.GCD (coprime, splitIntoCoprimes)+import Math.NumberTheory.Prefactored+import Math.NumberTheory.TestUtils++isValid :: (Eq a, Bits a, Integral a) => Prefactored a -> Bool+isValid pref+  = abs n == abs (product (map (uncurry (^)) fs))+  && and [ coprime g h | ((g, _) : gs) <- tails fs, (h, _) <- gs ]+  where+    n  = prefValue   pref+    fs = prefFactors pref++fromValueProperty :: Integer -> Bool+fromValueProperty n = isValid pref && prefValue pref == n+  where+    pref = fromValue n++fromFactorsProperty :: [(Integer, Power Word)] -> Bool+fromFactorsProperty fs' = isValid pref && abs (prefValue pref) == abs (product (map (uncurry (^)) fs))+  where+    fs   = map (second getPower) fs'+    pref = fromFactors (splitIntoCoprimes fs)++plusProperty :: Integer -> Integer -> Bool+plusProperty x y = isValid z && prefValue z == x + y+  where+    z = fromValue x + fromValue y++minusProperty :: Integer -> Integer -> Bool+minusProperty x y = isValid z && prefValue z == x - y+  where+    z = fromValue x - fromValue y++minusNaturalProperty :: Natural -> Natural -> Bool+minusNaturalProperty x y = x < y || (isValid z && prefValue z == x - y)+  where+    z = fromValue x - fromValue y++multiplyProperty :: Integer -> Integer -> Bool+multiplyProperty x y = isValid z && prefValue z == x * y+  where+    z = fromValue x * fromValue y++negateProperty :: Integer -> Bool+negateProperty x = isValid z && prefValue z == negate x+  where+    z = negate (fromValue x)++absSignumProperty :: Integer -> Bool+absSignumProperty x = isValid z && prefValue z == x+  where+    z = abs (fromValue x) * signum (fromValue x)++fromIntegerProperty :: Integer -> Bool+fromIntegerProperty n = isValid pref && prefValue pref == n+  where+    pref = fromInteger n++testSuite :: TestTree+testSuite = testGroup "Prefactored"+  [ testSmallAndQuick "fromValue"   fromValueProperty+  , testSmallAndQuick "fromFactors" fromFactorsProperty+  , testGroup "Num instance"+    [ testSmallAndQuick "plus"         plusProperty+    , testSmallAndQuick "minus"        minusProperty+    , testSmallAndQuick "minusNatural" minusNaturalProperty+    , testSmallAndQuick "multiply"     multiplyProperty+    , testSmallAndQuick "negate"       negateProperty+    , testSmallAndQuick "absSignum"    absSignumProperty+    , testSmallAndQuick "fromInteger"  fromIntegerProperty+    ]+  ]

@@ -17,6 +17,7 @@ import Test.Tasty import Test.Tasty.HUnit +import Control.Monad (zipWithM_) import Data.List (nub, sort)  import Math.NumberTheory.Primes.Factorisation@@ -46,6 +47,15 @@   , (16757651897802863152387219654541878166,[(2,1),(23,1),(277,1),(505353699591289,1),(2602436338718275457,1)])   ] +lazyCases :: [(Integer, [(Integer, Int)])]+lazyCases =+  [ ( 14145130711+    * 10000000000000000000000000000000000000121+    * 100000000000000000000000000000000000000000000000447+    , [(14145130711, 1)]+    )+  ]+ factoriseProperty1 :: Assertion factoriseProperty1 = assertEqual "0" [] (factorise 1) @@ -64,8 +74,11 @@ factoriseProperty5 (Positive n) = product (map (uncurry (^)) (factorise n)) == n  factoriseProperty6 :: (Integer, [(Integer, Int)]) -> Assertion-factoriseProperty6 (n, fs) = assertEqual (show n) fs (factorise n)+factoriseProperty6 (n, fs) = assertEqual (show n) (sort fs) (sort (factorise n)) +factoriseProperty7 :: (Integer, [(Integer, Int)]) -> Assertion+factoriseProperty7 (n, fs) = zipWithM_ (assertEqual (show n)) fs (factorise n)+ testSuite :: TestTree testSuite = testGroup "Factorisation"   [ testGroup "factorise" $@@ -76,4 +89,6 @@     , testSmallAndQuick "factorback"                     factoriseProperty5     ] ++     map (\x -> testCase ("special case " ++ show (fst x)) (factoriseProperty6 x)) specialCases+    +++    map (\x -> testCase ("laziness " ++ show (fst x)) (factoriseProperty7 x)) lazyCases   ]

@@ -1,67 +0,0 @@--- |--- Module:      Math.NumberTheory.Primes.HeapTests--- Copyright:   (c) 2016 Andrew Lelechenko--- Licence:     MIT--- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>--- Stability:   Provisional------ Tests for Math.NumberTheory.Primes.Heap-----{-# LANGUAGE CPP #-}--{-# OPTIONS_GHC -fno-warn-type-defaults #-}--module Math.NumberTheory.Primes.HeapTests-  ( testSuite-  ) where--import Prelude hiding (words)--import Test.Tasty-import Test.Tasty.HUnit--#if MIN_VERSION_base(4,8,0)-#else-import Data.Word-#endif--import Math.NumberTheory.Primes.Heap-import Math.NumberTheory.Primes.Testing-import Math.NumberTheory.TestUtils---- | Check that 'primes' over different integral types matches with 'isPrime'.-primesProperty1 :: Assertion-primesProperty1 = do-  assertEqual "ints  == integers" (trim ints)  (trim integers)-  assertEqual "words == integers" (trim words) (trim integers)-  assertEqual "naive == integers" (trim naive) (trim integers)-  where-    trim :: Integral a => [a] -> [Integer]-    trim = map toInteger . take 100000--    ints     = primes :: [Int]-    words    = primes :: [Word]-    integers = primes :: [Integer]-    naive    = filter isPrime [1..] :: [Integer]---- | Check that 'sieveFrom' over different integral types matches with 'isPrime'.-sieveFromProperty1 :: NonNegative Integer -> Bool-sieveFromProperty1 (NonNegative lowBound)-  =  trim ints  == trim integers-  && trim words == trim integers-  && trim naive == trim integers-  where-    trim :: Integral a => [a] -> [Integer]-    trim = map toInteger . take 1000--    ints     = sieveFrom (fromInteger lowBound) :: [Int]-    words    = sieveFrom (fromInteger lowBound) :: [Word]-    integers = sieveFrom lowBound               :: [Integer]-    naive    = filter isPrime [lowBound..]      :: [Integer]--testSuite :: TestTree-testSuite = testGroup "Heap"-  [ testCase          "primes"    primesProperty1-  , testSmallAndQuick "sieveFrom" sieveFromProperty1-  ]

@@ -11,6 +11,7 @@ {-# LANGUAGE CPP #-}  {-# OPTIONS_GHC -fno-warn-type-defaults #-}+{-# OPTIONS_GHC -fno-warn-deprecations  #-}  module Math.NumberTheory.Primes.SieveTests   ( testSuite

@@ -8,6 +8,7 @@ -- Tests for Math.NumberTheory.Primes -- +{-# OPTIONS_GHC -fno-warn-deprecations  #-} {-# OPTIONS_GHC -fno-warn-type-defaults #-}  module Math.NumberTheory.PrimesTests

@@ -47,7 +47,7 @@ import Test.Tasty.SmallCheck as SC import Test.Tasty.QuickCheck as QC hiding (Positive, NonNegative, generate, getNonNegative) -import Test.SmallCheck.Series (Positive(..), NonNegative(..), Serial(..), Series, generate)+import Test.SmallCheck.Series (Positive(..), NonNegative(..), Serial(..), Series, generate, (\/))  #if !(MIN_VERSION_base(4,8,0)) import Control.Applicative@@ -59,6 +59,8 @@ import Numeric.Natural  import Math.NumberTheory.GaussianIntegers (GaussianInteger(..))+import Math.NumberTheory.Moduli.PrimitiveRoot (CyclicGroup(..))+import Math.NumberTheory.UniqueFactorisation (UniqueFactorisation, Prime, unPrime)  import Math.NumberTheory.TestUtils.MyCompose import Math.NumberTheory.TestUtils.Wrappers@@ -88,6 +90,39 @@ instance Monad m => Serial m GaussianInteger where   series = cons2 (:+) +-------------------------------------------------------------------------------+-- Cyclic group++instance (Eq a, Num a, UniqueFactorisation a, Arbitrary a) => Arbitrary (CyclicGroup a) where+  arbitrary = frequency+    [ (1, pure CG2)+    , (1, pure CG4)+    , (9, CGOddPrimePower+      <$> (arbitrary :: Gen (PrimeWrapper a)) `suchThatMap` isOddPrime+      <*> (getPower <$> arbitrary))+    , (9, CGDoubleOddPrimePower+      <$> (arbitrary :: Gen (PrimeWrapper a)) `suchThatMap` isOddPrime+      <*> (getPower <$> arbitrary))+    ]++instance (Monad m, Eq a, Num a, UniqueFactorisation a, Serial m a) => Serial m (CyclicGroup a) where+  series = pure CG2+        \/ pure CG4+        \/ (CGOddPrimePower+           <$> (series :: Series m (PrimeWrapper a)) `suchThatMapSerial` isOddPrime+           <*> (getPower <$> series))+        \/ (CGDoubleOddPrimePower+           <$> (series :: Series m (PrimeWrapper a)) `suchThatMapSerial` isOddPrime+           <*> (getPower <$> series))++isOddPrime+  :: forall a. (Eq a, Num a, UniqueFactorisation a)+  => PrimeWrapper a+  -> Maybe (Prime a)+isOddPrime (PrimeWrapper p) = if (unPrime p :: a) == 2 then Nothing else Just p++-------------------------------------------------------------------------------+ -- https://www.cs.ox.ac.uk/projects/utgp/school/andres.pdf, p. 21 -- :k Compose = (k1 -> Constraint) -> (k2 -> k1) -> (k2 -> Constraint) class    (f (g x)) => (f `Compose` g) x@@ -114,7 +149,7 @@  testIntegralProperty   :: forall wrapper bool. (TestableIntegral wrapper, SC.Testable IO bool, QC.Testable bool)-  => String -> (forall a. (Integral a, Bits a) => wrapper a -> bool) -> TestTree+  => String -> (forall a. (Integral a, Bits a, UniqueFactorisation a, Show a) => wrapper a -> bool) -> TestTree testIntegralProperty name f = testGroup name   [ SC.testProperty "smallcheck Int"     (f :: wrapper Int     -> bool)   , SC.testProperty "smallcheck Word"    (f :: wrapper Word    -> bool)@@ -129,7 +164,7 @@  testSameIntegralProperty   :: forall wrapper1 wrapper2 bool. (TestableIntegral wrapper1, TestableIntegral wrapper2, SC.Testable IO bool, QC.Testable bool)-  => String -> (forall a. (Integral a, Bits a) => wrapper1 a -> wrapper2 a -> bool) -> TestTree+  => String -> (forall a. (Integral a, Bits a, UniqueFactorisation a, Show a) => wrapper1 a -> wrapper2 a -> bool) -> TestTree testSameIntegralProperty name f = testGroup name   [ SC.testProperty "smallcheck Int"     (f :: wrapper1 Int     -> wrapper2 Int     -> bool)   , SC.testProperty "smallcheck Word"    (f :: wrapper1 Word    -> wrapper2 Word    -> bool)@@ -144,7 +179,7 @@  testIntegral2Property   :: forall wrapper1 wrapper2 bool. (TestableIntegral wrapper1, TestableIntegral wrapper2, SC.Testable IO bool, QC.Testable bool)-  => String -> (forall a1 a2. (Integral a1, Integral a2, Bits a1, Bits a2) => wrapper1 a1 -> wrapper2 a2 -> bool) -> TestTree+  => String -> (forall a1 a2. (Integral a1, Integral a2, Bits a1, Bits a2, UniqueFactorisation a1, UniqueFactorisation a2, Show a1, Show a2) => wrapper1 a1 -> wrapper2 a2 -> bool) -> TestTree testIntegral2Property name f = testGroup name   [ SC.testProperty "smallcheck Int Int"         (f :: wrapper1 Int     -> wrapper2 Int     -> bool)   , SC.testProperty "smallcheck Int Word"        (f :: wrapper1 Int     -> wrapper2 Word    -> bool)

@@ -13,11 +13,14 @@ {-# LANGUAGE DeriveFoldable             #-} {-# LANGUAGE DeriveFunctor              #-} {-# LANGUAGE DeriveTraversable          #-}+{-# LANGUAGE FlexibleContexts           #-} {-# LANGUAGE FlexibleInstances          #-} {-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE MultiParamTypeClasses      #-} {-# LANGUAGE ScopedTypeVariables        #-} {-# LANGUAGE StandaloneDeriving         #-}+{-# LANGUAGE TypeFamilies               #-}+{-# LANGUAGE UndecidableInstances       #-}  {-# OPTIONS_GHC -fno-warn-orphans       #-} {-# OPTIONS_GHC -fno-warn-type-defaults #-}@@ -35,7 +38,7 @@ import Test.Tasty.QuickCheck as QC hiding (Positive, NonNegative, generate, getNonNegative, getPositive) import Test.SmallCheck.Series (Positive(..), NonNegative(..), Serial(..), Series) -import Math.NumberTheory.Primes (isPrime, nthPrime)+import Math.NumberTheory.UniqueFactorisation  ------------------------------------------------------------------------------- -- AnySign@@ -167,19 +170,40 @@ ------------------------------------------------------------------------------- -- Prime -newtype Prime = Prime { getPrime :: Integer }-  deriving (Eq, Ord, Show)+newtype PrimeWrapper a = PrimeWrapper { getPrime :: Prime a } -instance Arbitrary Prime where-  arbitrary = do-    n <- arbitrary-    return $ Prime $ head $ filter isPrime [abs n ..]+deriving instance Eq   (Prime a) => Eq   (PrimeWrapper a)+deriving instance Ord  (Prime a) => Ord  (PrimeWrapper a)+deriving instance Show (Prime a) => Show (PrimeWrapper a) -instance Monad m => Serial m Prime where-  series = Prime . nthPrime <$> series `suchThatSerial` (> 0)+instance (Arbitrary a, UniqueFactorisation a) => Arbitrary (PrimeWrapper a) where+  arbitrary = PrimeWrapper <$> (arbitrary :: Gen a) `suchThatMap` isPrime +instance (Monad m, Serial m a, UniqueFactorisation a) => Serial m (PrimeWrapper a) where+  series = PrimeWrapper <$> (series :: Series m a) `suchThatMapSerial` isPrime+ -------------------------------------------------------------------------------+-- UniqueFactorisation++type instance Prime (Large a) = Prime a++instance UniqueFactorisation a => UniqueFactorisation (Large a) where+  unPrime p = Large (unPrime p)+  factorise (Large x) = factorise x+  isPrime (Large x) = isPrime x++type instance Prime (Huge a) = Prime a++instance UniqueFactorisation a => UniqueFactorisation (Huge a) where+  unPrime p = Huge (unPrime p)+  factorise (Huge x) = factorise x+  isPrime (Huge x) = isPrime x++------------------------------------------------------------------------------- -- Utils  suchThatSerial :: Series m a -> (a -> Bool) -> Series m a suchThatSerial s p = s >>= \x -> if p x then pure x else empty++suchThatMapSerial :: Series m a -> (a -> Maybe b) -> Series m b+suchThatMapSerial s p = s >>= maybe empty pure . p

@@ -26,7 +26,7 @@  import Math.NumberTheory.GaussianIntegers hiding (factorise) import Math.NumberTheory.UniqueFactorisation-import Math.NumberTheory.TestUtils hiding (Prime)+import Math.NumberTheory.TestUtils  import Numeric.Natural 

@@ -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.JacobiTests as ModuliJacobi+import qualified Math.NumberTheory.Moduli.PrimitiveRootTests as ModuliPrimitiveRoot import qualified Math.NumberTheory.Moduli.SqrtTests as ModuliSqrt  import qualified Math.NumberTheory.MoebiusInversionTests as MoebiusInversion@@ -17,18 +18,22 @@ import qualified Math.NumberTheory.Powers.CubesTests as Cubes import qualified Math.NumberTheory.Powers.FourthTests as Fourth import qualified Math.NumberTheory.Powers.GeneralTests as General+import qualified Math.NumberTheory.Powers.ModularTests as Modular import qualified Math.NumberTheory.Powers.SquaresTests as Squares +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.HeapTests as Heap import qualified Math.NumberTheory.Primes.SieveTests as Sieve import qualified Math.NumberTheory.Primes.TestingTests as Testing  import qualified Math.NumberTheory.GaussianIntegersTests as Gaussian  import qualified Math.NumberTheory.ArithmeticFunctionsTests as ArithmeticFunctions+import qualified Math.NumberTheory.ArithmeticFunctions.MertensTests as Mertens+import qualified Math.NumberTheory.ArithmeticFunctions.SieveBlockTests as SieveBlock import qualified Math.NumberTheory.UniqueFactorisationTests as UniqueFactorisation import qualified Math.NumberTheory.ZetaTests as Zeta import qualified Math.NumberTheory.CurvesTests as Curves@@ -42,6 +47,7 @@     [ Cubes.testSuite     , Fourth.testSuite     , General.testSuite+    , Modular.testSuite     , Squares.testSuite     ]   , testGroup "GCD"@@ -56,17 +62,20 @@     [ ModuliChinese.testSuite     , ModuliClass.testSuite     , ModuliJacobi.testSuite+    , ModuliPrimitiveRoot.testSuite     , ModuliSqrt.testSuite     ]   , testGroup "MoebiusInversion"     [ MoebiusInversion.testSuite     , MoebiusInversionInt.testSuite     ]+  , testGroup "Prefactored"+    [ Prefactored.testSuite+    ]   , testGroup "Primes"     [ Primes.testSuite     , Counting.testSuite     , Factorisation.testSuite-    , Heap.testSuite     , Sieve.testSuite     , Testing.testSuite     ]@@ -75,6 +84,8 @@     ]   , testGroup "ArithmeticFunctions"     [ ArithmeticFunctions.testSuite+    , Mertens.testSuite+    , SieveBlock.testSuite     ]   , testGroup "UniqueFactorisation"     [ UniqueFactorisation.testSuite