diff --git a/Changes b/Changes --- a/Changes +++ b/Changes @@ -1,5 +1,19 @@ +0.4.3.0: + This release supports GHC 7.6, 7.8, 7.10 and 8.0. + + Add Math.NumberTheory.ArithmeticFunctions with brand-new machinery + for arithmetic functions: divisors, tau, sigma, totient, jordan, + moebius, liouville, smallOmega, bigOmega, carmichael, expMangoldt (#30). + Old implementations (exposed via Math.NumberTheory.Primes.Factorisation + and Math.NumberTheory.Powers.Integer) are deprecated and will be removed + in the next major release. + + Add Karatsuba sqrt algorithm, improving performance on large integers (#6). + + Fix incorrect indexing of FactorSieve (#35). + 0.4.2.0: - This release supports GHC 7.6, 7.8 and 8.0. + This release supports GHC 7.6, 7.8, 7.10 and 8.0. Add new cabal flag check-bounds, which replaces all unsafe array functions with safe ones. diff --git a/Math/NumberTheory/ArithmeticFunctions.hs b/Math/NumberTheory/ArithmeticFunctions.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/ArithmeticFunctions.hs @@ -0,0 +1,20 @@ +-- | +-- Module: Math.NumberTheory.ArithmeticFunctions +-- Copyright: (c) 2016 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- Stability: Provisional +-- Portability: Non-portable (GHC extensions) +-- +-- This module provides an interface for defining and manipulating +-- arithmetic functions. It also defines several most widespreaded +-- arithmetic functions. +-- + +module Math.NumberTheory.ArithmeticFunctions + ( module Math.NumberTheory.ArithmeticFunctions.Class + , module Math.NumberTheory.ArithmeticFunctions.Standard + ) where + +import Math.NumberTheory.ArithmeticFunctions.Class +import Math.NumberTheory.ArithmeticFunctions.Standard diff --git a/Math/NumberTheory/ArithmeticFunctions/Class.hs b/Math/NumberTheory/ArithmeticFunctions/Class.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/ArithmeticFunctions/Class.hs @@ -0,0 +1,102 @@ +-- | +-- Module: Math.NumberTheory.ArithmeticFunctions.Class +-- Copyright: (c) 2016 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- Stability: Provisional +-- Portability: Non-portable (GHC extensions) +-- +-- Generic type for arithmetic functions over arbitrary unique +-- factorisation domains. +-- + +{-# LANGUAGE CPP #-} +{-# LANGUAGE GADTs #-} + +{-# OPTIONS_HADDOCK hide #-} + +module Math.NumberTheory.ArithmeticFunctions.Class + ( ArithmeticFunction(..) + , runFunction + ) where + +import Control.Applicative +import Data.Semigroup + +#if MIN_VERSION_base(4,8,0) +#else +import Data.Word +#endif + +import Math.NumberTheory.UniqueFactorisation + +-- | A typical arithmetic function operates on the canonical factorisation of +-- a number into prime's powers and consists of two rules. The first one +-- determines the values of the function on the powers of primes. The second +-- one determines how to combine these values into final result. +-- +-- In the following definition the first argument is the function on prime's +-- powers, the monoid instance determines a rule of combination (typically +-- 'Product' or 'Sum'), and the second argument is convenient for unwrapping +-- (typically, 'getProduct' or 'getSum'). +data ArithmeticFunction n a where + ArithmeticFunction + :: Monoid m + => (Prime n -> Word -> m) + -> (m -> a) + -> ArithmeticFunction n a + +-- | Convert to function. The value on 0 is undefined. +runFunction :: UniqueFactorisation n => ArithmeticFunction n a -> n -> a +runFunction (ArithmeticFunction f g) + = g + . mconcat + . map (uncurry f) + . factorise + +instance Functor (ArithmeticFunction n) where + fmap f (ArithmeticFunction g h) = ArithmeticFunction g (f . h) + +instance Applicative (ArithmeticFunction n) where + pure x + = ArithmeticFunction (\_ _ -> ()) (const x) + (ArithmeticFunction f1 g1) <*> (ArithmeticFunction f2 g2) + = ArithmeticFunction (\p k -> (f1 p k, f2 p k)) (\(a1, a2) -> g1 a1 (g2 a2)) + +instance Semigroup a => Semigroup (ArithmeticFunction n a) where + (<>) = liftA2 (<>) + +instance Monoid a => Monoid (ArithmeticFunction n a) where + mempty = pure mempty + mappend = liftA2 mappend + +-- | Factorisation is expensive, so it is better to avoid doing it twice. +-- Write 'runFunction (f + g) n' instead of 'runFunction f n + runFunction g n'. +instance Num a => Num (ArithmeticFunction n a) where + fromInteger = pure . fromInteger + negate = fmap negate + signum = fmap signum + abs = fmap abs + (+) = liftA2 (+) + (-) = liftA2 (-) + (*) = liftA2 (*) + +instance Fractional a => Fractional (ArithmeticFunction n a) where + fromRational = pure . fromRational + recip = fmap recip + (/) = liftA2 (/) + +instance Floating a => Floating (ArithmeticFunction n a) where + pi = pure pi + exp = fmap exp + log = fmap log + sin = fmap sin + cos = fmap cos + asin = fmap asin + acos = fmap acos + atan = fmap atan + sinh = fmap sinh + cosh = fmap cosh + asinh = fmap asinh + acosh = fmap acosh + atanh = fmap atanh diff --git a/Math/NumberTheory/ArithmeticFunctions/Standard.hs b/Math/NumberTheory/ArithmeticFunctions/Standard.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/ArithmeticFunctions/Standard.hs @@ -0,0 +1,299 @@ +-- | +-- Module: Math.NumberTheory.ArithmeticFunctions.Standard +-- Copyright: (c) 2016 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- Stability: Provisional +-- Portability: Non-portable (GHC extensions) +-- +-- Textbook arithmetic functions. +-- + +{-# LANGUAGE CPP #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE TypeFamilies #-} +{-# LANGUAGE ViewPatterns #-} + +{-# OPTIONS_HADDOCK hide #-} + +module Math.NumberTheory.ArithmeticFunctions.Standard + ( -- * Multiplicative functions + multiplicative + , divisors, divisorsA + , divisorsSmall, divisorsSmallA + , tau, tauA + , sigma, sigmaA + , totient, totientA + , jordan, jordanA + , moebius, moebiusA + , liouville, liouvilleA + -- * Additive functions + , additive + , smallOmega, smallOmegaA + , bigOmega, bigOmegaA + -- * Misc + , carmichael, carmichaelA + , expMangoldt, expMangoldtA + ) where + +import Data.IntSet (IntSet) +import qualified Data.IntSet as IS +import Data.Set (Set) +import qualified Data.Set as S +import Data.Semigroup + +import Math.NumberTheory.ArithmeticFunctions.Class +import Math.NumberTheory.UniqueFactorisation + +import Numeric.Natural + +#if MIN_VERSION_base(4,8,0) +#else +import Data.Foldable +import Data.Word +#endif + +#if MIN_VERSION_base(4,7,0) +import Data.Coerce +#else +import Unsafe.Coerce + +coerce :: a -> b +coerce = unsafeCoerce +#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 + +divisors :: (UniqueFactorisation n, Num n, Ord n) => n -> Set n +divisors = runFunction divisorsA +{-# SPECIALIZE divisors :: Natural -> Set Natural #-} +{-# SPECIALIZE divisors :: Integer -> Set Integer #-} + +-- | The set of all (positive) divisors of an argument. +divisorsA :: forall n. (UniqueFactorisation n, Num n, Ord n) => ArithmeticFunction n (Set n) +divisorsA = ArithmeticFunction (\((unPrime :: Prime n -> n) -> p) k -> SetProduct $ divisorsHelper p k) (S.insert 1 . getSetProduct) + +divisorsHelper :: Num n => n -> Word -> Set n +divisorsHelper _ 0 = S.empty +divisorsHelper p 1 = S.singleton p +divisorsHelper p a = S.fromDistinctAscList $ p : p * p : map (p ^) [3 .. wordToInt a] +{-# INLINE divisorsHelper #-} + +divisorsSmall :: (UniqueFactorisation n, Prime n ~ Prime Int) => n -> IntSet +divisorsSmall = runFunction divisorsSmallA + +-- | Same as 'divisors', but with better performance on cost of type restriction. +divisorsSmallA :: forall n. (Prime n ~ Prime Int) => ArithmeticFunction n IntSet +divisorsSmallA = ArithmeticFunction (\p k -> IntSetProduct $ divisorsHelperSmall (unPrime p) k) (IS.insert 1 . getIntSetProduct) + +divisorsHelperSmall :: Int -> Word -> IntSet +divisorsHelperSmall _ 0 = IS.empty +divisorsHelperSmall p 1 = IS.singleton p +divisorsHelperSmall p a = IS.fromDistinctAscList $ p : p * p : map (p ^) [3 .. wordToInt a] +{-# INLINE divisorsHelperSmall #-} + +tau :: (UniqueFactorisation n, Num a) => n -> a +tau = runFunction tauA + +-- | The number of (positive) divisors of an argument. +-- +-- > tauA = multiplicative (\_ k -> k + 1) +tauA :: Num a => ArithmeticFunction n a +tauA = multiplicative $ const (fromIntegral . succ) + +sigma :: (UniqueFactorisation n, Integral n) => Word -> n -> n +sigma = runFunction . sigmaA + +-- | The sum of the @k@-th powers of (positive) divisors of an argument. +-- +-- > sigmaA = multiplicative (\p k -> sum $ map (p ^) [0..k]) +-- > sigmaA 0 = tauA +sigmaA :: forall n. (UniqueFactorisation n, Integral n) => Word -> ArithmeticFunction n n +sigmaA 0 = tauA +sigmaA 1 = multiplicative $ \((unPrime :: Prime n -> n) -> p) -> sigmaHelper p +sigmaA a = multiplicative $ \((unPrime :: Prime n -> n) -> p) -> sigmaHelper (p ^ wordToInt a) + +sigmaHelper :: Integral n => n -> Word -> n +sigmaHelper pa 1 = pa + 1 +sigmaHelper pa 2 = pa * pa + pa + 1 +sigmaHelper pa k = (pa ^ wordToInt (k + 1) - 1) `quot` (pa - 1) +{-# INLINE sigmaHelper #-} + +totient :: (UniqueFactorisation n, Integral 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 = multiplicative $ \((unPrime :: Prime n -> n) -> p) -> jordanHelper p + +jordan :: (UniqueFactorisation n, Integral 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 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 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 = runFunction moebiusA + +-- | Calculates the Moebius function of an argument. +moebiusA :: Num a => ArithmeticFunction n a +moebiusA = ArithmeticFunction (const f) runMoebius + where + f 1 = MoebiusN + f 0 = MoebiusP + f _ = MoebiusZ + +liouville :: (UniqueFactorisation n, Num a) => n -> a +liouville = runFunction liouvilleA + +-- | Calculates the Liouville function of an argument. +liouvilleA :: Num a => ArithmeticFunction n a +liouvilleA = ArithmeticFunction (const $ Xor . odd) runXor + +carmichael :: (UniqueFactorisation n, Integral n) => n -> n +carmichael = runFunction carmichaelA +{- The specializations reflects available specializations of lcm. -} +{-# SPECIALIZE carmichael :: Int -> Int #-} +{-# SPECIALIZE carmichael :: Integer -> Integer #-} + +-- | Calculates the Carmichael function for a positive integer, that is, +-- the (smallest) exponent of the group of units in @ℤ/(n)@. +carmichaelA :: forall n. (UniqueFactorisation n, Integral n) => ArithmeticFunction n n +carmichaelA = ArithmeticFunction (\((unPrime :: Prime n -> n) -> p) k -> LCM $ f p k) getLCM + where + f 2 1 = 1 + f 2 2 = 2 + f 2 k = 2 ^ wordToInt (k - 2) + f p 1 = p - 1 + f p 2 = (p - 1) * p + f p k = (p - 1) * p ^ wordToInt (k - 1) + +-- | Create an additive function from the function on prime's powers. See examples below. +additive :: Num a => (Prime n -> Word -> a) -> ArithmeticFunction n a +additive f = ArithmeticFunction ((Sum .) . f) getSum + +smallOmega :: (UniqueFactorisation n, Num a) => n -> a +smallOmega = runFunction smallOmegaA + +-- | Number of distinct prime factors. +-- +-- > smallOmegaA = additive (\_ _ -> 1) +smallOmegaA :: Num a => ArithmeticFunction n a +smallOmegaA = additive (\_ _ -> 1) + +bigOmega :: UniqueFactorisation n => n -> Word +bigOmega = runFunction bigOmegaA + +-- | Number of prime factors, counted with multiplicity. +-- +-- > bigOmegaA = additive (\_ k -> k) +bigOmegaA :: ArithmeticFunction n Word +bigOmegaA = additive $ const id + +expMangoldt :: (UniqueFactorisation n, Num n) => n -> n +expMangoldt = runFunction expMangoldtA + +-- | 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 + | MangoldtOne a + | MangoldtMany + +runMangoldt :: Num a => Mangoldt a -> a +runMangoldt m = case m of + MangoldtZero -> 1 + MangoldtOne a -> a + MangoldtMany -> 1 + +instance Semigroup (Mangoldt a) where + MangoldtZero <> a = a + a <> MangoldtZero = a + _ <> _ = MangoldtMany + +instance Monoid (Mangoldt a) where + mempty = MangoldtZero + mappend = (<>) + +newtype LCM a = LCM { getLCM :: a } + +instance Integral a => Semigroup (LCM a) where + (<>) = coerce (lcm :: a -> a -> a) + +instance Integral a => Monoid (LCM a) where + mempty = LCM 1 + mappend = (<>) + +newtype Xor = Xor { _getXor :: Bool } + +runXor :: Num a => Xor -> a +runXor m = case m of + Xor False -> 1 + Xor True -> -1 + +instance Semigroup Xor where + (<>) = coerce ((/=) :: Bool -> Bool -> Bool) + +instance Monoid Xor where + mempty = Xor False + mappend = (<>) + +newtype SetProduct a = SetProduct { getSetProduct :: Set a } + +instance (Num a, Ord a) => Semigroup (SetProduct a) where + SetProduct s1 <> SetProduct s2 = SetProduct $ s1 <> s2 <> foldMap (\n -> S.mapMonotonic (* n) s2) s1 + +instance (Num a, Ord a) => Monoid (SetProduct a) where + mempty = SetProduct mempty + mappend = (<>) + +newtype IntSetProduct = IntSetProduct { getIntSetProduct :: IntSet } + +instance Semigroup IntSetProduct where + IntSetProduct s1 <> IntSetProduct s2 = IntSetProduct $ IS.unions $ s1 : s2 : map (\n -> IS.map (* n) s2) (IS.toAscList s1) + +instance Monoid IntSetProduct where + mempty = IntSetProduct mempty + mappend = (<>) diff --git a/Math/NumberTheory/GCD/LowLevel.hs b/Math/NumberTheory/GCD/LowLevel.hs --- a/Math/NumberTheory/GCD/LowLevel.hs +++ b/Math/NumberTheory/GCD/LowLevel.hs @@ -10,7 +10,8 @@ -- Normally, accessing these via the higher level interface of "Math.NumberTheory.GCD" -- should be sufficient. -- -{-# LANGUAGE CPP, MagicHash, UnboxedTuples #-} +{-# LANGUAGE MagicHash #-} +{-# LANGUAGE UnboxedTuples #-} module Math.NumberTheory.GCD.LowLevel ( -- * Specialised GCDs gcdInt @@ -27,10 +28,6 @@ ) where import GHC.Base - -#if __GLASGOW_HASKELL__ < 705 -import GHC.Word (Word(..)) -- Moved to GHC.Types -#endif import Math.NumberTheory.Utils diff --git a/Math/NumberTheory/GaussianIntegers.hs b/Math/NumberTheory/GaussianIntegers.hs --- a/Math/NumberTheory/GaussianIntegers.hs +++ b/Math/NumberTheory/GaussianIntegers.hs @@ -11,6 +11,7 @@ -- {-# LANGUAGE BangPatterns #-} + module Math.NumberTheory.GaussianIntegers ( GaussianInteger((:+)), ι, @@ -187,7 +188,7 @@ where s = a .^ div e 2 --- |Compute the prime factorization of a Gaussian integer. This is unique up to units (+/- 1, +/- i). +-- |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" diff --git a/Math/NumberTheory/Logarithms.hs b/Math/NumberTheory/Logarithms.hs --- a/Math/NumberTheory/Logarithms.hs +++ b/Math/NumberTheory/Logarithms.hs @@ -30,14 +30,11 @@ import GHC.Base -#if __GLASGOW_HASKELL__ < 705 -import GHC.Word (Word(..)) -- Moved to GHC.Types -#endif - import Data.Bits import Data.Array.Unboxed -import Math.NumberTheory.Logarithms.Internal +import GHC.Integer.Logarithms + import Math.NumberTheory.Powers.Integer import Math.NumberTheory.Unsafe #if __GLASGOW_HASKELL__ < 707 diff --git a/Math/NumberTheory/Logarithms/Internal.hs b/Math/NumberTheory/Logarithms/Internal.hs deleted file mode 100644 --- a/Math/NumberTheory/Logarithms/Internal.hs +++ /dev/null @@ -1,155 +0,0 @@ --- | --- Module: Math.NumberTheory.Logarithms.Internal --- Copyright: (c) 2011 Daniel Fischer --- Licence: MIT --- Maintainer: Daniel Fischer --- Stability: Provisional --- Portability: Non-portable (GHC extensions) --- --- Low level stuff for integer logarithms. -{-# LANGUAGE CPP, MagicHash, UnboxedTuples #-} -{-# OPTIONS_HADDOCK hide #-} -module Math.NumberTheory.Logarithms.Internal - ( -- * Functions - integerLog2# - , wordLog2# - ) where - -#if __GLASGOW_HASKELL__ >= 702 - --- Stuff is already there -import GHC.Integer.Logarithms - -#else - --- We have to define it here -#include "MachDeps.h" - -import GHC.Base -import GHC.Integer.GMP.Internals - -#if (WORD_SIZE_IN_BITS != 32) && (WORD_SIZE_IN_BITS != 64) -#error Only word sizes 32 and 64 are supported. -#endif - - -#if WORD_SIZE_IN_BITS == 32 - -#define WSHIFT 5 -#define MMASK 31 - -#else - -#define WSHIFT 6 -#define MMASK 63 - -#endif - -{- - -Reference implementation only, the algorithm in M.NT.Logarithms is better. - --- | Calculate the integer logarithm for an arbitrary base. --- The base must be greater than 1, the second argument, the number --- whose logarithm is sought; should be positive, otherwise the --- result is meaningless. --- --- > base ^ integerLogBase# base m <= m < base ^ (integerLogBase# base m + 1) --- --- for @base > 1@ and @m > 0@. -integerLogBase# :: Integer -> Integer -> Int# -integerLogBase# b m = case step b of - (# _, e #) -> e - where - step pw = - if m `ltInteger` pw - then (# m, 0# #) - else case step (pw `timesInteger` pw) of - (# q, e #) -> - if q `ltInteger` pw - then (# q, 2# *# e #) - else (# q `quotInteger` pw, 2# *# e +# 1# #) --} - --- | Calculate the integer base 2 logarithm of an 'Integer'. --- The calculation is much more efficient than for the general case. --- --- The argument must be strictly positive, that condition is /not/ checked. -integerLog2# :: Integer -> Int# -integerLog2# (S# i) = wordLog2# (int2Word# i) -integerLog2# (J# s ba) = check (s -# 1#) - where - check i = case indexWordArray# ba i of - 0## -> check (i -# 1#) - w -> wordLog2# w +# (uncheckedIShiftL# i WSHIFT#) - --- | This function calculates the integer base 2 logarithm of a 'Word#'. --- @'wordLog2#' 0## = -1#@. -{-# INLINE wordLog2# #-} -wordLog2# :: Word# -> Int# -wordLog2# w = - case leadingZeros of - BA lz -> - let zeros u = indexInt8Array# lz (word2Int# u) in -#if WORD_SIZE_IN_BITS == 64 - case uncheckedShiftRL# w 56# of - a -> - if a `neWord#` 0## - then 64# -# zeros a - else - case uncheckedShiftRL# w 48# of - b -> - if b `neWord#` 0## - then 56# -# zeros b - else - case uncheckedShiftRL# w 40# of - c -> - if c `neWord#` 0## - then 48# -# zeros c - else - case uncheckedShiftRL# w 32# of - d -> - if d `neWord#` 0## - then 40# -# zeros d - else -#endif - case uncheckedShiftRL# w 24# of - e -> - if e `neWord#` 0## - then 32# -# zeros e - else - case uncheckedShiftRL# w 16# of - f -> - if f `neWord#` 0## - then 24# -# zeros f - else - case uncheckedShiftRL# w 8# of - g -> - if g `neWord#` 0## - then 16# -# zeros g - else 8# -# zeros w - --- Lookup table -data BA = BA ByteArray# - -leadingZeros :: BA -leadingZeros = - let mkArr s = - case newByteArray# 256# s of - (# s1, mba #) -> - case writeInt8Array# mba 0# 9# s1 of - s2 -> - let fillA lim val idx st = - if idx ==# 256# - then st - else if idx <# lim - then case writeInt8Array# mba idx val st of - nx -> fillA lim val (idx +# 1#) nx - else fillA (2# *# lim) (val -# 1#) idx st - in case fillA 2# 8# 1# s2 of - s3 -> case unsafeFreezeByteArray# mba s3 of - (# _, ba #) -> ba - in case mkArr realWorld# of - b -> BA b - -#endif diff --git a/Math/NumberTheory/Powers/Cubes.hs b/Math/NumberTheory/Powers/Cubes.hs --- a/Math/NumberTheory/Powers/Cubes.hs +++ b/Math/NumberTheory/Powers/Cubes.hs @@ -24,15 +24,12 @@ import Data.Array.ST import Data.Bits -#if __GLASGOW_HASKELL__ < 705 -import Data.Word -#endif import GHC.Base import GHC.Integer import GHC.Integer.GMP.Internals +import GHC.Integer.Logarithms (integerLog2#) -import Math.NumberTheory.Logarithms.Internal (integerLog2#) import Math.NumberTheory.Unsafe #if __GLASGOW_HASKELL__ < 707 import Math.NumberTheory.Utils (isTrue#) diff --git a/Math/NumberTheory/Powers/Fourth.hs b/Math/NumberTheory/Powers/Fourth.hs --- a/Math/NumberTheory/Powers/Fourth.hs +++ b/Math/NumberTheory/Powers/Fourth.hs @@ -23,16 +23,13 @@ import GHC.Base import GHC.Integer import GHC.Integer.GMP.Internals +import GHC.Integer.Logarithms (integerLog2#) import Data.Array.Unboxed import Data.Array.ST import Data.Bits -#if __GLASGOW_HASKELL__ < 705 -import Data.Word -#endif -import Math.NumberTheory.Logarithms.Internal (integerLog2#) import Math.NumberTheory.Unsafe #if __GLASGOW_HASKELL__ < 707 import Math.NumberTheory.Utils (isTrue#) diff --git a/Math/NumberTheory/Powers/General.hs b/Math/NumberTheory/Powers/General.hs --- a/Math/NumberTheory/Powers/General.hs +++ b/Math/NumberTheory/Powers/General.hs @@ -25,16 +25,13 @@ import GHC.Base import GHC.Integer import GHC.Integer.GMP.Internals +import GHC.Integer.Logarithms (integerLog2#) import Data.Bits -#if __GLASGOW_HASKELL__ < 705 -import Data.Word -#endif import Data.List (foldl') import qualified Data.Set as Set import Math.NumberTheory.Logarithms (integerLogBase') -import Math.NumberTheory.Logarithms.Internal (integerLog2#) import Math.NumberTheory.Utils (shiftToOddCount , splitOff #if __GLASGOW_HASKELL__ < 707 diff --git a/Math/NumberTheory/Powers/Integer.hs b/Math/NumberTheory/Powers/Integer.hs --- a/Math/NumberTheory/Powers/Integer.hs +++ b/Math/NumberTheory/Powers/Integer.hs @@ -16,11 +16,8 @@ ) where import GHC.Base -#if __GLASGOW_HASKELL__ < 705 -import GHC.Word -#endif +import GHC.Integer.Logarithms (wordLog2#) -import Math.NumberTheory.Logarithms.Internal ( wordLog2# ) #if __GLASGOW_HASKELL__ < 707 import Math.NumberTheory.Utils (isTrue#) #endif diff --git a/Math/NumberTheory/Powers/Squares.hs b/Math/NumberTheory/Powers/Squares.hs --- a/Math/NumberTheory/Powers/Squares.hs +++ b/Math/NumberTheory/Powers/Squares.hs @@ -8,11 +8,13 @@ -- -- Functions dealing with squares. Efficient calculation of integer square roots -- and efficient testing for squareness. -{-# LANGUAGE MagicHash, BangPatterns, CPP, FlexibleContexts #-} +{-# LANGUAGE MagicHash, BangPatterns, PatternGuards, CPP, FlexibleContexts #-} module Math.NumberTheory.Powers.Squares ( -- * Square root calculation integerSquareRoot , integerSquareRoot' + , integerSquareRootRem + , integerSquareRootRem' , exactSquareRoot -- * Tests for squares , isSquare @@ -23,24 +25,17 @@ #include "MachDeps.h" -import GHC.Base -import GHC.Integer -import GHC.Integer.GMP.Internals - import Data.Array.Unboxed import Data.Array.ST import Data.Bits -#if __GLASGOW_HASKELL__ < 705 +#if __GLASGOW_HASKELL__ < 709 import Data.Word -- Moved to GHC.Types #endif -import Math.NumberTheory.Logarithms.Internal (integerLog2#) import Math.NumberTheory.Unsafe -#if __GLASGOW_HASKELL__ < 707 -import Math.NumberTheory.Utils (isTrue#) -#endif +import Math.NumberTheory.Powers.Squares.Internal -- | Calculate the integer square root of a nonnegative number @n@, -- that is, the largest integer @r@ with @r*r <= n@. @@ -58,13 +53,40 @@ -- that is, the largest integer @r@ with @r*r <= n@. -- The precondition @n >= 0@ is not checked. {-# RULES -"integerSquareRoot'/Int" integerSquareRoot' = isqrtInt' -"integerSquareRoot'/Word" integerSquareRoot' = isqrtWord +"integerSquareRoot'/Int" integerSquareRoot' = isqrtInt' +"integerSquareRoot'/Word" integerSquareRoot' = isqrtWord +"integerSquareRoot'/Integer" integerSquareRoot' = isqrtInteger #-} {-# INLINE [1] integerSquareRoot' #-} integerSquareRoot' :: Integral a => a -> a integerSquareRoot' = isqrtA +-- | Calculate the integer square root of a nonnegative number as well as +-- the difference of that number with the square of that root, that is if +-- @(s,r) = integerSquareRootRem n@ then @s^2 <= n == s^2+r < (s+1)^2@. +{-# SPECIALISE integerSquareRootRem :: + Int -> (Int, Int), + Word -> (Word, Word), + Integer -> (Integer, Integer) + #-} +integerSquareRootRem :: Integral a => a -> (a, a) +integerSquareRootRem n + | n < 0 = error "integerSquareRootRem: negative argument" + | otherwise = integerSquareRootRem' n + +-- | Calculate the integer square root of a nonnegative number as well as +-- the difference of that number with the square of that root, that is if +-- @(s,r) = integerSquareRootRem' n@ then @s^2 <= n == s^2+r < (s+1)^2@. +-- The precondition @n >= 0@ is not checked. +{-# RULES +"integerSquareRootRem'/Integer" integerSquareRootRem' = karatsubaSqrt + #-} +{-# INLINE [1] integerSquareRootRem' #-} +integerSquareRootRem' :: Integral a => a -> (a, a) +integerSquareRootRem' n = (s, n - s * s) + where + s = integerSquareRoot' n + -- | Returns 'Nothing' if the argument is not a square, -- @'Just' r@ if @r*r == n@ and @r >= 0@. Avoids the expensive calculation -- of the square root if @n@ is recognized as a non-square @@ -77,11 +99,10 @@ #-} exactSquareRoot :: Integral a => a -> Maybe a exactSquareRoot n - | n < 0 = Nothing - | isPossibleSquare n && r*r == n = Just r - | otherwise = Nothing - where - r = integerSquareRoot' n + | n >= 0 + , isPossibleSquare n + , (r, 0) <- integerSquareRootRem' n = Just r + | otherwise = Nothing -- | Test whether the argument is a square. -- After a number is found to be positive, first 'isPossibleSquare' @@ -104,7 +125,10 @@ Integer -> Bool #-} isSquare' :: Integral a => a -> Bool -isSquare' n = isPossibleSquare n && let r = integerSquareRoot' n in r*r == n +isSquare' n + | isPossibleSquare n + , (_, 0) <- integerSquareRootRem' n = True + | otherwise = False -- | Test whether a non-negative number may be a square. -- Non-negativity is not checked, passing negative arguments may @@ -157,63 +181,6 @@ ----------------------------------------------------------------------------- -- Auxiliary Stuff --- Find approximation to square root in 'Integer', then --- find the integer square root by the integer variant --- of Heron's method. Takes only a handful of steps --- unless the input is really large. -{-# SPECIALISE isqrtA :: Integer -> Integer #-} -isqrtA :: Integral a => a -> a -isqrtA 0 = 0 -isqrtA n = heron n (fromInteger . appSqrt . fromIntegral $ n) - --- Heron's method for integers. First make one step to ensure --- the value we're working on is @>= r@, then we have --- @k == r@ iff @k <= step k@. -{-# SPECIALISE heron :: Integer -> Integer -> Integer #-} -heron :: Integral a => a -> a -> a -heron n a = go (step a) - where - step k = (k + n `quot` k) `quot` 2 - go k - | m < k = go m - | otherwise = k - where - m = step k - --- threshold for shifting vs. direct fromInteger --- we shift when we expect more than 256 bits -#if WORD_SIZE_IN_BITS == 64 -#define THRESH 5 -#else -#define THRESH 9 -#endif - --- Find a fairly good approximation to the square root. --- At most one off for small Integers, about 48 bits should be correct --- for large Integers. -appSqrt :: Integer -> Integer -appSqrt (S# i#) = S# (double2Int# (sqrtDouble# (int2Double# i#))) -#if __GLASGOW_HASKELL__ < 709 -appSqrt n@(J# s# _) - | isTrue# (s# <# THRESH#) = floor (sqrt $ fromInteger n :: Double) -#else -appSqrt n@(Jp# bn#) - | isTrue# ((sizeofBigNat# bn#) <# THRESH#) = - floor (sqrt $ fromInteger n :: Double) -#endif - | otherwise = case integerLog2# n of - l# -> case uncheckedIShiftRA# l# 1# -# 47# of - h# -> case shiftRInteger n (2# *# h#) of - m -> case floor (sqrt $ fromInteger m :: Double) of - r -> shiftLInteger r h# -#if __GLASGOW_HASKELL__ >= 709 --- There's already a check for negative in integerSquareRoot, --- but integerSquareRoot' is exported directly too. -appSqrt _ = error "integerSquareRoot': negative argument" -#endif - --- Auxiliaries - -- Make an array indicating whether a remainder is a square remainder. sqRemArray :: Int -> UArray Int Bool sqRemArray md = runSTUArray $ do @@ -250,7 +217,7 @@ sr325 :: UArray Int Bool sr325 = sqRemArray 325 --- Specialisations for Int and Word +-- Specialisations for Int, Word, and Integer -- For @n <= 2^64@, the result of -- @@ -286,3 +253,6 @@ where !r = (fromIntegral :: Int -> Word) . (truncate :: Double -> Int) . sqrt $ fromIntegral n +{-# INLINE isqrtInteger #-} +isqrtInteger :: Integer -> Integer +isqrtInteger = fst . karatsubaSqrt diff --git a/Math/NumberTheory/Powers/Squares/Internal.hs b/Math/NumberTheory/Powers/Squares/Internal.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Powers/Squares/Internal.hs @@ -0,0 +1,144 @@ +-- | +-- Module: Math.NumberTheory.Powers.Squares.Internal +-- Copyright: (c) 2016 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- Stability: Provisional +-- Portability: Non-portable (GHC extensions) +-- +-- Internal functions dealing with square roots. End-users should not import this module. + +{-# LANGUAGE MagicHash #-} +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE PatternGuards #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE FlexibleContexts #-} + +module Math.NumberTheory.Powers.Squares.Internal + ( karatsubaSqrt + , isqrtA + ) where + +#include "MachDeps.h" + +import Data.Bits + +import GHC.Base +import GHC.Integer +import GHC.Integer.GMP.Internals +import GHC.Integer.Logarithms (integerLog2#) + +import Math.NumberTheory.Logarithms (integerLog2) +#if __GLASGOW_HASKELL__ < 707 +import Math.NumberTheory.Utils (isTrue#) +#endif + +-- Find approximation to square root in 'Integer', then +-- find the integer square root by the integer variant +-- of Heron's method. Takes only a handful of steps +-- unless the input is really large. +{-# SPECIALISE isqrtA :: Integer -> Integer #-} +isqrtA :: Integral a => a -> a +isqrtA 0 = 0 +isqrtA n = heron n (fromInteger . appSqrt . fromIntegral $ n) + +-- Heron's method for integers. First make one step to ensure +-- the value we're working on is @>= r@, then we have +-- @k == r@ iff @k <= step k@. +{-# SPECIALISE heron :: Integer -> Integer -> Integer #-} +heron :: Integral a => a -> a -> a +heron n a = go (step a) + where + step k = (k + n `quot` k) `quot` 2 + go k + | m < k = go m + | otherwise = k + where + m = step k + +-- threshold for shifting vs. direct fromInteger +-- we shift when we expect more than 256 bits +#if WORD_SIZE_IN_BITS == 64 +#define THRESH 5 +#else +#define THRESH 9 +#endif + +-- Find a fairly good approximation to the square root. +-- At most one off for small Integers, about 48 bits should be correct +-- for large Integers. +appSqrt :: Integer -> Integer +appSqrt (S# i#) = S# (double2Int# (sqrtDouble# (int2Double# i#))) +#if __GLASGOW_HASKELL__ < 709 +appSqrt n@(J# s# _) + | isTrue# (s# <# THRESH#) = floor (sqrt $ fromInteger n :: Double) +#else +appSqrt n@(Jp# bn#) + | isTrue# ((sizeofBigNat# bn#) <# THRESH#) = + floor (sqrt $ fromInteger n :: Double) +#endif + | otherwise = case integerLog2# n of + l# -> case uncheckedIShiftRA# l# 1# -# 47# of + h# -> case shiftRInteger n (2# *# h#) of + m -> case floor (sqrt $ fromInteger m :: Double) of + r -> shiftLInteger r h# +#if __GLASGOW_HASKELL__ >= 709 +-- There's already a check for negative in integerSquareRoot, +-- but integerSquareRoot' is exported directly too. +appSqrt _ = error "integerSquareRoot': negative argument" +#endif + + +-- Integer square root with remainder, using the Karatsuba Square Root +-- algorithm from +-- Paul Zimmermann. Karatsuba Square Root. [Research Report] RR-3805, 1999, +-- pp.8. + +karatsubaSqrt :: Integer -> (Integer, Integer) +karatsubaSqrt 0 = (0, 0) +karatsubaSqrt n + | lgN < 2300 = + let s = isqrtA n in (s, n - s * s) + | otherwise = + if lgN .&. 2 /= 0 then + karatsubaStep k (karatsubaSplit k n) + else + -- before we split n into 4 part we must ensure that the first part + -- is at least 2^k/4, since this doesn't happen here we scale n by + -- multiplying it by 4 + let n' = n `unsafeShiftL` 2 + (s, r) = karatsubaStep k (karatsubaSplit k n') + r' | s .&. 1 == 0 = r + | otherwise = r + double s - 1 + in (s `unsafeShiftR` 1, r' `unsafeShiftR` 2) + where + k = lgN `unsafeShiftR` 2 + 1 + lgN = integerLog2 n + +karatsubaStep :: Int -> (Integer, Integer, Integer, Integer) -> (Integer, Integer) +karatsubaStep k (a3, a2, a1, a0) + | r >= 0 = (s, r) + | otherwise = (s - 1, r + double s - 1) + where + r = cat u a0 - q * q + s = s' `unsafeShiftL` k + q + (q, u) = cat r' a1 `quotRem` double s' + (s', r') = karatsubaSqrt (cat a3 a2) + cat x y = x `unsafeShiftL` k .|. y + {-# INLINE cat #-} + +karatsubaSplit :: Int -> Integer -> (Integer, Integer, Integer, Integer) +karatsubaSplit k n0 = (a3, a2, a1, a0) + where + a3 = n3 + n3 = n2 `unsafeShiftR` k + a2 = n2 .&. m + n2 = n1 `unsafeShiftR` k + a1 = n1 .&. m + n1 = n0 `unsafeShiftR` k + a0 = n0 .&. m + m = 1 `unsafeShiftL` k - 1 + +double :: Bits a => a -> a +double x = x `unsafeShiftL` 1 +{-# INLINE double #-} diff --git a/Math/NumberTheory/Primes/Counting/Impl.hs b/Math/NumberTheory/Primes/Counting/Impl.hs --- a/Math/NumberTheory/Primes/Counting/Impl.hs +++ b/Math/NumberTheory/Primes/Counting/Impl.hs @@ -9,9 +9,7 @@ -- Number of primes not exceeding @n@, @π(n)@, and @n@-th prime. -- {-# LANGUAGE CPP, BangPatterns, FlexibleContexts #-} -#if __GLASGOW_HASKELL__ >= 700 {-# OPTIONS_GHC -fspec-constr-count=24 #-} -#endif {-# OPTIONS_HADDOCK hide #-} module Math.NumberTheory.Primes.Counting.Impl ( primeCount @@ -31,9 +29,6 @@ import Math.NumberTheory.Unsafe import Data.Array.ST -#if !MIN_VERSION_array(0,5,0) - hiding (unsafeThaw) -#endif import Control.Monad.ST import Data.Bits import Data.Int diff --git a/Math/NumberTheory/Primes/Factorisation.hs b/Math/NumberTheory/Primes/Factorisation.hs --- a/Math/NumberTheory/Primes/Factorisation.hs +++ b/Math/NumberTheory/Primes/Factorisation.hs @@ -13,6 +13,9 @@ -- on the canonical factorisation, these require that the number be positive -- and in the case of the Carmichael function that the list of prime factors -- with their multiplicities is ascending. + +{-# OPTIONS_GHC -fno-warn-warnings-deprecations #-} + module Math.NumberTheory.Primes.Factorisation ( -- * Factorisation functions -- $algorithm @@ -73,6 +76,8 @@ import Math.NumberTheory.Primes.Factorisation.Montgomery import Math.NumberTheory.Primes.Factorisation.TrialDivision import Math.NumberTheory.Primes.Sieve.Misc + +{-# DEPRECATED totient, φ, carmichael, λ, moebius, μ, divisors, tau, τ, divisorCount, divisorSum, sigma, σ, divisorPowerSum "Use 'Math.NumberTheory.ArithmeticFunctions'" #-} -- $algorithm -- 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 @@ -43,9 +43,6 @@ #include "MachDeps.h" import GHC.Base -#if __GLASGOW_HASKELL__ < 705 -import GHC.Word -- Moved to GHC.Types -#endif import System.Random import Control.Monad.State.Strict @@ -55,8 +52,9 @@ import Data.Bits import Data.Maybe +import GHC.Integer.Logarithms + import Math.NumberTheory.Logarithms -import Math.NumberTheory.Logarithms.Internal import Math.NumberTheory.Powers.General (highestPower, largePFPower) import Math.NumberTheory.Powers.Squares (integerSquareRoot') import Math.NumberTheory.Primes.Sieve.Eratosthenes diff --git a/Math/NumberTheory/Primes/Factorisation/Utils.hs b/Math/NumberTheory/Primes/Factorisation/Utils.hs --- a/Math/NumberTheory/Primes/Factorisation/Utils.hs +++ b/Math/NumberTheory/Primes/Factorisation/Utils.hs @@ -27,6 +27,8 @@ import Math.NumberTheory.Powers.Integer +{-# DEPRECATED totientFromCanonical, carmichaelFromCanonical, moebiusFromCanonical, divisorsFromCanonical, tauFromCanonical, divisorSumFromCanonical, sigmaFromCanonical "Use 'Math.NumberTheory.ArithmeticFunctions'" #-} + -- | Totient of a prime power. ppTotient :: (Integer,Int) -> Integer ppTotient (p,1) = p-1 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 @@ -16,9 +16,7 @@ -- This module is mainly intended for comparison and verification. {-# LANGUAGE BangPatterns, CPP, MonoLocalBinds #-} {-# OPTIONS_GHC -funbox-strict-fields #-} -#if __GLASGOW_HASKELL__ >= 700 {-# OPTIONS_GHC -fno-float-in -fno-spec-constr -fno-full-laziness #-} -#endif module Math.NumberTheory.Primes.Heap (primes, sieveFrom) where import Data.Array.Unboxed @@ -58,10 +56,6 @@ {-# SPECIALISE push :: Word -> Word -> Int -> Hipp Word -> Hipp Word #-} {-# SPECIALISE push :: Integer -> Integer -> Int -> Hipp Integer -> Hipp Integer #-} push :: Integral a => a -> a -> Int -> Hipp a -> Hipp a --- GHC 7 does not like the old code, so it gets a new implementation. --- That is faster than what it does with the old code, but still slower --- than what GHC 6 did with it. :( -#if __GLASGOW_HASKELL__ >= 700 push !c !p !w = go where less = (< c) @@ -69,21 +63,12 @@ | less hc = H hc hp hw (go r) l | otherwise = H c p w (push hc hp hw r) l go _ = H c p w E E -#else -push c p w (H hc hp hw l r) - | c < hc = H c p w (push hc hp hw r) l - | otherwise = H hc hp hw (push c p w r) l -push c p w E = H c p w E E -#endif -- bubble down increased top to regain heap invariant {-# SPECIALISE bubble :: Hipp Int -> Hipp Int #-} {-# SPECIALISE bubble :: Hipp Word -> Hipp Word #-} {-# SPECIALISE bubble :: Hipp Integer -> Hipp Integer #-} bubble :: Integral a => Hipp a -> Hipp a --- Again, GHC 6 fared better, so new code for GHC 7, still --- not as good as GHC 6 was. -#if __GLASGOW_HASKELL__ >= 700 bubble h@(H c p w l r) = case r of E -> case l of @@ -100,19 +85,9 @@ | rc < c -> H rc rp rw l (mkHipp c p w rl rr) | otherwise -> h _ -> error "Heap invariant violated, left smaller than right!" -#else -bubble h@(H c p w l@(H lc lp lw ll lr) r@(H rc rp rw rl rr)) - | c <= lc && c <= rc = h - | lc < rc = H lc lp lw (bubble (H c p w ll lr)) r - | otherwise = H rc rp rw l (bubble (H c p w rl rr)) -bubble h@(H c p w (H lc lp lw _ _) _) - | c <= lc = h - | otherwise = H lc lp lw (H c p w E E) E -#endif bubble h = h --- join two heaps and composite-data, GHC 7 doesn't do well on the old bubble. -#if __GLASGOW_HASKELL__ >= 700 +-- join two heaps and composite-data {-# SPECIALISE mkHipp :: Int -> Int -> Int -> Hipp Int -> Hipp Int -> Hipp Int, Integer -> Integer -> Int -> Hipp Integer -> Hipp Integer -> Hipp Integer, @@ -138,8 +113,6 @@ | less rc -> H rc rp rw l (go rl rr) | otherwise -> H c p w l r _ -> error "Heap invariant violated, left smaller than right!" --- {-# INLINE mkHipp #-} -#endif -- increase the top of the heap and re-heap {-# SPECIALISE inc :: Hipp Int -> Hipp Int #-} diff --git a/Math/NumberTheory/Primes/Sieve/Eratosthenes.hs b/Math/NumberTheory/Primes/Sieve/Eratosthenes.hs --- a/Math/NumberTheory/Primes/Sieve/Eratosthenes.hs +++ b/Math/NumberTheory/Primes/Sieve/Eratosthenes.hs @@ -9,9 +9,7 @@ -- Sieve -- {-# LANGUAGE CPP, BangPatterns, FlexibleContexts #-} -#if __GLASGOW_HASKELL__ >= 700 {-# OPTIONS_GHC -fspec-constr-count=8 #-} -#endif {-# OPTIONS_HADDOCK hide #-} module Math.NumberTheory.Primes.Sieve.Eratosthenes ( primes @@ -36,9 +34,6 @@ import Control.Monad.ST import Data.Array.ST -#if !MIN_VERSION_array(0,5,0) - hiding (unsafeFreeze, unsafeThaw, castSTUArray) -#endif import Control.Monad (when) import Data.Bits #if __GLASGOW_HASKELL__ < 709 || WORD_SIZE_IN_BITS == 32 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 @@ -6,10 +6,12 @@ -- Stability: Provisional -- Portability: Non-portable (GHC extensions) -- -{-# LANGUAGE CPP, BangPatterns, ScopedTypeVariables, MonoLocalBinds, FlexibleContexts #-} -#if __GLASGOW_HASKELL__ >= 700 +{-# LANGUAGE BangPatterns #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE MonoLocalBinds #-} +{-# LANGUAGE ScopedTypeVariables #-} {-# OPTIONS_GHC -fspec-constr-count=8 #-} -#endif +{-# OPTIONS_GHC -fno-warn-warnings-deprecations #-} {-# OPTIONS_HADDOCK hide #-} module Math.NumberTheory.Primes.Sieve.Misc ( -- * Types @@ -46,6 +48,27 @@ import Math.NumberTheory.Unsafe import Math.NumberTheory.Utils +{- +IMPORTANT NOTICE: Not all sieves use the same layout! + +FactorSieve: + + To remain as efficient as possible, FactorSieve omits only even numbers. + To relate an odd number x to its index i: + + i = (x `div` 2) - 1 + x = i * 2 + 3 + +TotientSieve, CarmichaelSieve: + + These sieves use a (2,3,5) wheel optimization, sacrificing performance to save + more memory. The only indices stored are those coprime to 2, 3, and 5. + To relate such an integer x to its index i: + + i = toIdx x + x = toPrim i +-} + -- | A compact store of smallest prime factors. data FactorSieve = FS {-# UNPACK #-} !Word {-# UNPACK #-} !(UArray Int Word16) @@ -148,8 +171,8 @@ j | j <= bound -> intLoop (fromIntegral (j `shiftR` 1) - 1) | otherwise -> tdLoop j (integerSquareRoot' j) (ix+1) where - p = toPrim ix - pix = unsafeAt sve $ fromIntegral p + p = fromIntegral $ 2 * ix + 3 + pix = unsafeAt sve ix curve n = stdGenFactorisation (Just (bound*(bound+2))) (mkStdGen $ fromIntegral n `xor` 0xdecaf00d) Nothing n -- | @'totientSieve' n@ creates a store of the totients of the numbers not exceeding @n@. @@ -277,6 +300,9 @@ curve tt n = tt `lcm` carmichaelFromCanonical (stdGenFactorisation (Just (bound*(bound+2))) (mkStdGen $ fromIntegral n `xor` 0xdecaf00d) Nothing n) +-- NOTE: This is a legacy implementation of FactorSieve which uses the +-- same (2,3,5) wheel optimization as the other sieves. +-- It is still used to generate the other sieves. spfSieve :: Word -> ST s (STUArray s Int Word) spfSieve bound = do let (octs,lidx) = idxPr bound diff --git a/Math/NumberTheory/Primes/Testing/Probabilistic.hs b/Math/NumberTheory/Primes/Testing/Probabilistic.hs --- a/Math/NumberTheory/Primes/Testing/Probabilistic.hs +++ b/Math/NumberTheory/Primes/Testing/Probabilistic.hs @@ -28,9 +28,6 @@ import GHC.Base -#if __GLASGOW_HASKELL__ < 705 -import GHC.Word -- Moved to GHC.Types -#endif import GHC.Integer.GMP.Internals -- | @'isPrime' n@ tests whether @n@ is a prime (negative or positive). diff --git a/Math/NumberTheory/UniqueFactorisation.hs b/Math/NumberTheory/UniqueFactorisation.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/UniqueFactorisation.hs @@ -0,0 +1,131 @@ +-- | +-- Module: Math.NumberTheory.UniqueFactorisation +-- Copyright: (c) 2016 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- Stability: Provisional +-- Portability: Non-portable (GHC extensions) +-- +-- An abstract type class for unique factorisation domains. +-- + +{-# LANGUAGE CPP #-} +{-# LANGUAGE TypeFamilies #-} + +module Math.NumberTheory.UniqueFactorisation + ( Prime + , UniqueFactorisation(..) + ) where + +import Control.Arrow + +#if MIN_VERSION_base(4,8,0) +#else +import Data.Word +#endif + +import Math.NumberTheory.Primes.Factorisation as F (factorise') +import Math.NumberTheory.GaussianIntegers as G + +import Numeric.Natural + +#if MIN_VERSION_base(4,7,0) +import Data.Coerce +#else +import Unsafe.Coerce + +coerce :: a -> b +coerce = unsafeCoerce +#endif + +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@: +-- +-- > abs n == abs (product (map (\(p, k) -> unPrime p ^ k) (factorise n))) +-- +-- The result of 'factorise' should not contain zero powers and should not change after multiplication of the argument by domain's unit. +class UniqueFactorisation a where + unPrime :: Prime a -> a + factorise :: a -> [(Prime a, Word)] + +instance UniqueFactorisation Int where + unPrime = coerce wordToInt + factorise m' = if m <= 1 + then [] + else map (coerce integerToWord *** intToWord) . F.factorise' . intToInteger $ m + where + m = abs m' + +instance UniqueFactorisation Word where + unPrime = coerce + factorise m = if m <= 1 + then [] + else map (coerce integerToWord *** intToWord) . F.factorise' . wordToInteger $ m + + +instance UniqueFactorisation Integer where + unPrime = coerce naturalToInteger + factorise m' = if m <= 1 + then [] + else map (coerce integerToNatural *** intToWord) . F.factorise' $ m + where + m = abs m' + +instance UniqueFactorisation Natural where + unPrime = coerce + factorise m = if m <= 1 + then [] + else map (coerce integerToNatural *** intToWord) . F.factorise' . naturalToInteger $ m + +newtype GaussianPrime = GaussianPrime { _unGaussianPrime :: G.GaussianInteger } + deriving (Eq, Show) + +instance UniqueFactorisation G.GaussianInteger where + unPrime = coerce + + 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/Unsafe.hs b/Math/NumberTheory/Unsafe.hs --- a/Math/NumberTheory/Unsafe.hs +++ b/Math/NumberTheory/Unsafe.hs @@ -34,9 +34,6 @@ , (!) ) import Data.Array.MArray -#if !MIN_VERSION_array(0,5,0) - hiding (unsafeFreeze, unsafeThaw) -#endif unsafeAt :: (IArray a e, Ix i) => a i e -> i -> e unsafeAt = (!) diff --git a/Math/NumberTheory/Utils.hs b/Math/NumberTheory/Utils.hs --- a/Math/NumberTheory/Utils.hs +++ b/Math/NumberTheory/Utils.hs @@ -28,9 +28,6 @@ #include "MachDeps.h" import GHC.Base -#if __GLASGOW_HASKELL__ < 705 -import GHC.Word -- Word and its constructor moved to GHC.Types -#endif import GHC.Integer import GHC.Integer.GMP.Internals @@ -178,24 +175,11 @@ -- | Number of 1-bits in a @'Word'@. bitCountWord :: Word -> Int -#if __GLASGOW_HASKELL__ >= 703 bitCountWord = popCount --- should yield a machine instruction -#else -bitCountWord w = case w - (shiftR w 1 .&. m5) of - !w1 -> case (w1 .&. m3) + (shiftR w1 2 .&. m3) of - !w2 -> case (w2 + shiftR w2 4) .&. mf of - !w3 -> fromIntegral (shiftR (w3 * m1) sd) -#endif -- | Number of 1-bits in an @'Int'@. bitCountInt :: Int -> Int -#if __GLASGOW_HASKELL__ >= 703 bitCountInt = popCount --- should yield a machine instruction -#else -bitCountInt (I# i#) = bitCountWord (W# (int2Word# i#)) -#endif -- | Number of trailing zeros in a @'Word#'@, wrong for @0@. {-# INLINE trailZeros# #-} @@ -214,11 +198,7 @@ -- Int -> Int -> (Int, Int), -- Word -> Word -> (Int, Word) -- #-} -#if __GLASGOW_HASKELL__ >= 700 {-# INLINABLE splitOff #-} -#else -{-# INLINE splitOff #-} -#endif splitOff :: Integral a => a -> a -> (Int, a) splitOff p n = go 0 n where diff --git a/arithmoi.cabal b/arithmoi.cabal --- a/arithmoi.cabal +++ b/arithmoi.cabal @@ -1,5 +1,5 @@ name : arithmoi -version : 0.4.2.0 +version : 0.4.3.0 cabal-version : >= 1.10 author : Daniel Fischer copyright : (c) 2011 Daniel Fischer @@ -21,10 +21,6 @@ powers (integer roots and tests, modular exponentiation), integer logarithms. - Note: Requires GHC >= 6.12 with the integer-gmp package - for efficiency. Portability is on the to-do list (with - low priority, however). - category : Math, Algorithms, Number Theory tested-with : GHC==7.6.3, GHC==7.8.4, GHC==7.10.3, GHC==8.0.1 @@ -38,15 +34,22 @@ library default-language: Haskell2010 - build-depends : base >= 4 && < 5 - , array >= 0.3 && < 0.6 - , ghc-prim < 0.6 - , integer-gmp < 1.1 - , containers >= 0.3 && < 0.6 - , random >= 1.0 && < 1.2 - , mtl >= 2.0 && < 2.3 + build-depends : base >= 4.6 && < 5 + , array >= 0.5 && < 0.6 + , ghc-prim < 0.6 + , integer-gmp < 1.1 + , containers >= 0.5 && < 0.6 + , random >= 1.0 && < 1.2 + , mtl >= 2.0 && < 2.3 + if impl(ghc < 7.10) + build-depends : nats >= 1 && <1.2 + if impl(ghc < 8.0) + build-depends : semigroups >= 0.8 - exposed-modules : Math.NumberTheory.Logarithms + exposed-modules : Math.NumberTheory.ArithmeticFunctions + Math.NumberTheory.ArithmeticFunctions.Class + Math.NumberTheory.ArithmeticFunctions.Standard + Math.NumberTheory.Logarithms Math.NumberTheory.Moduli Math.NumberTheory.MoebiusInversion Math.NumberTheory.MoebiusInversion.Int @@ -56,6 +59,7 @@ Math.NumberTheory.GCD.LowLevel Math.NumberTheory.Powers Math.NumberTheory.Powers.Squares + Math.NumberTheory.Powers.Squares.Internal Math.NumberTheory.Powers.Cubes Math.NumberTheory.Powers.Fourth Math.NumberTheory.Powers.General @@ -68,9 +72,9 @@ Math.NumberTheory.Primes.Testing Math.NumberTheory.Primes.Testing.Certificates Math.NumberTheory.Primes.Heap + Math.NumberTheory.UniqueFactorisation other-modules : Math.NumberTheory.Utils Math.NumberTheory.Unsafe - Math.NumberTheory.Logarithms.Internal Math.NumberTheory.Primes.Counting.Impl Math.NumberTheory.Primes.Counting.Approximate Math.NumberTheory.Primes.Factorisation.Montgomery @@ -96,8 +100,12 @@ benchmark criterion build-depends: base - ,arithmoi - ,criterion + , arithmoi + , criterion + , containers + , random + other-modules: Math.NumberTheory.ArithmeticFunctionsBench + , Math.NumberTheory.PowersBench hs-source-dirs: benchmark main-is: Bench.hs type: exitcode-stdio-1.0 @@ -109,15 +117,22 @@ ghc-options: -Wall main-is: Test.hs default-language: Haskell2010 - build-depends: base >= 4 && < 5 - , arithmoi >= 0.4 && < 0.5 - , tasty >= 0.10 && < 0.12 - , tasty-smallcheck >= 0.8 && < 0.9 - , tasty-quickcheck >= 0.8 && < 0.9 - , tasty-hunit >= 0.9 && < 0.10 - , QuickCheck >= 2.8 && < 2.9 - , smallcheck >= 1.1 && < 1.2 - other-modules : Math.NumberTheory.GaussianIntegersTests + build-depends: base >= 4.6 && < 5 + , containers >= 0.5 && < 0.6 + , arithmoi >= 0.4 && < 0.5 + , tasty >= 0.10 && < 0.12 + , tasty-smallcheck >= 0.8 && < 0.9 + , tasty-quickcheck >= 0.8 && < 0.9 + , tasty-hunit >= 0.9 && < 0.10 + , QuickCheck >= 2.7.6 && < 2.10 + , smallcheck >= 1.1 && < 1.2 + , transformers >= 0.3 + , transformers-compat >= 0.4 + if impl(ghc < 7.10) + build-depends : nats >= 1 && <1.2 + + other-modules : Math.NumberTheory.ArithmeticFunctionsTests + , Math.NumberTheory.GaussianIntegersTests , Math.NumberTheory.GCDTests , Math.NumberTheory.GCD.LowLevelTests , Math.NumberTheory.LogarithmsTests @@ -135,3 +150,6 @@ , Math.NumberTheory.Primes.HeapTests , Math.NumberTheory.Primes.SieveTests , Math.NumberTheory.TestUtils + , Math.NumberTheory.TestUtils.Wrappers + , Math.NumberTheory.TestUtils.Compose + , Math.NumberTheory.UniqueFactorisationTests diff --git a/benchmark/Bench.hs b/benchmark/Bench.hs --- a/benchmark/Bench.hs +++ b/benchmark/Bench.hs @@ -1,3 +1,13 @@ +{-# OPTIONS_GHC -fno-warn-warnings-deprecations #-} + module Main where -main = return () +import Criterion.Main + +import Math.NumberTheory.ArithmeticFunctionsBench as ArithmeticFunctions +import Math.NumberTheory.PowersBench as Powers + +main = defaultMain + [ ArithmeticFunctions.benchSuite + , Powers.benchSuite + ] diff --git a/benchmark/Math/NumberTheory/ArithmeticFunctionsBench.hs b/benchmark/Math/NumberTheory/ArithmeticFunctionsBench.hs new file mode 100644 --- /dev/null +++ b/benchmark/Math/NumberTheory/ArithmeticFunctionsBench.hs @@ -0,0 +1,34 @@ +{-# OPTIONS_GHC -fno-warn-warnings-deprecations #-} + +module Math.NumberTheory.ArithmeticFunctionsBench + ( benchSuite + ) where + +import Criterion.Main +import Data.Set (Set) + +import Math.NumberTheory.ArithmeticFunctions as A +import Math.NumberTheory.Primes.Factorisation as F + +compareFunctions :: String -> (Integer -> Integer) -> (Integer -> Integer) -> Benchmark +compareFunctions name old new = bgroup name + [ bench "old" $ nf (map old) [1..100000] + , bench "new" $ nf (map new) [1..100000] + ] + +compareSetFunctions :: String -> (Integer -> Set Integer) -> (Integer -> Set Integer) -> Benchmark +compareSetFunctions name old new = bgroup name + [ bench "old" $ nf (map old) [1..100000] + , bench "new" $ nf (map new) [1..100000] + ] + +benchSuite = bgroup "ArithmeticFunctions" + [ compareSetFunctions "divisors" F.divisors A.divisors + , bench "divisors/int" $ nf (map A.divisorsSmall) [1 :: Int .. 100000] + , compareFunctions "totient" F.totient A.totient + , compareFunctions "carmichael" F.carmichael A.carmichael + , compareFunctions "moebius" F.moebius A.moebius + , compareFunctions "tau" F.tau A.tau + , compareFunctions "sigma 1" (F.sigma 1) (A.sigma 1) + , compareFunctions "sigma 2" (F.sigma 2) (A.sigma 2) + ] diff --git a/benchmark/Math/NumberTheory/PowersBench.hs b/benchmark/Math/NumberTheory/PowersBench.hs new file mode 100644 --- /dev/null +++ b/benchmark/Math/NumberTheory/PowersBench.hs @@ -0,0 +1,30 @@ +{-# OPTIONS_GHC -fno-warn-warnings-deprecations #-} + +module Math.NumberTheory.PowersBench + ( benchSuite + ) where + +import Criterion.Main +import System.Random + +import Math.NumberTheory.Logarithms (integerLog2) +import Math.NumberTheory.Powers.Squares.Internal + +genInteger :: Int -> Int -> Integer +genInteger salt bits + = head + . dropWhile ((< bits) . integerLog2) + . scanl (\a r -> a * 2^31 + abs r) 1 + . randoms + . mkStdGen + $ salt + bits + +compareRoots :: Int -> Benchmark +compareRoots bits = bgroup ("sqrt" ++ show bits) + [ bench "new" $ nf (fst . karatsubaSqrt) n + , bench "old" $ nf isqrtA n + ] + where + n = genInteger 0 bits + +benchSuite = bgroup "Powers" $ map compareRoots [2300, 2400 .. 2600] diff --git a/test-suite/Math/NumberTheory/ArithmeticFunctionsTests.hs b/test-suite/Math/NumberTheory/ArithmeticFunctionsTests.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/NumberTheory/ArithmeticFunctionsTests.hs @@ -0,0 +1,279 @@ +-- | +-- Module: Math.NumberTheory.ArithmeticFunctionsTests +-- Copyright: (c) 2016 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- Stability: Provisional +-- +-- Tests for Math.NumberTheory.ArithmeticFunctions +-- + +{-# LANGUAGE CPP #-} + +{-# OPTIONS_GHC -fno-warn-type-defaults #-} + +module Math.NumberTheory.ArithmeticFunctionsTests + ( testSuite + ) where + +import Test.Tasty +import Test.Tasty.HUnit + +#if MIN_VERSION_base(4,8,0) +#else +import Prelude hiding (sum, all, elem) +import Data.Foldable +#endif + +import qualified Data.Set as S +import qualified Data.IntSet as IS + +import Math.NumberTheory.ArithmeticFunctions +import Math.NumberTheory.Primes.Factorisation +import Math.NumberTheory.TestUtils + +import Numeric.Natural + +oeisAssertion :: (Eq a, Show a) => String -> ArithmeticFunction Natural a -> [a] -> Assertion +oeisAssertion name f baseline = assertEqual name baseline (map (runFunction f) [1 .. fromIntegral (length baseline)]) + +-- | tau(n) equals to a number of divisors. +divisorsProperty1 :: Natural -> Bool +divisorsProperty1 n = S.size (runFunction divisorsA n) == runFunction tauA n + +-- | sigma(n) equals to a number of divisors. +divisorsProperty2 :: Natural -> Bool +divisorsProperty2 n = sum (runFunction divisorsA n) == runFunction (sigmaA 1) n + +-- | All divisors of n truly divides n. +divisorsProperty3 :: Natural -> Bool +divisorsProperty3 n = all (\d -> n `mod` d == 0) (runFunction divisorsA n) + +-- | All divisors of n truly divides n. +divisorsProperty4 :: Int -> Bool +divisorsProperty4 n = S.toAscList (runFunction divisorsA n) == IS.toAscList (runFunction divisorsSmallA n) + +-- | tau matches baseline from OEIS. +tauOeis :: Assertion +tauOeis = oeisAssertion "A000005" tauA + [ 1, 2, 2, 3, 2, 4, 2, 4, 3, 4, 2, 6, 2, 4, 4, 5, 2, 6, 2, 6, 4, 4, 2, 8 + , 3, 4, 4, 6, 2, 8, 2, 6, 4, 4, 4, 9, 2, 4, 4, 8, 2, 8, 2, 6, 6, 4, 2, 10 + , 3, 6, 4, 6, 2, 8, 4, 8, 4, 4, 2, 12, 2, 4, 6, 7, 4, 8, 2, 6, 4, 8, 2 + , 12, 2, 4, 6, 6, 4, 8, 2, 10, 5, 4, 2, 12, 4, 4, 4, 8, 2, 12, 4, 6, 4, 4 + , 4, 12, 2, 6, 6, 9, 2, 8, 2, 8 + ] + +-- | sigma_0 coincides with tau by definition +sigmaProperty1 :: Natural -> Bool +sigmaProperty1 n = runFunction tauA n == runFunction (sigmaA 0) n + +-- | value of totient is bigger than argument +sigmaProperty2 :: Natural -> Bool +sigmaProperty2 n = n <= 1 || runFunction (sigmaA 1) n > n + +-- | sigma_1 matches baseline from OEIS. +sigma1Oeis :: Assertion +sigma1Oeis = oeisAssertion "A000203" (sigmaA 1) + [ 1, 3, 4, 7, 6, 12, 8, 15, 13, 18, 12, 28, 14, 24, 24, 31, 18, 39, 20 + , 42, 32, 36, 24, 60, 31, 42, 40, 56, 30, 72, 32, 63, 48, 54, 48, 91, 38 + , 60, 56, 90, 42, 96, 44, 84, 78, 72, 48, 124, 57, 93, 72, 98, 54, 120 + , 72, 120, 80, 90, 60, 168, 62, 96, 104, 127, 84, 144, 68, 126, 96, 144 + ] + +-- | sigma_2 matches baseline from OEIS. +sigma2Oeis :: Assertion +sigma2Oeis = oeisAssertion "A001157" (sigmaA 2) + [ 1, 5, 10, 21, 26, 50, 50, 85, 91, 130, 122, 210, 170, 250, 260, 341, 290 + , 455, 362, 546, 500, 610, 530, 850, 651, 850, 820, 1050, 842, 1300, 962 + , 1365, 1220, 1450, 1300, 1911, 1370, 1810, 1700, 2210, 1682, 2500, 1850 + , 2562, 2366, 2650, 2210, 3410, 2451, 3255 + ] + +-- | value of totient if even, except totient(1) and totient(2) +totientProperty1 :: Natural -> Bool +totientProperty1 n = n <= 2 || even (runFunction totientA n) + +-- | value of totient is smaller than argument +totientProperty2 :: Natural -> Bool +totientProperty2 n = n <= 1 || runFunction totientA n < n + +totientSieve100 :: TotientSieve +totientSieve100 = totientSieve 100 + +-- | totient matches sieveTotient +totientProperty3 :: Natural -> Bool +totientProperty3 n = n < 1 + || fromIntegral (runFunction totientA n) + == sieveTotient totientSieve100 (fromIntegral n) + +-- | totient matches baseline from OEIS. +totientOeis :: Assertion +totientOeis = oeisAssertion "A000010" totientA + [ 1, 1, 2, 2, 4, 2, 6, 4, 6, 4, 10, 4, 12, 6, 8, 8, 16, 6, 18, 8, 12, 10 + , 22, 8, 20, 12, 18, 12, 28, 8, 30, 16, 20, 16, 24, 12, 36, 18, 24, 16, 40 + , 12, 42, 20, 24, 22, 46, 16, 42, 20, 32, 24, 52, 18, 40, 24, 36, 28, 58 + , 16, 60, 30, 36, 32, 48, 20, 66, 32, 44 + ] + +-- | jordan_0 is zero for argument > 1 +jordanProperty1 :: Natural -> Bool +jordanProperty1 n = n <= 1 || runFunction (jordanA 0) n == 0 + +-- | jordan_1 coincides with totient by definition +jordanProperty2 :: Natural -> Bool +jordanProperty2 n = runFunction totientA n == runFunction (jordanA 1) n + +-- | jordan_2 matches baseline from OEIS. +jordan2Oeis :: Assertion +jordan2Oeis = oeisAssertion "A007434" (jordanA 2) + [ 1, 3, 8, 12, 24, 24, 48, 48, 72, 72, 120, 96, 168, 144, 192, 192, 288 + , 216, 360, 288, 384, 360, 528, 384, 600, 504, 648, 576, 840, 576, 960 + , 768, 960, 864, 1152, 864, 1368, 1080, 1344, 1152, 1680, 1152, 1848, 1440 + , 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)) + +-- | 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 + ] + +-- | liouville values are [-1, 1] +liouvilleProperty1 :: Natural -> Bool +liouvilleProperty1 n = runFunction liouvilleA n `elem` [-1, 1] + +-- | moebius is zero or equal to liouville +liouvilleProperty2 :: Natural -> Bool +liouvilleProperty2 n = m == 0 || l == m + where + l = runFunction liouvilleA n + m = runFunction moebiusA n + +-- | liouville matches baseline from OEIS. +liouvilleOeis :: Assertion +liouvilleOeis = oeisAssertion "A008836" liouvilleA + [ 1, -1, -1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, 1, 1, 1, -1, -1, -1, -1, 1, 1 + , -1, 1, 1, 1, -1, -1, -1, -1, -1, -1, 1, 1, 1, 1, -1, 1, 1, 1, -1, -1, -1, -1 + , -1, 1, -1, -1, 1, -1, 1, -1, -1, 1, 1, 1, 1, 1, -1, 1, -1, 1, -1, 1, 1, -1 + , -1, -1, 1, -1, -1, -1, -1, 1, -1, -1, 1, -1, -1, -1, 1, 1, -1, 1, 1, 1, 1, 1 + , -1, 1, 1, -1, 1, 1, 1, 1, -1, -1, -1, 1, -1 + ] + +-- | carmichaeil divides totient +carmichaelProperty1 :: Natural -> Bool +carmichaelProperty1 n = runFunction totientA n `mod` runFunction carmichaelA n == 0 + +carmichaelSieve100 :: CarmichaelSieve +carmichaelSieve100 = carmichaelSieve 100 + +-- | carmichael matches sieveCarmichael +carmichaelProperty2 :: Natural -> Bool +carmichaelProperty2 n = n < 1 + || fromIntegral (runFunction carmichaelA n) + == sieveCarmichael carmichaelSieve100 (fromIntegral n) + +-- | carmichael matches baseline from OEIS. +carmichaelOeis :: Assertion +carmichaelOeis = oeisAssertion "A002322" carmichaelA + [ 1, 1, 2, 2, 4, 2, 6, 2, 6, 4, 10, 2, 12, 6, 4, 4, 16, 6, 18, 4, 6, 10, 22, 2 + , 20, 12, 18, 6, 28, 4, 30, 8, 10, 16, 12, 6, 36, 18, 12, 4, 40, 6, 42, 10, 12 + , 22, 46, 4, 42, 20, 16, 12, 52, 18, 20, 6, 18, 28, 58, 4, 60, 30, 6, 16, 12 + , 10, 66, 16, 22, 12, 70, 6, 72, 36, 20, 18, 30, 12, 78, 4, 54 + ] + +-- | smallOmega is smaller than bigOmega +omegaProperty1 :: Natural -> Bool +omegaProperty1 n = runFunction smallOmegaA n <= runFunction bigOmegaA n + +-- | smallOmega matches baseline from OEIS. +smallOmegaOeis :: Assertion +smallOmegaOeis = oeisAssertion "A001221" smallOmegaA + [ 0, 1, 1, 1, 1, 2, 1, 1, 1, 2, 1, 2, 1, 2, 2, 1, 1, 2, 1, 2, 2, 2, 1, 2, 1, 2 + , 1, 2, 1, 3, 1, 1, 2, 2, 2, 2, 1, 2, 2, 2, 1, 3, 1, 2, 2, 2, 1, 2, 1, 2, 2, 2 + , 1, 2, 2, 2, 2, 2, 1, 3, 1, 2, 2, 1, 2, 3, 1, 2, 2, 3, 1, 2, 1, 2, 2, 2, 2, 3 + , 1, 2, 1, 2, 1, 3, 2, 2, 2, 2, 1, 3, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 1, 3, 1, 2 + , 3, 2, 1, 2, 1, 3, 2 + ] + +-- | bigOmega matches baseline from OEIS. +bigOmegaOeis :: Assertion +bigOmegaOeis = oeisAssertion "A001222" bigOmegaA + [ 0, 1, 1, 2, 1, 2, 1, 3, 2, 2, 1, 3, 1, 2, 2, 4, 1, 3, 1, 3, 2, 2, 1, 4, 2, 2 + , 3, 3, 1, 3, 1, 5, 2, 2, 2, 4, 1, 2, 2, 4, 1, 3, 1, 3, 3, 2, 1, 5, 2, 3, 2, 3 + , 1, 4, 2, 4, 2, 2, 1, 4, 1, 2, 3, 6, 2, 3, 1, 3, 2, 3, 1, 5, 1, 2, 3, 3, 2, 3 + , 1, 5, 4, 2, 1, 4, 2, 2, 2, 4, 1, 4, 2, 3, 2, 2, 2, 6, 1, 3, 3, 4, 1, 3, 1, 4 + , 3, 2, 1, 5, 1, 3, 2 + ] + +-- | expMangoldt matches baseline from OEIS. +mangoldtOeis :: Assertion +mangoldtOeis = oeisAssertion "A014963" expMangoldtA + [ 1, 2, 3, 2, 5, 1, 7, 2, 3, 1, 11, 1, 13, 1, 1, 2, 17, 1, 19, 1, 1, 1, 23, 1 + , 5, 1, 3, 1, 29, 1, 31, 2, 1, 1, 1, 1, 37, 1, 1, 1, 41, 1, 43, 1, 1, 1, 47, 1 + , 7, 1, 1, 1, 53, 1, 1, 1, 1, 1, 59, 1, 61, 1, 1, 2, 1, 1, 67, 1, 1, 1, 71, 1 + , 73, 1, 1, 1, 1, 1, 79, 1, 3, 1, 83, 1, 1, 1, 1, 1, 89, 1, 1, 1, 1, 1, 1 + ] + +testSuite :: TestTree +testSuite = testGroup "ArithmeticFunctions" + [ testGroup "Divisors" + [ testSmallAndQuick "length . divisors = tau" divisorsProperty1 + , testSmallAndQuick "sum . divisors = sigma_1" divisorsProperty2 + , testSmallAndQuick "matches definition" divisorsProperty3 + , testSmallAndQuick "divisors = divisorsSmall" divisorsProperty4 + ] + , testGroup "Tau" + [ testCase "OEIS" tauOeis + ] + , testGroup "Sigma" + [ testSmallAndQuick "sigma_0 = tau" sigmaProperty1 + , testSmallAndQuick "sigma_1 n > n" sigmaProperty2 + , testCase "OEIS sigma_1" sigma1Oeis + , testCase "OEIS sigma_2" sigma2Oeis + ] + , testGroup "Totient" + [ testSmallAndQuick "totient is even" totientProperty1 + , testSmallAndQuick "totient n < n" totientProperty2 + , testSmallAndQuick "matches sieveTotient" totientProperty3 + , testCase "OEIS" totientOeis + ] + , testGroup "Jordan" + [ testSmallAndQuick "jordan_0 = [== 1]" jordanProperty1 + , testSmallAndQuick "jordan_1 = totient" jordanProperty2 + , testCase "OEIS jordan_2" jordan2Oeis + ] + , testGroup "Moebius" + [ testSmallAndQuick "moebius values" moebiusProperty1 + , testCase "OEIS" moebiusOeis + , testCase "Lazy" moebiusLazy + ] + , testGroup "Liouville" + [ testSmallAndQuick "liouville values" liouvilleProperty1 + , testSmallAndQuick "liouville matches moebius" liouvilleProperty2 + , testCase "OEIS" liouvilleOeis + ] + , testGroup "Carmichael" + [ testSmallAndQuick "carmichael divides totient" carmichaelProperty1 + , testSmallAndQuick "matches sieveCarmichael" carmichaelProperty2 + , testCase "OEIS" carmichaelOeis + ] + , testGroup "Omegas" + [ testSmallAndQuick "smallOmega <= bigOmega" omegaProperty1 + , testCase "OEIS smallOmega" smallOmegaOeis + , testCase "OEIS bigOmega" bigOmegaOeis + ] + , testGroup "Mangoldt" + [ testCase "OEIS" mangoldtOeis + ] + ] diff --git a/test-suite/Math/NumberTheory/ModuliTests.hs b/test-suite/Math/NumberTheory/ModuliTests.hs --- a/test-suite/Math/NumberTheory/ModuliTests.hs +++ b/test-suite/Math/NumberTheory/ModuliTests.hs @@ -8,8 +8,8 @@ -- Tests for Math.NumberTheory.Moduli -- -{-# LANGUAGE CPP #-} -{-# LANGUAGE ViewPatterns #-} +{-# LANGUAGE CPP #-} +{-# LANGUAGE ViewPatterns #-} {-# OPTIONS_GHC -fno-warn-type-defaults #-} @@ -23,51 +23,51 @@ import Data.Bits import Data.List (tails, nub) import Data.Maybe +import Data.Functor.Compose import Math.NumberTheory.Moduli import Math.NumberTheory.TestUtils -toOdd :: Num a => a -> a -toOdd n = n * 2 + 1 - unwrapPP :: (Prime, Power Int) -> (Integer, Int) unwrapPP (Prime p, Power e) = (p, e) -- | Check that 'jacobi' matches 'jacobi''. -jacobiProperty1 :: (Integral a, Bits a) => AnySign a -> NonNegative a -> Bool -jacobiProperty1 (AnySign a) (NonNegative (toOdd -> n)) = n == 1 && j == 1 || n > 1 && j == j' +jacobiProperty1 :: (Integral a, Bits a) => AnySign a -> (Compose Positive Odd) a -> Bool +jacobiProperty1 (AnySign a) (Compose (Positive (Odd n))) = n == 1 && j == 1 || n > 1 && j == j' where j = jacobi a n j' = jacobi' a n -- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 2 -jacobiProperty2 :: (Integral a, Bits a) => AnySign a -> NonNegative a -> Bool -jacobiProperty2 (AnySign a) (NonNegative (toOdd -> n)) = jacobi a n == jacobi (a + n) n +jacobiProperty2 :: (Integral a, Bits a) => AnySign a -> (Compose Positive Odd) a -> Bool +jacobiProperty2 (AnySign a) (Compose (Positive (Odd n))) + = a + n < a -- check overflow + || jacobi a n == jacobi (a + n) n -- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 3 -jacobiProperty3 :: (Integral a, Bits a) => AnySign a -> NonNegative a -> Bool -jacobiProperty3 (AnySign a) (NonNegative (toOdd -> n)) = j == 0 && g /= 1 || abs j == 1 && g == 1 +jacobiProperty3 :: (Integral a, Bits a) => AnySign a -> (Compose Positive Odd) a -> Bool +jacobiProperty3 (AnySign a) (Compose (Positive (Odd n))) = j == 0 && g /= 1 || abs j == 1 && g == 1 where j = jacobi a n g = gcd a n -- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 4 -jacobiProperty4 :: (Integral a, Bits a) => AnySign a -> AnySign a -> NonNegative a -> Bool -jacobiProperty4 (AnySign a) (AnySign b) (NonNegative (toOdd -> n)) = jacobi (a * b) n == jacobi a n * jacobi b n +jacobiProperty4 :: (Integral a, Bits a) => AnySign a -> AnySign a -> (Compose Positive Odd) a -> Bool +jacobiProperty4 (AnySign a) (AnySign b) (Compose (Positive (Odd n))) = jacobi (a * b) n == jacobi a n * jacobi b n -jacobiProperty4_Integer :: AnySign Integer -> AnySign Integer -> NonNegative Integer -> Bool +jacobiProperty4_Integer :: AnySign Integer -> AnySign Integer -> (Compose Positive Odd) Integer -> Bool jacobiProperty4_Integer = jacobiProperty4 -- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 5 -jacobiProperty5 :: (Integral a, Bits a) => AnySign a -> NonNegative a -> NonNegative a -> Bool -jacobiProperty5 (AnySign a) (NonNegative (toOdd -> m)) (NonNegative (toOdd -> n)) = jacobi a (m * n) == jacobi a m * jacobi a n +jacobiProperty5 :: (Integral a, Bits a) => AnySign a -> (Compose Positive Odd) a -> (Compose Positive Odd) a -> Bool +jacobiProperty5 (AnySign a) (Compose (Positive (Odd m))) (Compose (Positive (Odd n))) = jacobi a (m * n) == jacobi a m * jacobi a n -jacobiProperty5_Integer :: AnySign Integer -> NonNegative Integer -> NonNegative Integer -> Bool +jacobiProperty5_Integer :: AnySign Integer -> (Compose Positive Odd) Integer -> (Compose Positive Odd) Integer -> Bool jacobiProperty5_Integer = jacobiProperty5 -- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 6 -jacobiProperty6 :: (Integral a, Bits a) => NonNegative a -> NonNegative a -> Bool -jacobiProperty6 (NonNegative (toOdd -> m)) (NonNegative (toOdd -> n)) = gcd m n /= 1 || jacobi m n * jacobi n m == (if m `mod` 4 == 1 || n `mod` 4 == 1 then 1 else -1) +jacobiProperty6 :: (Integral a, Bits a) => (Compose Positive Odd) a -> (Compose Positive Odd) a -> Bool +jacobiProperty6 (Compose (Positive (Odd m))) (Compose (Positive (Odd n))) = gcd m n /= 1 || jacobi m n * jacobi n m == (if m `mod` 4 == 1 || n `mod` 4 == 1 then 1 else -1) -- | Check that 'invertMod' inverts numbers modulo. invertModProperty :: AnySign Integer -> Positive Integer -> Bool @@ -98,6 +98,10 @@ powerModProperty3 :: (Integral a, Bits a) => AnySign a -> AnySign a -> AnySign Integer -> Positive Integer -> Bool powerModProperty3 (AnySign e1) (AnySign e2) (AnySign b) (Positive m) = (e1 < 0 || e2 < 0) && isNothing (invertMod b m) + || e2 >= 0 && e1 + e2 < e1 -- check overflow + || e1 >= 0 && e1 + e2 < e2 -- check overflow + || e2 <= 0 && e1 + e2 > e1 -- check overflow + || e1 <= 0 && e1 + e2 > e2 -- check overflow || pm1 * pm2 `mod` m == pm12 where pm1 = powerMod b e1 m 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 @@ -19,7 +19,7 @@ import Test.Tasty.QuickCheck as QC hiding (Positive) import Math.NumberTheory.MoebiusInversion.Int -import Math.NumberTheory.Primes.Factorisation +import Math.NumberTheory.ArithmeticFunctions import Math.NumberTheory.TestUtils totientSumProperty :: Positive Int -> Bool 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 @@ -19,7 +19,7 @@ import Test.Tasty.QuickCheck as QC hiding (Positive) import Math.NumberTheory.MoebiusInversion -import Math.NumberTheory.Primes.Factorisation +import Math.NumberTheory.ArithmeticFunctions import Math.NumberTheory.TestUtils totientSumProperty :: Positive Int -> Bool diff --git a/test-suite/Math/NumberTheory/Powers/SquaresTests.hs b/test-suite/Math/NumberTheory/Powers/SquaresTests.hs --- a/test-suite/Math/NumberTheory/Powers/SquaresTests.hs +++ b/test-suite/Math/NumberTheory/Powers/SquaresTests.hs @@ -49,6 +49,10 @@ integerSquareRootProperty_Word :: NonNegative Word -> Bool integerSquareRootProperty_Word = integerSquareRootProperty +-- | Specialized to trigger 'isqrtInteger'. +integerSquareRootProperty_Integer :: NonNegative Integer -> Bool +integerSquareRootProperty_Integer = integerSquareRootProperty + -- | Check that 'integerSquareRoot' returns the largest integer @m@ with @m*m <= n@, where @n@ has form @k@^2-1. integerSquareRootProperty2 :: Integral a => NonNegative a -> Bool integerSquareRootProperty2 (NonNegative k) = n < 0 @@ -65,6 +69,10 @@ integerSquareRootProperty2_Word :: NonNegative Word -> Bool integerSquareRootProperty2_Word = integerSquareRootProperty2 +-- | Specialized to trigger 'isqrtInteger'. +integerSquareRootProperty2_Integer :: NonNegative Integer -> Bool +integerSquareRootProperty2_Integer = integerSquareRootProperty2 + #if WORD_SIZE_IN_BITS == 64 -- | Check that 'integerSquareRoot' of 2^62-1 is 2^31-1, not 2^31. @@ -130,14 +138,16 @@ testSuite :: TestTree testSuite = testGroup "Squares" - [ testGroup "integerSquareRoor" + [ testGroup "integerSquareRoot" [ testIntegralProperty "generic" integerSquareRootProperty , testSmallAndQuick "generic Int" integerSquareRootProperty_Int , testSmallAndQuick "generic Word" integerSquareRootProperty_Word + , testSmallAndQuick "generic Integer" integerSquareRootProperty_Integer - , testIntegralProperty "almost square" integerSquareRootProperty2 - , testSmallAndQuick "almost square Int" integerSquareRootProperty2_Int - , testSmallAndQuick "almost square Word" integerSquareRootProperty2_Word + , testIntegralProperty "almost square" integerSquareRootProperty2 + , testSmallAndQuick "almost square Int" integerSquareRootProperty2_Int + , testSmallAndQuick "almost square Word" integerSquareRootProperty2_Word + , testSmallAndQuick "almost square Integer" integerSquareRootProperty2_Integer #if WORD_SIZE_IN_BITS == 64 , testCase "maxBound / 2 :: Int" integerSquareRootSpecialCase1_Int 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 @@ -31,10 +31,16 @@ sieveFactorSpecialCase1 :: Assertion -sieveFactorSpecialCase1 = assertEqual "sieveFactor" [(29, 1), (73, 1)] $ sieveFactor (factorSieve 2048) (29*73) +sieveFactorSpecialCase1 = do + assertEqual "sieveFactor 2048" [(29, 1), (73, 1)] $ sieveFactor (factorSieve 2048) (29*73) + assertEqual "sieveFactor 15" [(3, 1), (5, 2)] $ sieveFactor (factorSieve 15) (75) +sieveFactorProperty :: Integer -> Bool +sieveFactorProperty n = (n==0) || factorise n == sieveFactor (factorSieve 25) n + testSuite :: TestTree testSuite = testGroup "Primes" [ testSmallAndQuick "primesSum" primesSumProperty - , testCase "sieveFactor" sieveFactorSpecialCase1 + , testCase "sieveFactor special" sieveFactorSpecialCase1 + , testSmallAndQuick "sieveFactor property" sieveFactorProperty ] 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 @@ -12,12 +12,8 @@ {-# LANGUAGE ConstraintKinds #-} {-# LANGUAGE CPP #-} {-# LANGUAGE DataKinds #-} -{-# LANGUAGE DeriveFoldable #-} -{-# LANGUAGE DeriveFunctor #-} -{-# LANGUAGE DeriveTraversable #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-} -{-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE KindSignatures #-} {-# LANGUAGE MultiParamTypeClasses #-} {-# LANGUAGE RankNTypes #-} @@ -36,11 +32,16 @@ {-# OPTIONS_GHC -fno-warn-type-defaults #-} module Math.NumberTheory.TestUtils - ( module Math.NumberTheory.TestUtils + ( module Math.NumberTheory.TestUtils.Wrappers , module Test.SmallCheck.Series , Large(..) + , testIntegralProperty + , testSameIntegralProperty + , testIntegral2Property + , testSmallAndQuick ) where +import Test.SmallCheck.Series (cons2) import Test.Tasty import Test.Tasty.SmallCheck as SC import Test.Tasty.QuickCheck as QC hiding (Positive, NonNegative, generate, getNonNegative) @@ -49,75 +50,41 @@ import Control.Applicative import Data.Bits -#if MIN_VERSION_base(4,8,0) -#else -import Data.Foldable (Foldable) -import Data.Traversable (Traversable) +#if !(MIN_VERSION_base(4,8,0)) import Data.Word #endif import GHC.Exts - -import Math.NumberTheory.Primes - -newtype AnySign a = AnySign { getAnySign :: a } - deriving (Eq, Ord, Read, Show, Num, Enum, Bounded, Integral, Real, Functor, Foldable, Traversable, Arbitrary) - -instance (Monad m, Serial m a) => Serial m (AnySign a) where - series = AnySign <$> series - -instance (Num a, Ord a, Arbitrary a) => Arbitrary (Positive a) where - arbitrary = Positive <$> (arbitrary `suchThat` (> 0)) - shrink (Positive x) = Positive <$> filter (> 0) (shrink x) - -instance (Num a, Ord a, Arbitrary a) => Arbitrary (NonNegative a) where - arbitrary = NonNegative <$> (arbitrary `suchThat` (>= 0)) - shrink (NonNegative x) = NonNegative <$> filter (>= 0) (shrink x) - -instance (Num a, Bounded a) => Bounded (Positive a) where - minBound = Positive 1 - maxBound = Positive (maxBound :: a) - -instance (Num a, Bounded a) => Bounded (NonNegative a) where - minBound = NonNegative 0 - maxBound = NonNegative (maxBound :: a) - -newtype Huge a = Huge { getHuge :: a } - deriving (Eq, Ord, Enum, Bounded, Show, Num, Real, Integral) - -instance (Num a, Arbitrary a) => Arbitrary (Huge a) where - arbitrary = do - Positive l <- arbitrary - ds <- vector l - return $ Huge $ foldl1 (\acc n -> acc * 2^63 + n) ds - -newtype Power a = Power { getPower :: a } - deriving (Eq, Ord, Enum, Bounded, Show, Num, Real, Integral) - -instance (Monad m, Num a, Ord a, Serial m a) => Serial m (Power a) where - series = Power <$> series `suchThatSerial` (> 0) +import Numeric.Natural -instance (Num a, Ord a, Integral a, Arbitrary a) => Arbitrary (Power a) where - arbitrary = Power <$> (getSmall <$> arbitrary) `suchThat` (> 0) - shrink (Power x) = Power <$> filter (> 0) (shrink x) +import Math.NumberTheory.GaussianIntegers (GaussianInteger(..)) -newtype Prime = Prime { getPrime :: Integer } - deriving (Eq, Ord, Show) +import Math.NumberTheory.TestUtils.Compose () +import Math.NumberTheory.TestUtils.Wrappers -instance Arbitrary Prime where - arbitrary = Prime <$> arbitrary `suchThat` (\p -> p > 0 && isPrime p) +instance Monad m => Serial m Word where + series = + generate (\d -> if d >= 0 then pure 0 else empty) <|> nats + where + nats = generate $ \d -> if d > 0 then [1 .. fromInteger (toInteger d)] else empty -instance Monad m => Serial m Prime where - series = Prime <$> series `suchThatSerial` (\p -> p > 0 && isPrime p) +#if !(MIN_VERSION_base(4,8,0)) && !(MIN_VERSION_QuickCheck(2,9,0)) +instance Arbitrary Natural where + arbitrary = fromInteger <$> (arbitrary `suchThat` (>= 0)) + shrink = map fromInteger . filter (>= 0) . shrink . toInteger +#endif -instance Monad m => Serial m Word where +instance Monad m => Serial m Natural where series = generate (\d -> if d >= 0 then pure 0 else empty) <|> nats where nats = generate $ \d -> if d > 0 then [1 .. fromInteger (toInteger d)] else empty -suchThatSerial :: Series m a -> (a -> Bool) -> Series m a -suchThatSerial s p = s >>= \x -> if p x then pure x else empty +instance Arbitrary GaussianInteger where + arbitrary = (:+) <$> arbitrary <*> arbitrary + shrink (x :+ y) = (:+) <$> shrink x <*> shrink y +instance Monad m => Serial m GaussianInteger where + series = cons2 (:+) -- https://www.cs.ox.ac.uk/projects/utgp/school/andres.pdf, p. 21 -- :k Compose = (k1 -> Constraint) -> (k2 -> k1) -> (k2 -> Constraint) @@ -146,8 +113,10 @@ type TestableIntegral wrapper = ( Matrix '[Arbitrary, Show, Serial IO] wrapper '[Int, Word, Integer] + , Matrix '[Arbitrary, Show] wrapper '[Large Int, Large Word, Huge Integer] , Matrix '[Bounded, Integral] wrapper '[Int, Word] , Num (wrapper Integer) + , Functor wrapper ) @@ -176,9 +145,9 @@ , QC.testProperty "quickcheck Int" (f :: wrapper1 Int -> wrapper2 Int -> bool) , QC.testProperty "quickcheck Word" (f :: wrapper1 Word -> wrapper2 Word -> bool) , QC.testProperty "quickcheck Integer" (f :: wrapper1 Integer -> wrapper2 Integer -> bool) - , QC.testProperty "quickcheck Large Int" (\(Large a) (Large b) -> (f :: wrapper1 Int -> wrapper2 Int -> bool) a b) - , QC.testProperty "quickcheck Large Word" (\(Large a) (Large b) -> (f :: wrapper1 Word -> wrapper2 Word -> bool) a b) - , QC.testProperty "quickcheck Huge Integer" (\(Huge a) (Huge b) -> (f :: wrapper1 Integer -> wrapper2 Integer -> bool) a b) + , QC.testProperty "quickcheck Large Int" (\a b -> (f :: wrapper1 Int -> wrapper2 Int -> bool) (getLarge <$> a) (getLarge <$> b)) + , QC.testProperty "quickcheck Large Word" (\a b -> (f :: wrapper1 Word -> wrapper2 Word -> bool) (getLarge <$> a) (getLarge <$> b)) + , QC.testProperty "quickcheck Huge Integer" (\a b -> (f :: wrapper1 Integer -> wrapper2 Integer -> bool) (getHuge <$> a) (getHuge <$> b)) ] testIntegral2Property @@ -205,15 +174,15 @@ , QC.testProperty "quickcheck Integer Word" (f :: wrapper1 Integer -> wrapper2 Word -> bool) , QC.testProperty "quickcheck Integer Integer" (f :: wrapper1 Integer -> wrapper2 Integer -> bool) - , QC.testProperty "quickcheck Large Int Int" ((f :: wrapper1 Int -> wrapper2 Int -> bool) . getLarge) - , QC.testProperty "quickcheck Large Int Word" ((f :: wrapper1 Int -> wrapper2 Word -> bool) . getLarge) - , QC.testProperty "quickcheck Large Int Integer" ((f :: wrapper1 Int -> wrapper2 Integer -> bool) . getLarge) - , QC.testProperty "quickcheck Large Word Int" ((f :: wrapper1 Word -> wrapper2 Int -> bool) . getLarge) - , QC.testProperty "quickcheck Large Word Word" ((f :: wrapper1 Word -> wrapper2 Word -> bool) . getLarge) - , QC.testProperty "quickcheck Large Word Integer" ((f :: wrapper1 Word -> wrapper2 Integer -> bool) . getLarge) - , QC.testProperty "quickcheck Huge Integer Int" ((f :: wrapper1 Integer -> wrapper2 Int -> bool) . getHuge) - , QC.testProperty "quickcheck Huge Integer Word" ((f :: wrapper1 Integer -> wrapper2 Word -> bool) . getHuge) - , QC.testProperty "quickcheck Huge Integer Integer" ((f :: wrapper1 Integer -> wrapper2 Integer -> bool) . getHuge) + , QC.testProperty "quickcheck Large Int Int" ((f :: wrapper1 Int -> wrapper2 Int -> bool) . fmap getLarge) + , QC.testProperty "quickcheck Large Int Word" ((f :: wrapper1 Int -> wrapper2 Word -> bool) . fmap getLarge) + , QC.testProperty "quickcheck Large Int Integer" ((f :: wrapper1 Int -> wrapper2 Integer -> bool) . fmap getLarge) + , QC.testProperty "quickcheck Large Word Int" ((f :: wrapper1 Word -> wrapper2 Int -> bool) . fmap getLarge) + , QC.testProperty "quickcheck Large Word Word" ((f :: wrapper1 Word -> wrapper2 Word -> bool) . fmap getLarge) + , QC.testProperty "quickcheck Large Word Integer" ((f :: wrapper1 Word -> wrapper2 Integer -> bool) . fmap getLarge) + , QC.testProperty "quickcheck Huge Integer Int" ((f :: wrapper1 Integer -> wrapper2 Int -> bool) . fmap getHuge) + , QC.testProperty "quickcheck Huge Integer Word" ((f :: wrapper1 Integer -> wrapper2 Word -> bool) . fmap getHuge) + , QC.testProperty "quickcheck Huge Integer Integer" ((f :: wrapper1 Integer -> wrapper2 Integer -> bool) . fmap getHuge) ] testSmallAndQuick diff --git a/test-suite/Math/NumberTheory/TestUtils/Compose.hs b/test-suite/Math/NumberTheory/TestUtils/Compose.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/NumberTheory/TestUtils/Compose.hs @@ -0,0 +1,47 @@ +-- | +-- Module: Math.NumberTheory.TestUtils.Compose +-- Copyright: (c) 2016 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- Stability: Provisional +-- Portability: Non-portable (GHC extensions) +-- +-- Utils to test Math.NumberTheory +-- + +{-# LANGUAGE CPP #-} +{-# LANGUAGE DeriveGeneric #-} +{-# LANGUAGE FlexibleContexts #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE StandaloneDeriving #-} +{-# LANGUAGE UndecidableInstances #-} + +{-# OPTIONS_GHC -fno-warn-orphans #-} + +module Math.NumberTheory.TestUtils.Compose where + +import Data.Functor.Compose +#if MIN_VERSION_transformers(0,5,0) +#else +import GHC.Generics +#endif + +import Test.Tasty.QuickCheck (Arbitrary) +import Test.SmallCheck.Series (Serial) + +deriving instance Num (f (g a)) => Num (Compose f g a) +deriving instance Enum (f (g a)) => Enum (Compose f g a) +deriving instance Bounded (f (g a)) => Bounded (Compose f g a) + +deriving instance (Ord (Compose f g a), Real (f (g a))) => Real (Compose f g a) +deriving instance (Ord (Compose f g a), Integral (f (g a))) => Integral (Compose f g a) + +deriving instance Arbitrary (f (g a)) => Arbitrary (Compose f g a) + +#if MIN_VERSION_transformers(0,5,0) +#else +deriving instance Generic (Compose f g a) +#endif +instance (Monad m, Serial m (f (g a))) => Serial m (Compose f g a) diff --git a/test-suite/Math/NumberTheory/TestUtils/Wrappers.hs b/test-suite/Math/NumberTheory/TestUtils/Wrappers.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/NumberTheory/TestUtils/Wrappers.hs @@ -0,0 +1,255 @@ +-- | +-- Module: Math.NumberTheory.TestUtils.Wrappers +-- Copyright: (c) 2016 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- Stability: Provisional +-- Portability: Non-portable (GHC extensions) +-- +-- Utils to test Math.NumberTheory +-- + +{-# LANGUAGE CPP #-} +{-# LANGUAGE DeriveFoldable #-} +{-# LANGUAGE DeriveFunctor #-} +{-# LANGUAGE DeriveTraversable #-} +{-# LANGUAGE FlexibleInstances #-} +{-# LANGUAGE GeneralizedNewtypeDeriving #-} +{-# LANGUAGE MultiParamTypeClasses #-} +{-# LANGUAGE ScopedTypeVariables #-} +{-# LANGUAGE StandaloneDeriving #-} + +{-# OPTIONS_GHC -fno-warn-orphans #-} +{-# OPTIONS_GHC -fno-warn-type-defaults #-} + +module Math.NumberTheory.TestUtils.Wrappers where + +import Control.Applicative +import Data.Functor.Classes +#if MIN_VERSION_base(4,8,0) +#else +import Data.Foldable (Foldable) +import Data.Traversable (Traversable) +#endif + +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) + +------------------------------------------------------------------------------- +-- AnySign + +newtype AnySign a = AnySign { getAnySign :: a } + deriving (Eq, Ord, Read, Show, Num, Enum, Bounded, Integral, Real, Functor, Foldable, Traversable, Arbitrary) + +instance (Monad m, Serial m a) => Serial m (AnySign a) where + series = AnySign <$> series + +instance Eq1 AnySign where +#if MIN_VERSION_transformers(0,5,0) + liftEq eq (AnySign a) (AnySign b) = a `eq` b +#else + (AnySign a) `eq1` (AnySign b) = a == b +#endif + +instance Ord1 AnySign where +#if MIN_VERSION_transformers(0,5,0) + liftCompare cmp (AnySign a) (AnySign b) = a `cmp` b +#else + (AnySign a) `compare1` (AnySign b) = a `compare` b +#endif + +instance Show1 AnySign where +#if MIN_VERSION_transformers(0,5,0) + liftShowsPrec shw _ p (AnySign a) = shw p a +#else + showsPrec1 p (AnySign a) = showsPrec p a +#endif + +------------------------------------------------------------------------------- +-- Positive from smallcheck + +deriving instance Functor Positive + +instance (Num a, Ord a, Arbitrary a) => Arbitrary (Positive a) where + arbitrary = Positive <$> (arbitrary `suchThat` (> 0)) + shrink (Positive x) = Positive <$> filter (> 0) (shrink x) + +instance (Num a, Bounded a) => Bounded (Positive a) where + minBound = Positive 1 + maxBound = Positive (maxBound :: a) + +instance Eq1 Positive where +#if MIN_VERSION_transformers(0,5,0) + liftEq eq (Positive a) (Positive b) = a `eq` b +#else + (Positive a) `eq1` (Positive b) = a == b +#endif + +instance Ord1 Positive where +#if MIN_VERSION_transformers(0,5,0) + liftCompare cmp (Positive a) (Positive b) = a `cmp` b +#else + (Positive a) `compare1` (Positive b) = a `compare` b +#endif + +instance Show1 Positive where +#if MIN_VERSION_transformers(0,5,0) + liftShowsPrec shw _ p (Positive a) = shw p a +#else + showsPrec1 p (Positive a) = showsPrec p a +#endif + +------------------------------------------------------------------------------- +-- NonNegative from smallcheck + +deriving instance Functor NonNegative + +instance (Num a, Ord a, Arbitrary a) => Arbitrary (NonNegative a) where + arbitrary = NonNegative <$> (arbitrary `suchThat` (>= 0)) + shrink (NonNegative x) = NonNegative <$> filter (>= 0) (shrink x) + +instance (Num a, Bounded a) => Bounded (NonNegative a) where + minBound = NonNegative 0 + maxBound = NonNegative (maxBound :: a) + +instance Eq1 NonNegative where +#if MIN_VERSION_transformers(0,5,0) + liftEq eq (NonNegative a) (NonNegative b) = a `eq` b +#else + (NonNegative a) `eq1` (NonNegative b) = a == b +#endif + +instance Ord1 NonNegative where +#if MIN_VERSION_transformers(0,5,0) + liftCompare cmp (NonNegative a) (NonNegative b) = a `cmp` b +#else + (NonNegative a) `compare1` (NonNegative b) = a `compare` b +#endif + +instance Show1 NonNegative where +#if MIN_VERSION_transformers(0,5,0) + liftShowsPrec shw _ p (NonNegative a) = shw p a +#else + showsPrec1 p (NonNegative a) = showsPrec p a +#endif + +------------------------------------------------------------------------------- +-- Huge + +newtype Huge a = Huge { getHuge :: a } + deriving (Eq, Ord, Read, Show, Num, Enum, Bounded, Integral, Real, Functor, Foldable, Traversable) + +instance (Num a, Arbitrary a) => Arbitrary (Huge a) where + arbitrary = do + Positive l <- arbitrary + ds <- vector l + return $ Huge $ foldl1 (\acc n -> acc * 2^63 + n) ds + +instance Eq1 Huge where +#if MIN_VERSION_transformers(0,5,0) + liftEq eq (Huge a) (Huge b) = a `eq` b +#else + (Huge a) `eq1` (Huge b) = a == b +#endif + +instance Ord1 Huge where +#if MIN_VERSION_transformers(0,5,0) + liftCompare cmp (Huge a) (Huge b) = a `cmp` b +#else + (Huge a) `compare1` (Huge b) = a `compare` b +#endif + +instance Show1 Huge where +#if MIN_VERSION_transformers(0,5,0) + liftShowsPrec shw _ p (Huge a) = shw p a +#else + showsPrec1 p (Huge a) = showsPrec p a +#endif + +------------------------------------------------------------------------------- +-- Power + +newtype Power a = Power { getPower :: a } + deriving (Eq, Ord, Read, Show, Num, Enum, Bounded, Integral, Real, Functor, Foldable, Traversable) + +instance (Monad m, Num a, Ord a, Serial m a) => Serial m (Power a) where + series = Power <$> series `suchThatSerial` (> 0) + +instance (Num a, Ord a, Integral a, Arbitrary a) => Arbitrary (Power a) where + arbitrary = Power <$> (getSmall <$> arbitrary) `suchThat` (> 0) + shrink (Power x) = Power <$> filter (> 0) (shrink x) + +instance Eq1 Power where +#if MIN_VERSION_transformers(0,5,0) + liftEq eq (Power a) (Power b) = a `eq` b +#else + (Power a) `eq1` (Power b) = a == b +#endif + +instance Ord1 Power where +#if MIN_VERSION_transformers(0,5,0) + liftCompare cmp (Power a) (Power b) = a `cmp` b +#else + (Power a) `compare1` (Power b) = a `compare` b +#endif + +instance Show1 Power where +#if MIN_VERSION_transformers(0,5,0) + liftShowsPrec shw _ p (Power a) = shw p a +#else + showsPrec1 p (Power a) = showsPrec p a +#endif + +------------------------------------------------------------------------------- +-- Odd + +newtype Odd a = Odd { getOdd :: a } + deriving (Eq, Ord, Read, Show, Num, Enum, Bounded, Integral, Real, Functor, Foldable, Traversable) + +instance (Monad m, Serial m a, Integral a) => Serial m (Odd a) where + series = Odd <$> series `suchThatSerial` odd + +instance (Integral a, Arbitrary a) => Arbitrary (Odd a) where + arbitrary = Odd <$> (arbitrary `suchThat` odd) + shrink (Odd x) = Odd <$> filter odd (shrink x) + +instance Eq1 Odd where +#if MIN_VERSION_transformers(0,5,0) + liftEq eq (Odd a) (Odd b) = a `eq` b +#else + (Odd a) `eq1` (Odd b) = a == b +#endif + +instance Ord1 Odd where +#if MIN_VERSION_transformers(0,5,0) + liftCompare cmp (Odd a) (Odd b) = a `cmp` b +#else + (Odd a) `compare1` (Odd b) = a `compare` b +#endif + +instance Show1 Odd where +#if MIN_VERSION_transformers(0,5,0) + liftShowsPrec shw _ p (Odd a) = shw p a +#else + showsPrec1 p (Odd a) = showsPrec p a +#endif + +------------------------------------------------------------------------------- +-- Prime + +newtype Prime = Prime { getPrime :: Integer } + deriving (Eq, Ord, Show) + +instance Arbitrary Prime where + arbitrary = Prime <$> arbitrary `suchThat` (\p -> p > 0 && isPrime p) + +instance Monad m => Serial m Prime where + series = Prime <$> series `suchThatSerial` (\p -> p > 0 && isPrime p) + +------------------------------------------------------------------------------- +-- Utils + +suchThatSerial :: Series m a -> (a -> Bool) -> Series m a +suchThatSerial s p = s >>= \x -> if p x then pure x else empty diff --git a/test-suite/Math/NumberTheory/UniqueFactorisationTests.hs b/test-suite/Math/NumberTheory/UniqueFactorisationTests.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/NumberTheory/UniqueFactorisationTests.hs @@ -0,0 +1,52 @@ +-- | +-- Module: Math.NumberTheory.UniqueFactorisationTests +-- Copyright: (c) 2016 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- Stability: Provisional +-- +-- Tests for Math.NumberTheory.ArithmeticFunctions +-- + +{-# LANGUAGE CPP #-} +{-# LANGUAGE ScopedTypeVariables #-} + +{-# OPTIONS_GHC -fno-warn-type-defaults #-} + +module Math.NumberTheory.UniqueFactorisationTests + ( testSuite + ) where + +import Test.Tasty + +#if MIN_VERSION_base(4,8,0) +#else +import Data.Word +#endif + +import Math.NumberTheory.GaussianIntegers hiding (factorise) +import Math.NumberTheory.UniqueFactorisation +import Math.NumberTheory.TestUtils hiding (Prime) + +import Numeric.Natural + +testRules :: forall a. (UniqueFactorisation a, Num a, Eq a) => a -> Bool +testRules n + = n == 0 + || all (\(p, _) -> unP p == abs (unP p)) fs + && abs n == abs (product (map (\(p, k) -> unP p ^ k) fs)) + where + fs = factorise n + + unP :: Prime a -> a + unP = unPrime + +testSuite :: TestTree +testSuite = testGroup "UniqueFactorisation" + [ testSmallAndQuick "Int" (testRules :: Int -> Bool) + , testSmallAndQuick "Word" (testRules :: Word -> Bool) + , testSmallAndQuick "Integer" (testRules :: Integer -> Bool) + , testSmallAndQuick "Natural" (testRules :: Natural -> Bool) + + , testSmallAndQuick "GaussianInteger" (testRules :: GaussianInteger -> Bool) + ] diff --git a/test-suite/Test.hs b/test-suite/Test.hs --- a/test-suite/Test.hs +++ b/test-suite/Test.hs @@ -25,6 +25,9 @@ import qualified Math.NumberTheory.GaussianIntegersTests as Gaussian +import qualified Math.NumberTheory.ArithmeticFunctionsTests as ArithmeticFunctions +import qualified Math.NumberTheory.UniqueFactorisationTests as UniqueFactorisation + main :: IO () main = defaultMain tests @@ -62,5 +65,11 @@ ] , testGroup "Gaussian" [ Gaussian.testSuite + ] + , testGroup "ArithmeticFunctions" + [ ArithmeticFunctions.testSuite + ] + , testGroup "UniqueFactorisation" + [ UniqueFactorisation.testSuite ] ]