diff --git a/Changes b/Changes --- a/Changes +++ b/Changes @@ -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: diff --git a/GHC/TypeNats/Compat.hs b/GHC/TypeNats/Compat.hs --- a/GHC/TypeNats/Compat.hs +++ b/GHC/TypeNats/Compat.hs @@ -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 diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -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, diff --git a/Math/NumberTheory/ArithmeticFunctions/Class.hs b/Math/NumberTheory/ArithmeticFunctions/Class.hs --- a/Math/NumberTheory/ArithmeticFunctions/Class.hs +++ b/Math/NumberTheory/ArithmeticFunctions/Class.hs @@ -21,7 +21,9 @@ ) where import Control.Applicative +#if __GLASGOW_HASKELL__ < 803 import Data.Semigroup +#endif #if MIN_VERSION_base(4,8,0) #else diff --git a/Math/NumberTheory/ArithmeticFunctions/Mertens.hs b/Math/NumberTheory/ArithmeticFunctions/Mertens.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/ArithmeticFunctions/Mertens.hs @@ -0,0 +1,78 @@ +-- | +-- Module: Math.NumberTheory.ArithmeticFunctions.Mertens +-- Copyright: (c) 2018 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- Stability: Provisional +-- Portability: Non-portable (GHC extensions) +-- +-- Values of . +-- + +{-# 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 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 diff --git a/Math/NumberTheory/ArithmeticFunctions/Moebius.hs b/Math/NumberTheory/ArithmeticFunctions/Moebius.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/ArithmeticFunctions/Moebius.hs @@ -0,0 +1,176 @@ +-- | +-- Module: Math.NumberTheory.ArithmeticFunctions.Moebius +-- Copyright: (c) 2018 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- Stability: Provisional +-- Portability: Non-portable (GHC extensions) +-- +-- Values of . +-- + +{-# 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 . +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 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 diff --git a/Math/NumberTheory/ArithmeticFunctions/SieveBlock.hs b/Math/NumberTheory/ArithmeticFunctions/SieveBlock.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/ArithmeticFunctions/SieveBlock.hs @@ -0,0 +1,132 @@ +-- | +-- Module: Math.NumberTheory.ArithmeticFunctions.SieveBlock +-- Copyright: (c) 2017 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- 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 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 diff --git a/Math/NumberTheory/ArithmeticFunctions/SieveBlock/Unboxed.hs b/Math/NumberTheory/ArithmeticFunctions/SieveBlock/Unboxed.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/ArithmeticFunctions/SieveBlock/Unboxed.hs @@ -0,0 +1,132 @@ +-- | +-- Module: Math.NumberTheory.ArithmeticFunctions.SieveBlock.Unboxed +-- Copyright: (c) 2017 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- 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 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 #-} diff --git a/Math/NumberTheory/ArithmeticFunctions/Standard.hs b/Math/NumberTheory/ArithmeticFunctions/Standard.hs --- a/Math/NumberTheory/ArithmeticFunctions/Standard.hs +++ b/Math/NumberTheory/ArithmeticFunctions/Standard.hs @@ -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 @ℤ/(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 diff --git a/Math/NumberTheory/GCD.hs b/Math/NumberTheory/GCD.hs --- a/Math/NumberTheory/GCD.hs +++ b/Math/NumberTheory/GCD.hs @@ -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) diff --git a/Math/NumberTheory/GaussianIntegers.hs b/Math/NumberTheory/GaussianIntegers.hs --- a/Math/NumberTheory/GaussianIntegers.hs +++ b/Math/NumberTheory/GaussianIntegers.hs @@ -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 -- 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 diff --git a/Math/NumberTheory/Moduli/Class.hs b/Math/NumberTheory/Moduli/Class.hs --- a/Math/NumberTheory/Moduli/Class.hs +++ b/Math/NumberTheory/Moduli/Class.hs @@ -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 :: diff --git a/Math/NumberTheory/Moduli/Jacobi.hs b/Math/NumberTheory/Moduli/Jacobi.hs --- a/Math/NumberTheory/Moduli/Jacobi.hs +++ b/Math/NumberTheory/Moduli/Jacobi.hs @@ -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 diff --git a/Math/NumberTheory/Moduli/PrimitiveRoot.hs b/Math/NumberTheory/Moduli/PrimitiveRoot.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Moduli/PrimitiveRoot.hs @@ -0,0 +1,148 @@ +-- | +-- Module: Math.NumberTheory.Moduli.PrimitiveRoot +-- Copyright: (c) 2017 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- 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 +-- 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 . +-- +-- > > 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) diff --git a/Math/NumberTheory/MoebiusInversion.hs b/Math/NumberTheory/MoebiusInversion.hs --- a/Math/NumberTheory/MoebiusInversion.hs +++ b/Math/NumberTheory/MoebiusInversion.hs @@ -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 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 diff --git a/Math/NumberTheory/MoebiusInversion/Int.hs b/Math/NumberTheory/MoebiusInversion/Int.hs --- a/Math/NumberTheory/MoebiusInversion/Int.hs +++ b/Math/NumberTheory/MoebiusInversion/Int.hs @@ -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 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 diff --git a/Math/NumberTheory/Powers.hs b/Math/NumberTheory/Powers.hs --- a/Math/NumberTheory/Powers.hs +++ b/Math/NumberTheory/Powers.hs @@ -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 diff --git a/Math/NumberTheory/Powers/Integer.hs b/Math/NumberTheory/Powers/Integer.hs deleted file mode 100644 --- a/Math/NumberTheory/Powers/Integer.hs +++ /dev/null @@ -1,44 +0,0 @@ --- | --- Module: Math.NumberTheory.Powers.Integer --- Copyright: (c) 2011-2014 Daniel Fischer --- Licence: MIT --- Maintainer: Daniel Fischer --- 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" #-} diff --git a/Math/NumberTheory/Powers/Modular.hs b/Math/NumberTheory/Powers/Modular.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Powers/Modular.hs @@ -0,0 +1,97 @@ +-- | +-- Module: Math.NumberTheory.Powers.Modular +-- Copyright: (c) 2017 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- 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 diff --git a/Math/NumberTheory/Prefactored.hs b/Math/NumberTheory/Prefactored.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Prefactored.hs @@ -0,0 +1,130 @@ +-- | +-- Module: Math.NumberTheory.Prefactored +-- Copyright: (c) 2017 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- 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 diff --git a/Math/NumberTheory/Primes/Factorisation/Montgomery.hs b/Math/NumberTheory/Primes/Factorisation/Montgomery.hs --- a/Math/NumberTheory/Primes/Factorisation/Montgomery.hs +++ b/Math/NumberTheory/Primes/Factorisation/Montgomery.hs @@ -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 diff --git a/Math/NumberTheory/Primes/Heap.hs b/Math/NumberTheory/Primes/Heap.hs --- a/Math/NumberTheory/Primes/Heap.hs +++ b/Math/NumberTheory/Primes/Heap.hs @@ -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 diff --git a/Math/NumberTheory/Primes/Sieve/Misc.hs b/Math/NumberTheory/Primes/Sieve/Misc.hs --- a/Math/NumberTheory/Primes/Sieve/Misc.hs +++ b/Math/NumberTheory/Primes/Sieve/Misc.hs @@ -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! diff --git a/Math/NumberTheory/Primes/Testing.hs b/Math/NumberTheory/Primes/Testing.hs --- a/Math/NumberTheory/Primes/Testing.hs +++ b/Math/NumberTheory/Primes/Testing.hs @@ -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'. diff --git a/Math/NumberTheory/Primes/Types.hs b/Math/NumberTheory/Primes/Types.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Primes/Types.hs @@ -0,0 +1,53 @@ +-- | +-- Module: Math.NumberTheory.Primes.Types +-- Copyright: (c) 2017 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- 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 diff --git a/Math/NumberTheory/SmoothNumbers.hs b/Math/NumberTheory/SmoothNumbers.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/SmoothNumbers.hs @@ -0,0 +1,130 @@ +-- | +-- Module: Math.NumberTheory.SmoothNumbers +-- Copyright: (c) 2018 Frederick Schneider +-- Licence: MIT +-- Maintainer: Frederick Schneider +-- Stability: Provisional +-- Portability: Non-portable (GHC extensions) +-- +-- A +-- 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 +-- +-- 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 +-- +-- 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 +-- +-- 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 + diff --git a/Math/NumberTheory/UniqueFactorisation.hs b/Math/NumberTheory/UniqueFactorisation.hs --- a/Math/NumberTheory/UniqueFactorisation.hs +++ b/Math/NumberTheory/UniqueFactorisation.hs @@ -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 diff --git a/Math/NumberTheory/Utils.hs b/Math/NumberTheory/Utils.hs --- a/Math/NumberTheory/Utils.hs +++ b/Math/NumberTheory/Utils.hs @@ -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# #-} diff --git a/Math/NumberTheory/Utils/FromIntegral.hs b/Math/NumberTheory/Utils/FromIntegral.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Utils/FromIntegral.hs @@ -0,0 +1,55 @@ +-- | +-- Module: Math.NumberTheory.Utils.FromIntegral +-- Copyright: (c) 2017 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- 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 #-} diff --git a/Math/NumberTheory/Utils/Hyperbola.hs b/Math/NumberTheory/Utils/Hyperbola.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Utils/Hyperbola.hs @@ -0,0 +1,92 @@ +-- | +-- Module: Math.NumberTheory.Utils.Hyperbola +-- Copyright: (c) 2018 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- 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 +-- 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) diff --git a/arithmoi.cabal b/arithmoi.cabal --- a/arithmoi.cabal +++ b/arithmoi.cabal @@ -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 diff --git a/benchmark/Bench.hs b/benchmark/Bench.hs --- a/benchmark/Bench.hs +++ b/benchmark/Bench.hs @@ -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 ] diff --git a/benchmark/Math/NumberTheory/ArithmeticFunctionsBench.hs b/benchmark/Math/NumberTheory/ArithmeticFunctionsBench.hs --- a/benchmark/Math/NumberTheory/ArithmeticFunctionsBench.hs +++ b/benchmark/Math/NumberTheory/ArithmeticFunctionsBench.hs @@ -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) diff --git a/benchmark/Math/NumberTheory/GCDBench.hs b/benchmark/Math/NumberTheory/GCDBench.hs new file mode 100644 --- /dev/null +++ b/benchmark/Math/NumberTheory/GCDBench.hs @@ -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 + ] diff --git a/benchmark/Math/NumberTheory/MertensBench.hs b/benchmark/Math/NumberTheory/MertensBench.hs new file mode 100644 --- /dev/null +++ b/benchmark/Math/NumberTheory/MertensBench.hs @@ -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 diff --git a/benchmark/Math/NumberTheory/PowersBench.hs b/benchmark/Math/NumberTheory/PowersBench.hs --- a/benchmark/Math/NumberTheory/PowersBench.hs +++ b/benchmark/Math/NumberTheory/PowersBench.hs @@ -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] diff --git a/benchmark/Math/NumberTheory/PrimesBench.hs b/benchmark/Math/NumberTheory/PrimesBench.hs --- a/benchmark/Math/NumberTheory/PrimesBench.hs +++ b/benchmark/Math/NumberTheory/PrimesBench.hs @@ -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) diff --git a/benchmark/Math/NumberTheory/RecurrenciesBench.hs b/benchmark/Math/NumberTheory/RecurrenciesBench.hs --- a/benchmark/Math/NumberTheory/RecurrenciesBench.hs +++ b/benchmark/Math/NumberTheory/RecurrenciesBench.hs @@ -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 diff --git a/benchmark/Math/NumberTheory/SieveBlockBench.hs b/benchmark/Math/NumberTheory/SieveBlockBench.hs new file mode 100644 --- /dev/null +++ b/benchmark/Math/NumberTheory/SieveBlockBench.hs @@ -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 + ] + ] diff --git a/test-suite/Math/NumberTheory/ArithmeticFunctions/MertensTests.hs b/test-suite/Math/NumberTheory/ArithmeticFunctions/MertensTests.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/NumberTheory/ArithmeticFunctions/MertensTests.hs @@ -0,0 +1,76 @@ +-- | +-- Module: Math.NumberTheory.ArithmeticFunctions.MertensTests +-- Copyright: (c) 2018 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- 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 + ] diff --git a/test-suite/Math/NumberTheory/ArithmeticFunctions/SieveBlockTests.hs b/test-suite/Math/NumberTheory/ArithmeticFunctions/SieveBlockTests.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/NumberTheory/ArithmeticFunctions/SieveBlockTests.hs @@ -0,0 +1,109 @@ +-- | +-- Module: Math.NumberTheory.ArithmeticFunctions.SieveBlockTests +-- Copyright: (c) 2016 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- 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 + ] diff --git a/test-suite/Math/NumberTheory/ArithmeticFunctionsTests.hs b/test-suite/Math/NumberTheory/ArithmeticFunctionsTests.hs --- a/test-suite/Math/NumberTheory/ArithmeticFunctionsTests.hs +++ b/test-suite/Math/NumberTheory/ArithmeticFunctionsTests.hs @@ -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" diff --git a/test-suite/Math/NumberTheory/GCDTests.hs b/test-suite/Math/NumberTheory/GCDTests.hs --- a/test-suite/Math/NumberTheory/GCDTests.hs +++ b/test-suite/Math/NumberTheory/GCDTests.hs @@ -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 + ] ] diff --git a/test-suite/Math/NumberTheory/GaussianIntegersTests.hs b/test-suite/Math/NumberTheory/GaussianIntegersTests.hs --- a/test-suite/Math/NumberTheory/GaussianIntegersTests.hs +++ b/test-suite/Math/NumberTheory/GaussianIntegersTests.hs @@ -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 -- 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 diff --git a/test-suite/Math/NumberTheory/Moduli/JacobiTests.hs b/test-suite/Math/NumberTheory/Moduli/JacobiTests.hs --- a/test-suite/Math/NumberTheory/Moduli/JacobiTests.hs +++ b/test-suite/Math/NumberTheory/Moduli/JacobiTests.hs @@ -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 diff --git a/test-suite/Math/NumberTheory/Moduli/PrimitiveRootTests.hs b/test-suite/Math/NumberTheory/Moduli/PrimitiveRootTests.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/NumberTheory/Moduli/PrimitiveRootTests.hs @@ -0,0 +1,118 @@ +-- | +-- Module: Math.NumberTheory.Moduli.PrimitiveRootTests +-- Copyright: (c) 2017 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- 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 + ] + ] diff --git a/test-suite/Math/NumberTheory/Moduli/SqrtTests.hs b/test-suite/Math/NumberTheory/Moduli/SqrtTests.hs --- a/test-suite/Math/NumberTheory/Moduli/SqrtTests.hs +++ b/test-suite/Math/NumberTheory/Moduli/SqrtTests.hs @@ -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) diff --git a/test-suite/Math/NumberTheory/MoebiusInversion/IntTests.hs b/test-suite/Math/NumberTheory/MoebiusInversion/IntTests.hs --- a/test-suite/Math/NumberTheory/MoebiusInversion/IntTests.hs +++ b/test-suite/Math/NumberTheory/MoebiusInversion/IntTests.hs @@ -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 diff --git a/test-suite/Math/NumberTheory/MoebiusInversionTests.hs b/test-suite/Math/NumberTheory/MoebiusInversionTests.hs --- a/test-suite/Math/NumberTheory/MoebiusInversionTests.hs +++ b/test-suite/Math/NumberTheory/MoebiusInversionTests.hs @@ -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 diff --git a/test-suite/Math/NumberTheory/Powers/ModularTests.hs b/test-suite/Math/NumberTheory/Powers/ModularTests.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/NumberTheory/Powers/ModularTests.hs @@ -0,0 +1,98 @@ +-- | +-- Module: Math.NumberTheory.Powers.ModularTests +-- Copyright: (c) 2017 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- 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 + ] + ] diff --git a/test-suite/Math/NumberTheory/PrefactoredTests.hs b/test-suite/Math/NumberTheory/PrefactoredTests.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/NumberTheory/PrefactoredTests.hs @@ -0,0 +1,101 @@ +-- | +-- Module: Math.NumberTheory.PrefactoredTests +-- Copyright: (c) 2017 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- 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 + ] + ] diff --git a/test-suite/Math/NumberTheory/Primes/FactorisationTests.hs b/test-suite/Math/NumberTheory/Primes/FactorisationTests.hs --- a/test-suite/Math/NumberTheory/Primes/FactorisationTests.hs +++ b/test-suite/Math/NumberTheory/Primes/FactorisationTests.hs @@ -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 ] diff --git a/test-suite/Math/NumberTheory/Primes/HeapTests.hs b/test-suite/Math/NumberTheory/Primes/HeapTests.hs deleted file mode 100644 --- a/test-suite/Math/NumberTheory/Primes/HeapTests.hs +++ /dev/null @@ -1,67 +0,0 @@ --- | --- Module: Math.NumberTheory.Primes.HeapTests --- Copyright: (c) 2016 Andrew Lelechenko --- Licence: MIT --- Maintainer: Andrew Lelechenko --- 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 - ] diff --git a/test-suite/Math/NumberTheory/Primes/SieveTests.hs b/test-suite/Math/NumberTheory/Primes/SieveTests.hs --- a/test-suite/Math/NumberTheory/Primes/SieveTests.hs +++ b/test-suite/Math/NumberTheory/Primes/SieveTests.hs @@ -11,6 +11,7 @@ {-# LANGUAGE CPP #-} {-# OPTIONS_GHC -fno-warn-type-defaults #-} +{-# OPTIONS_GHC -fno-warn-deprecations #-} module Math.NumberTheory.Primes.SieveTests ( testSuite diff --git a/test-suite/Math/NumberTheory/PrimesTests.hs b/test-suite/Math/NumberTheory/PrimesTests.hs --- a/test-suite/Math/NumberTheory/PrimesTests.hs +++ b/test-suite/Math/NumberTheory/PrimesTests.hs @@ -8,6 +8,7 @@ -- Tests for Math.NumberTheory.Primes -- +{-# OPTIONS_GHC -fno-warn-deprecations #-} {-# OPTIONS_GHC -fno-warn-type-defaults #-} module Math.NumberTheory.PrimesTests diff --git a/test-suite/Math/NumberTheory/TestUtils.hs b/test-suite/Math/NumberTheory/TestUtils.hs --- a/test-suite/Math/NumberTheory/TestUtils.hs +++ b/test-suite/Math/NumberTheory/TestUtils.hs @@ -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) diff --git a/test-suite/Math/NumberTheory/TestUtils/Wrappers.hs b/test-suite/Math/NumberTheory/TestUtils/Wrappers.hs --- a/test-suite/Math/NumberTheory/TestUtils/Wrappers.hs +++ b/test-suite/Math/NumberTheory/TestUtils/Wrappers.hs @@ -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 diff --git a/test-suite/Math/NumberTheory/UniqueFactorisationTests.hs b/test-suite/Math/NumberTheory/UniqueFactorisationTests.hs --- a/test-suite/Math/NumberTheory/UniqueFactorisationTests.hs +++ b/test-suite/Math/NumberTheory/UniqueFactorisationTests.hs @@ -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 diff --git a/test-suite/Test.hs b/test-suite/Test.hs --- a/test-suite/Test.hs +++ b/test-suite/Test.hs @@ -9,6 +9,7 @@ import qualified Math.NumberTheory.Moduli.ChineseTests as ModuliChinese import qualified Math.NumberTheory.Moduli.ClassTests as ModuliClass import qualified Math.NumberTheory.Moduli.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