diff --git a/Changes b/Changes --- a/Changes +++ b/Changes @@ -1,3 +1,40 @@ +0.5.0.0: + This release supports GHC 7.8, 7.10 and 8.0. GHC 7.6 is no longer supported. + + Breaking changes: + + Remove deprecated interface to arithmetic functions (divisors, tau, + sigma, totient, jordan, moebius, liouville, smallOmega, bigOmega, + carmichael, expMangoldt). New interface is exposed via + Math.NumberTheory.ArithmeticFunctions (#30). + + Deprecate integerPower and integerWordPower from + Math.NumberTheory.Powers.Integer. Use (^) instead (#51). + + Math.NumberTheory.Logarithms has been moved to the separate package + integer-logarithms (#51). + + Rename Math.NumberTheory.Lucas to Math.NumberTheory.Recurrencies.Linear. + + New functions: + + Add basic combinatorial sequences: binomial coefficients, Stirling + numbers of both kinds, Eulerian numbers of both kinds, Bernoulli + numbers (#39). E. g., + + > take 10 $ Math.NumberTheory.Recurrencies.Bilinear.bernoulli + [1 % 1,(-1) % 2,1 % 6,0 % 1,(-1) % 30,0 % 1,1 % 42,0 % 1,(-1) % 30,0 % 1] + + Add the Riemann zeta function on non-negative integer arguments (#44). + E. g., + + > take 5 $ Math.NumberTheory.Zeta.zetas 1e-15 + [-0.5,Infinity,1.6449340668482262,1.2020569031595945,1.0823232337111381] + + Improvements: + + Speed up isPrime twice; rework millerRabinV and isStrongFermatPP (#22, #25). + 0.4.3.0: This release supports GHC 7.6, 7.8, 7.10 and 8.0. diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ -Copyright (c) 2011 Daniel Fischer +Copyright (c) 2011 Daniel Fischer, 2016-2017 Andrew Lelechenko, Carter Schonwald Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, diff --git a/Math/NumberTheory/ArithmeticFunctions/Standard.hs b/Math/NumberTheory/ArithmeticFunctions/Standard.hs --- a/Math/NumberTheory/ArithmeticFunctions/Standard.hs +++ b/Math/NumberTheory/ArithmeticFunctions/Standard.hs @@ -36,6 +36,7 @@ , expMangoldt, expMangoldtA ) where +import Data.Coerce import Data.IntSet (IntSet) import qualified Data.IntSet as IS import Data.Set (Set) @@ -51,15 +52,6 @@ #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 diff --git a/Math/NumberTheory/Logarithms.hs b/Math/NumberTheory/Logarithms.hs deleted file mode 100644 --- a/Math/NumberTheory/Logarithms.hs +++ /dev/null @@ -1,195 +0,0 @@ --- | --- Module: Math.NumberTheory.Logarithms --- Copyright: (c) 2011 Daniel Fischer --- Licence: MIT --- Maintainer: Daniel Fischer --- Stability: Provisional --- Portability: Non-portable (GHC extensions) --- --- Integer Logarithms. For efficiency, the internal representation of 'Integer's --- from integer-gmp is used. --- -{-# LANGUAGE CPP, MagicHash #-} -module Math.NumberTheory.Logarithms - ( -- * Integer logarithms with input checks - integerLogBase - , integerLog2 - , integerLog10 - - , intLog2 - , wordLog2 - - -- * Integer logarithms without input checks - , integerLogBase' - , integerLog2' - , integerLog10' - - , intLog2' - , wordLog2' - ) where - -import GHC.Base - -import Data.Bits -import Data.Array.Unboxed - -import GHC.Integer.Logarithms - -import Math.NumberTheory.Powers.Integer -import Math.NumberTheory.Unsafe -#if __GLASGOW_HASKELL__ < 707 -import Math.NumberTheory.Utils (isTrue#) -#endif - --- | Calculate the integer logarithm for an arbitrary base. --- The base must be greater than 1, the second argument, the number --- whose logarithm is sought, must be positive, otherwise an error is thrown. --- If @base == 2@, the specialised version is called, which is more --- efficient than the general algorithm. --- --- Satisfies: --- --- > base ^ integerLogBase base m <= m < base ^ (integerLogBase base m + 1) --- --- for @base > 1@ and @m > 0@. -integerLogBase :: Integer -> Integer -> Int -integerLogBase b n - | n < 1 = error "Math.NumberTheory.Logarithms.integerLogBase: argument must be positive." - | n < b = 0 - | b == 2 = integerLog2' n - | b < 2 = error "Math.NumberTheory.Logarithms.integerLogBase: base must be greater than one." - | otherwise = integerLogBase' b n - --- | Calculate the integer logarithm of an 'Integer' to base 2. --- The argument must be positive, otherwise an error is thrown. -integerLog2 :: Integer -> Int -integerLog2 n - | n < 1 = error "Math.NumberTheory.Logarithms.integerLog2: argument must be positive" - | otherwise = I# (integerLog2# n) - --- | Calculate the integer logarithm of an 'Int' to base 2. --- The argument must be positive, otherwise an error is thrown. -intLog2 :: Int -> Int -intLog2 (I# i#) - | isTrue# (i# <# 1#) = error "Math.NumberTheory.Logarithms.intLog2: argument must be positive" - | otherwise = I# (wordLog2# (int2Word# i#)) - --- | Calculate the integer logarithm of a 'Word' to base 2. --- The argument must be positive, otherwise an error is thrown. -wordLog2 :: Word -> Int -wordLog2 (W# w#) - | isTrue# (w# `eqWord#` 0##) = error "Math.NumberTheory.Logarithms.wordLog2: argument must not be 0." - | otherwise = I# (wordLog2# w#) - --- | Same as 'integerLog2', but without checks, saves a little time when --- called often for known good input. -integerLog2' :: Integer -> Int -integerLog2' n = I# (integerLog2# n) - --- | Same as 'intLog2', but without checks, saves a little time when --- called often for known good input. -intLog2' :: Int -> Int -intLog2' (I# i#) = I# (wordLog2# (int2Word# i#)) - --- | Same as 'wordLog2', but without checks, saves a little time when --- called often for known good input. -wordLog2' :: Word -> Int -wordLog2' (W# w#) = I# (wordLog2# w#) - --- | Calculate the integer logarithm of an 'Integer' to base 10. --- The argument must be positive, otherwise an error is thrown. -integerLog10 :: Integer -> Int -integerLog10 n - | n < 1 = error "Math.NumberTheory.Logarithms.integerLog10: argument must be positive" - | otherwise = integerLog10' n - --- | Same as 'integerLog10', but without a check for a positive --- argument. Saves a little time when called often for known good --- input. -integerLog10' :: Integer -> Int -integerLog10' n - | n < 10 = 0 - | n < 100 = 1 - | otherwise = ex + integerLog10' (n `quot` integerPower 10 ex) - where - ln = I# (integerLog2# n) - -- u/v is a good approximation of log 2/log 10 - u = 1936274 - v = 6432163 - -- so ex is a good approximation to integerLogBase 10 n - ex = fromInteger ((u * fromIntegral ln) `quot` v) - --- | Same as 'integerLogBase', but without checks, saves a little time when --- called often for known good input. -integerLogBase' :: Integer -> Integer -> Int -integerLogBase' b n - | n < b = 0 - | ln-lb < lb = 1 -- overflow safe version of ln < 2*lb, implies n < b*b - | b < 33 = let bi = fromInteger b - ix = 2*bi-4 - -- u/v is a good approximation of log 2/log b - u = logArr `unsafeAt` ix - v = logArr `unsafeAt` (ix+1) - -- hence ex is a rather good approximation of integerLogBase b n - -- most of the time, it will already be exact - ex = fromInteger ((fromIntegral u * fromIntegral ln) `quot` fromIntegral v) - in case u of - 1 -> ln `quot` v -- a power of 2, easy - _ -> ex + integerLogBase' b (n `quot` integerPower b ex) - | otherwise = let -- shift b so that 16 <= bi < 32 - bi = fromInteger (b `shiftR` (lb-4)) - -- we choose an approximation of log 2 / log (bi+1) to - -- be sure we underestimate - ix = 2*bi-2 - -- u/w is a reasonably good approximation to log 2/log b - -- it is too small, but not by much, so the recursive call - -- should most of the time be caught by one of the first - -- two guards unless n is huge, but then it'd still be - -- a call with a much smaller second argument. - u = fromIntegral $ logArr `unsafeAt` ix - v = fromIntegral $ logArr `unsafeAt` (ix+1) - w = v + u*fromIntegral (lb-4) - ex = fromInteger ((u * fromIntegral ln) `quot` w) - in ex + integerLogBase' b (n `quot` integerPower b ex) - where - lb = integerLog2' b - ln = integerLog2' n - --- Lookup table for logarithms of 2 <= k <= 32 --- In each row "x , y", x/y is a good rational approximation of log 2 / log k. --- For the powers of 2, it is exact, otherwise x/y < log 2/log k, since we don't --- want to overestimate integerLogBase b n = floor $ (log 2/log b)*logBase 2 n. -logArr :: UArray Int Int -logArr = listArray (0, 61) - [ 1 , 1, - 190537 , 301994, - 1 , 2, - 1936274 , 4495889, - 190537 , 492531, - 91313 , 256348, - 1 , 3, - 190537 , 603988, - 1936274 , 6432163, - 1686227 , 5833387, - 190537 , 683068, - 5458 , 20197, - 91313 , 347661, - 416263 , 1626294, - 1 , 4, - 32631 , 133378, - 190537 , 794525, - 163451 , 694328, - 1936274 , 8368437, - 1454590 , 6389021, - 1686227 , 7519614, - 785355 , 3552602, - 190537 , 873605, - 968137 , 4495889, - 5458 , 25655, - 190537 , 905982, - 91313 , 438974, - 390321 , 1896172, - 416263 , 2042557, - 709397 , 3514492, - 1 , 5 - ] diff --git a/Math/NumberTheory/Lucas.hs b/Math/NumberTheory/Lucas.hs deleted file mode 100644 --- a/Math/NumberTheory/Lucas.hs +++ /dev/null @@ -1,99 +0,0 @@ --- | --- Module: Math.NumberTheory.Lucas --- Copyright: (c) 2011 Daniel Fischer --- Licence: MIT --- Maintainer: Daniel Fischer --- Stability: Provisional --- Portability: Non-portable (GHC extensions) --- --- Efficient calculation of Lucas sequences. -{-# LANGUAGE CPP #-} -module Math.NumberTheory.Lucas - ( fibonacci - , fibonacciPair - , lucas - , lucasPair - , generalLucas - ) where - -#include "MachDeps.h" - -import Data.Bits - --- | @'fibonacci' k@ calculates the @k@-th Fibonacci number in --- /O/(@log (abs k)@) steps. The index may be negative. This --- is efficient for calculating single Fibonacci numbers (with --- large index), but for computing many Fibonacci numbers in --- close proximity, it is better to use the simple addition --- formula starting from an appropriate pair of successive --- Fibonacci numbers. -fibonacci :: Int -> Integer -fibonacci = fst . fibonacciPair - --- | @'fibonacciPair' k@ returns the pair @(F(k), F(k+1))@ of the @k@-th --- Fibonacci number and its successor, thus it can be used to calculate --- the Fibonacci numbers from some index on without needing to compute --- the previous. The pair is efficiently calculated --- in /O/(@log (abs k)@) steps. The index may be negative. -fibonacciPair :: Int -> (Integer, Integer) -fibonacciPair n - | n < 0 = let (f,g) = fibonacciPair (-(n+1)) in if testBit n 0 then (g, -f) else (-g, f) - | n == 0 = (0, 1) - | otherwise = look (WORD_SIZE_IN_BITS - 2) - where - look k - | testBit n k = go (k-1) 0 1 - | otherwise = look (k-1) - go k g f - | k < 0 = (f, f+g) - | testBit n k = go (k-1) (f*(f+shiftL g 1)) ((f+g)*shiftL f 1 + g*g) - | otherwise = go (k-1) (f*f+g*g) (f*(f+shiftL g 1)) - --- | @'lucas' k@ computes the @k@-th Lucas number. Very similar --- to @'fibonacci'@. -lucas :: Int -> Integer -lucas = fst . lucasPair - --- | @'lucasPair' k@ computes the pair @(L(k), L(k+1))@ of the @k@-th --- Lucas number and its successor. Very similar to @'fibonacciPair'@. -lucasPair :: Int -> (Integer, Integer) -lucasPair n - | n < 0 = let (f,g) = lucasPair (-(n+1)) in if testBit n 0 then (-g, f) else (g, -f) - | n == 0 = (2, 1) - | otherwise = look (WORD_SIZE_IN_BITS - 2) - where - look k - | testBit n k = go (k-1) 0 1 - | otherwise = look (k-1) - go k g f - | k < 0 = (shiftL g 1 + f,g+3*f) - | otherwise = go (k-1) g' f' - where - (f',g') - | testBit n k = (shiftL (f*(f+g)) 1 + g*g,f*(shiftL g 1 + f)) - | otherwise = (f*(shiftL g 1 + f),f*f+g*g) - - --- | @'generalLucas' p q k@ calculates the quadruple @(U(k), U(k+1), V(k), V(k+1))@ --- where @U(i)@ is the Lucas sequence of the first kind and @V(i)@ the Lucas --- sequence of the second kind for the parameters @p@ and @q@, where @p^2-4q /= 0@. --- Both sequences satisfy the recurrence relation @A(j+2) = p*A(j+1) - q*A(j)@, --- the starting values are @U(0) = 0, U(1) = 1@ and @V(0) = 2, V(1) = p@. --- The Fibonacci numbers form the Lucas sequence of the first kind for the --- parameters @p = 1, q = -1@ and the Lucas numbers form the Lucas sequence of --- the second kind for these parameters. --- Here, the index must be non-negative, since the terms of the sequence for --- negative indices are in general not integers. -generalLucas :: Integer -> Integer -> Int -> (Integer, Integer, Integer, Integer) -generalLucas p q k - | k < 0 = error "generalLucas: negative index" - | k == 0 = (0,1,2,p) - | otherwise = look (WORD_SIZE_IN_BITS - 2) - where - look i - | testBit k i = go (i-1) 1 p p q - | otherwise = look (i-1) - go i un un1 vn qn - | i < 0 = (un, un1, vn, p*un1 - shiftL (q*un) 1) - | testBit k i = go (i-1) (un1*vn-qn) ((p*un1-q*un)*vn - p*qn) ((p*un1 - (2*q)*un)*vn - p*qn) (qn*qn*q) - | otherwise = go (i-1) (un*vn) (un1*vn-qn) (vn*vn - 2*qn) (qn*qn) 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 @@ -31,9 +31,6 @@ import GHC.Integer.Logarithms (integerLog2#) import Math.NumberTheory.Unsafe -#if __GLASGOW_HASKELL__ < 707 -import Math.NumberTheory.Utils (isTrue#) -#endif -- | Calculate the integer cube root of an integer @n@, -- that is the largest integer @r@ such that @r^3 <= n@. 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 @@ -31,9 +31,6 @@ import Data.Bits import Math.NumberTheory.Unsafe -#if __GLASGOW_HASKELL__ < 707 -import Math.NumberTheory.Utils (isTrue#) -#endif -- | Calculate the integer fourth root of a nonnegative number, -- that is, the largest integer @r@ with @r^4 <= n@. 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 @@ -34,9 +34,6 @@ import Math.NumberTheory.Logarithms (integerLogBase') import Math.NumberTheory.Utils (shiftToOddCount , splitOff -#if __GLASGOW_HASKELL__ < 707 - , isTrue# -#endif ) import qualified Math.NumberTheory.Powers.Squares as P2 import qualified Math.NumberTheory.Powers.Cubes as P3 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 @@ -9,17 +9,16 @@ -- Potentially faster power function for 'Integer' base and 'Int' -- or 'Word' exponent. -- -{-# LANGUAGE MagicHash, BangPatterns, CPP #-} +{-# LANGUAGE CPP #-} module Math.NumberTheory.Powers.Integer + {-# DEPRECATED "It is no faster than (^)" #-} ( integerPower , integerWordPower ) where -import GHC.Base -import GHC.Integer.Logarithms (wordLog2#) - -#if __GLASGOW_HASKELL__ < 707 -import Math.NumberTheory.Utils (isTrue#) +#if MIN_VERSION_base(4,8,0) +#else +import Data.Word #endif -- | Power of an 'Integer' by the left-to-right repeated squaring algorithm. @@ -36,23 +35,10 @@ -- /Warning:/ No check for the negativity of the exponent is performed, -- a negative exponent is interpreted as a large positive exponent. integerPower :: Integer -> Int -> Integer -integerPower b (I# e#) = power b (int2Word# e#) +integerPower = (^) +{-# DEPRECATED integerPower "Use (^) instead" #-} -- | Same as 'integerPower', but for exponents of type 'Word'. integerWordPower :: Integer -> Word -> Integer -integerWordPower b (W# w#) = power b w# - -power :: Integer -> Word# -> Integer -power b w# - | isTrue# (w# `eqWord#` 0##) = 1 - | isTrue# (w# `eqWord#` 1##) = b - | otherwise = go (wordLog2# w# -# 1#) b (b*b) - where - go 0# l h = if isTrue# ((w# `and#` 1##) `eqWord#` 0##) then l*l else (l*h) - go i# l h - | w# `hasBit#` i# = go (i# -# 1#) (l*h) (h*h) - | otherwise = go (i# -# 1#) (l*l) (l*h) - --- | A raw version of testBit for 'Word#'. -hasBit# :: Word# -> Int# -> Bool -hasBit# w# i# = isTrue# (((w# `uncheckedShiftRL#` i#) `and#` 1##) `neWord#` 0##) +integerWordPower = (^) +{-# DEPRECATED integerWordPower "Use (^) instead" #-} diff --git a/Math/NumberTheory/Powers/Squares/Internal.hs b/Math/NumberTheory/Powers/Squares/Internal.hs --- a/Math/NumberTheory/Powers/Squares/Internal.hs +++ b/Math/NumberTheory/Powers/Squares/Internal.hs @@ -29,9 +29,6 @@ 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 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 @@ -14,8 +14,6 @@ -- 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 @@ -38,47 +36,19 @@ -- *** Single curve worker , montgomeryFactorisation -- * Totients - , totient - , φ , TotientSieve , totientSieve , sieveTotient - , totientFromCanonical -- * Carmichael function - , carmichael - , λ , CarmichaelSieve , carmichaelSieve , sieveCarmichael - , carmichaelFromCanonical - -- * Moebius function - , moebius - , μ - , moebiusFromCanonical - -- * Divisors - , divisors - , tau - , τ - , divisorCount - , divisorSum - , sigma - , σ - , divisorPowerSum - , divisorsFromCanonical - , tauFromCanonical - , divisorSumFromCanonical - , sigmaFromCanonical ) where -import Data.Set (Set, singleton) - -import Math.NumberTheory.Primes.Factorisation.Utils 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 -- -- Factorisation of 'Integer's by the elliptic curve algorithm after Montgomery. @@ -95,95 +65,3 @@ -- -- Given enough time, the algorithm should be able to factor numbers of 100-120 digits, but it -- is best suited for numbers of up to 50-60 digits. - --- | 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)@. -totient :: Integer -> Integer -totient n - | n < 1 = error "Totient only defined for positive numbers" - | n == 1 = 1 - | otherwise = totientFromCanonical (factorise' n) - --- | Alias of 'totient' for people who prefer Greek letters. -φ :: Integer -> Integer -φ = totient - --- | Calculates the Carmichael function for a positive integer, that is, --- the (smallest) exponent of the group of units in @ℤ/(n)@. -carmichael :: Integer -> Integer -carmichael n - | n < 1 = error "Carmichael function only defined for positive numbers" - | n == 1 = 1 - | otherwise = carmichaelFromCanonical (factorise' n) - --- | Alias of 'carmichael' for people who prefer Greek letters. -λ :: Integer -> Integer -λ = carmichael - --- | Calculates the Moebius function for a positive integer. -moebius :: Integer -> Integer -moebius n - | n < 1 = error "Carmichael function only defined for positive numbers" - | n == 1 = 1 - | otherwise = moebiusFromCanonical (factorise' n) - --- | Alias of 'moebius' for people who prefer Greek letters. -μ :: Integer -> Integer -μ = moebius - --- | @'divisors' n@ is the set of all (positive) divisors of @n@. --- @'divisors' 0@ is an error because we can't create the set of all 'Integer's. -divisors :: Integer -> Set Integer -divisors n - | n < 0 = divisors (-n) - | n == 0 = error "Can't create set of divisors of 0" - | n == 1 = singleton 1 - | otherwise = divisorsFromCanonical (factorise' n) - --- | @'tau' n@ is the number of (positive) divisors of @n@. --- @'tau' 0@ is an error because @0@ has infinitely many divisors. -tau :: Integer -> Integer -tau n - | n < 0 = tau (-n) - | n == 0 = error "0 has infinitely many divisors" - | n == 1 = 1 - | otherwise = tauFromCanonical (factorise' n) - --- | Alias for 'tau'. -divisorCount :: Integer -> Integer -divisorCount = tau - --- | The sum of all (positive) divisors of a positive number @n@, --- calculated from its prime factorisation. -divisorSum :: Integer -> Integer -divisorSum n - | n < 1 = error "divisor sum only defined for positive numbers" - | n == 1 = 1 - | otherwise = divisorSumFromCanonical (factorise' n) - --- | Alias for 'sigma'. -divisorPowerSum :: Int -> Integer -> Integer -divisorPowerSum = sigma - --- | @'sigma' k n@ is the sum of the @k@-th powers of the --- (positive) divisors of @n@. @k@ must be non-negative and @n@ positive. --- For @k == 0@, it is the divisor count (@d^0 = 1@). -sigma :: Int -> Integer -> Integer -sigma 0 n = tau n -sigma 1 n = divisorSum n -sigma k n - | k < 0 = error "sigma: exponent must be non-negative" - | n < 1 = error "sigma: n must be positive" - | n == 1 = 1 - | otherwise = sigmaFromCanonical k (factorise' n) - --- | Alias for 'sigma' for people preferring Greek letters. -σ :: Int -> Integer -> Integer -σ 0 = divisorCount -σ 1 = divisorSum -σ k = divisorPowerSum k - --- | Alias for 'tau' for people preferring Greek letters. -τ :: Integer -> Integer -τ = tau diff --git a/Math/NumberTheory/Primes/Factorisation/Utils.hs b/Math/NumberTheory/Primes/Factorisation/Utils.hs deleted file mode 100644 --- a/Math/NumberTheory/Primes/Factorisation/Utils.hs +++ /dev/null @@ -1,89 +0,0 @@ --- | --- Module: Math.NumberTheory.Primes.Factorisation.Utils --- Copyright: (c) 2011 Daniel Fischer --- Licence: MIT --- Maintainer: Daniel Fischer --- Stability: Provisional --- Portability: Non-portable (GHC extensions) --- --- Some utilities related to factorisation, defined here to avoid import cycles. -{-# LANGUAGE BangPatterns #-} -{-# OPTIONS_HADDOCK hide #-} -module Math.NumberTheory.Primes.Factorisation.Utils - ( ppTotient - , totientFromCanonical - , carmichaelFromCanonical - , moebiusFromCanonical - , divisorsFromCanonical - , tauFromCanonical - , divisorSumFromCanonical - , sigmaFromCanonical - ) where - -import Data.Set (Set) -import qualified Data.Set as Set -import Data.Bits -import Data.List - -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 -ppTotient (p,k) = (p-1)*(integerPower p (k-1)) -- slightly faster than (^) usually - --- | Calculate the totient from the canonical factorisation. -totientFromCanonical :: [(Integer,Int)] -> Integer -totientFromCanonical = product . map ppTotient - --- | Calculate the Carmichael function from the factorisation. --- Requires that the list of prime factors is strictly ascending. -carmichaelFromCanonical :: [(Integer,Int)] -> Integer -carmichaelFromCanonical = go2 - where - go2 ((2,k):ps) = let acc = case k of - 1 -> 1 - 2 -> 2 - _ -> 1 `shiftL` (k-2) - in go acc ps - go2 ps = go 1 ps - go !acc ((p,1):pps) = go (lcm acc (p-1)) pps - go acc ((p,k):pps) = go ((lcm acc (p-1))*integerPower p (k-1)) pps - go acc [] = acc - --- | Calculate the Moebius function from the canonical factorisation. -moebiusFromCanonical :: [(a, Int)] -> Integer -moebiusFromCanonical = go 1 - where - go acc [] = acc - go acc ((_, 1) : xs) = go (negate acc) xs - go acc ((_, 0) : xs) = go acc xs -- Should not really happen - go _ _ = 0 -- Short circuit for powers > 1 - --- | The set of divisors, efficiently calculated from the canonical factorisation. -divisorsFromCanonical :: [(Integer,Int)] -> Set Integer -divisorsFromCanonical = foldl' step (Set.singleton 1) - where - step st (p,k) = Set.unions (st:[Set.mapMonotonic (*pp) st | pp <- take k (iterate (*p) p) ]) - --- | The number of divisors, efficiently calculated from the canonical factorisation. -tauFromCanonical :: [(a,Int)] -> Integer -tauFromCanonical pps = product [fromIntegral k + 1 | (_,k) <- pps] - --- | The sum of all divisors, efficiently calculated from the canonical factorisation. -divisorSumFromCanonical :: [(Integer,Int)] -> Integer -divisorSumFromCanonical = product . map ppDivSum - -ppDivSum :: (Integer,Int) -> Integer -ppDivSum (p,1) = p+1 -ppDivSum (p,k) = (p^(k+1)-1) `quot` (p-1) - --- | The sum of the powers (with fixed exponent) of all divisors, --- efficiently calculated from the canonical factorisation. -sigmaFromCanonical :: Int -> [(Integer,Int)] -> Integer -sigmaFromCanonical k = product . map (ppDivPowerSum k) - -ppDivPowerSum :: Int -> (Integer,Int) -> Integer -ppDivPowerSum k (p,m) = (p^(k*(m+1)) - 1) `quot` (p^k - 1) 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 @@ -11,7 +11,6 @@ {-# LANGUAGE MonoLocalBinds #-} {-# LANGUAGE ScopedTypeVariables #-} {-# OPTIONS_GHC -fspec-constr-count=8 #-} -{-# OPTIONS_GHC -fno-warn-warnings-deprecations #-} {-# OPTIONS_HADDOCK hide #-} module Math.NumberTheory.Primes.Sieve.Misc ( -- * Types @@ -44,7 +43,6 @@ import Math.NumberTheory.Powers.Squares (integerSquareRoot') import Math.NumberTheory.Primes.Sieve.Indexing import Math.NumberTheory.Primes.Factorisation.Montgomery -import Math.NumberTheory.Primes.Factorisation.Utils import Math.NumberTheory.Unsafe import Math.NumberTheory.Utils @@ -232,6 +230,15 @@ pix = unsafeAt sve ix curve tt n = tt * totientFromCanonical (stdGenFactorisation (Just (bound*(bound+2))) (mkStdGen $ fromIntegral n `xor` 0xdecaf00d) Nothing n) +-- | Calculate the totient from the canonical factorisation. +totientFromCanonical :: [(Integer,Int)] -> Integer +totientFromCanonical = product . map ppTotient + +-- | Totient of a prime power. +ppTotient :: (Integer, Int) -> Integer +ppTotient (p, 1) = p - 1 +ppTotient (p, k) = (p - 1) * p ^ (k - 1) + -- | @'carmichaelSieve' n@ creates a store of values of the Carmichael function -- for numbers not exceeding @n@. -- Like a 'TotientSieve', a 'CarmichaelSieve' only stores values for numbers coprime to @30@ @@ -299,6 +306,20 @@ pix = unsafeAt sve ix curve tt n = tt `lcm` carmichaelFromCanonical (stdGenFactorisation (Just (bound*(bound+2))) (mkStdGen $ fromIntegral n `xor` 0xdecaf00d) Nothing n) +-- | Calculate the Carmichael function from the factorisation. +-- Requires that the list of prime factors is strictly ascending. +carmichaelFromCanonical :: [(Integer, Int)] -> Integer +carmichaelFromCanonical = go2 + where + go2 ((2, k) : ps) = let acc = case k of + 1 -> 1 + 2 -> 2 + _ -> 1 `shiftL` (k-2) + in go acc ps + go2 ps = go 1 ps + go !acc ((p, 1) : pps) = go (lcm acc (p - 1)) pps + go acc ((p, k) : pps) = go ((lcm acc (p - 1)) * p ^ (k - 1)) pps + go acc [] = acc -- NOTE: This is a legacy implementation of FactorSieve which uses the -- same (2,3,5) wheel optimization as the other sieves. 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 @@ -1,6 +1,6 @@ -- | -- Module: Math.NumberTheory.Primes.Testing.Probabilistic --- Copyright: (c) 2011 Daniel Fischer +-- Copyright: (c) 2011 Daniel Fischer, 2017 Andrew Lelechenko -- Licence: MIT -- Maintainer: Daniel Fischer -- Stability: Provisional @@ -20,59 +20,47 @@ #include "MachDeps.h" -import Math.NumberTheory.Moduli -import Math.NumberTheory.Utils -import Math.NumberTheory.Powers.Squares - import Data.Bits - import GHC.Base - import GHC.Integer.GMP.Internals --- | @'isPrime' n@ tests whether @n@ is a prime (negative or positive). --- First, trial division by the primes less than @1200@ is performed. --- If that hasn't determined primality or compositeness, a Baillie PSW --- test is performed. +import Math.NumberTheory.Moduli +import Math.NumberTheory.Utils +import Math.NumberTheory.Powers.Squares + +-- | @isPrime n@ tests whether @n@ is a prime (negative or positive). +-- It is a combination of trial division and Baillie-PSW test. -- --- Since the Baillie PSW test may not be perfect, it is possible that --- some large composites are wrongly deemed prime, however, no composites --- passing the test are known and none exist below @2^64@. +-- If @isPrime n@ returns @False@ then @n@ is definitely composite. +-- There is a theoretical possibility that @isPrime n@ is @True@, +-- but in fact @n@ is not prime. However, no such numbers are known +-- and none exist below @2^64@. If you have found one, please report it, +-- because it is a major discovery. isPrime :: Integer -> Bool isPrime n | n < 0 = isPrime (-n) | n < 2 = False | n < 4 = True - | otherwise = go smallPrimes - where - go (p:ps) - | p*p > n = True - | otherwise = case n `rem` p of - 0 -> False - _ -> go ps - go [] = bailliePSW n + | otherwise = millerRabinV 0 n -- trial division test + && bailliePSW n --- | A Miller-Rabin like probabilistic primality test with preceding --- trial division. While the classic Miller-Rabin test uses --- randomly chosen bases, @'millerRabinV' k n@ uses the @k@ --- smallest primes as bases if trial division has not reached --- a conclusive result. (Only the primes up to @1200@ are --- available in this module, so the maximal effective @k@ is @196@.) +-- | Miller-Rabin probabilistic primality test. It consists of the trial +-- division test and several rounds of the strong Fermat test with different +-- bases. The choice of trial divisors and bases are +-- implementation details and may change in future silently. +-- +-- First argument stands for the number of rounds of strong Fermat test. +-- If it is 0, only trial division test is performed. +-- +-- If @millerRabinV k n@ returns @False@ then @n@ is definitely composite. +-- Otherwise @n@ may appear composite with probability @1/4^k@. millerRabinV :: Int -> Integer -> Bool -millerRabinV k n - | n < 0 = millerRabinV k (-n) - | n < 2 = False - | n < 4 = True - | otherwise = go smallPrimes - where - go (p:ps) - | p*p > n = True - | otherwise = (n `rem` p /= 0) && go ps - go [] = all (isStrongFermatPP n) (take k smallPrimes) +millerRabinV (I# k) n = case testPrimeInteger n k of + 0# -> False + _ -> True --- | @'isStrongFermatPP' n b@ tests whether @n@ is a strong Fermat --- probable prime for base @b@, where @n > 2@ and @1 < b < n@. --- The conditions on the arguments are not checked. +-- | @'isStrongFermatPP' n b@ tests whether non-negative @n@ is +-- a strong Fermat probable prime for base @b@. -- -- Apart from primes, also some composite numbers have the tested -- property, but those are rare. Very rare are composite numbers @@ -88,21 +76,19 @@ -- @1/4@, so five to ten tests give a reasonable level of certainty -- in general. -- --- Some notes about the choice of bases: @b@ is a strong Fermat base --- for @n@ if and only if @n-b@ is, hence one needs only test @b <= (n-1)/2@. --- If @b@ is a strong Fermat base for @n@, then so is @b^k `mod` n@ for --- all @k > 1@, hence one needs not test perfect powers, since their --- base yields a stronger condition. Finally, if @a@ and @b@ are strong --- Fermat bases for @n@, then @a*b@ is in most cases a strong Fermat --- base for @n@, it can only fail to be so if @n `mod` 4 == 1@ and --- the strong Fermat condition is reached at the same step for @a@ as for @b@, --- so primes are the most powerful bases to test. +-- Please consult +-- for the best choice of bases. isStrongFermatPP :: Integer -> Integer -> Bool -isStrongFermatPP n b = a == 1 || go t a +isStrongFermatPP n b + | n < 0 = error "isStrongFermatPP: negative argument" + | n <= 1 = False + | n == 2 = True + | b `mod` n == 0 = True + | otherwise = a == 1 || go t a where m = n-1 (t,u) = shiftToOddCount m - a = powerModInteger' b u n + a = powerModInteger' (b `mod` n) u n go 0 _ = False go k x = x == m || go (k-1) ((x*x) `rem` n) @@ -126,21 +112,21 @@ isFermatPP n b = powerModInteger' b (n-1) n == 1 -- | Primality test after Baillie, Pomerance, Selfridge and Wagstaff. --- The Baillie PSW test consists of a strong Fermat probable primality +-- The Baillie-PSW test consists of a strong Fermat probable primality -- test followed by a (strong) Lucas primality test. This implementation -- assumes that the number @n@ to test is odd and larger than @3@. -- Even and small numbers have to be handled before. Also, before -- applying this test, trial division by small primes should be performed --- to identify many composites cheaply (although the Baillie PSW test is +-- to identify many composites cheaply (although the Baillie-PSW test is -- rather fast, about the same speed as a strong Fermat test for four or -- five bases usually, it is, for large numbers, much more costly than -- trial division by small primes, the primes less than @1000@, say, so -- eliminating numbers with small prime factors beforehand is more efficient). -- --- The Baillie PSW test is very reliable, so far no composite numbers +-- The Baillie-PSW test is very reliable, so far no composite numbers -- passing it are known, and it is known (Gilchrist 2010) that no --- Baillie PSW pseudoprimes exist below @2^64@. However, a heuristic argument --- by Pomerance indicates that there are likely infinitely many Baillie PSW +-- Baillie-PSW pseudoprimes exist below @2^64@. However, a heuristic argument +-- by Pomerance indicates that there are likely infinitely many Baillie-PSW -- pseudoprimes. On the other hand, according to -- there is -- reason to believe that there are none with less than several @@ -229,11 +215,3 @@ -- Listed as a precondition of lucasTest testLucas _ _ _ = error "lucasTest: negative argument" #endif - -smallPrimes :: [Integer] -smallPrimes = 2:3:5:prs - where - prs = 7:11:filter isPr (takeWhile (< 1200) . scanl (+) 13 $ cycle [4,2,4,6,2,6,4,2]) - isPr n = td n prs - td n (p:ps) = (p*p > n) || (n `rem` p /= 0 && td n ps) - td _ [] = True diff --git a/Math/NumberTheory/Recurrencies/Bilinear.hs b/Math/NumberTheory/Recurrencies/Bilinear.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Recurrencies/Bilinear.hs @@ -0,0 +1,197 @@ +-- | +-- Module: Math.NumberTheory.Recurrencies.Bilinear +-- Copyright: (c) 2016 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- Stability: Provisional +-- Portability: Non-portable (GHC extensions) +-- +-- Bilinear recurrent sequences and Bernoulli numbers, +-- roughly covering Ch. 5-6 of /Concrete Mathematics/ +-- by R. L. Graham, D. E. Knuth and O. Patashnik. +-- +-- #memory# __Note on memory leaks and memoization.__ +-- Top-level definitions in this module are polymorphic, so the results of computations are not retained in memory. +-- Make them monomorphic to take advantages of memoization. Compare +-- +-- > > :set +s +-- > > binomial !! 1000 !! 1000 :: Integer +-- > 1 +-- > (0.01 secs, 1,385,512 bytes) +-- > > binomial !! 1000 !! 1000 :: Integer +-- > 1 +-- > (0.01 secs, 1,381,616 bytes) +-- +-- against +-- +-- > > let binomial' = binomial :: [[Integer]] +-- > > binomial' !! 1000 !! 1000 :: Integer +-- > 1 +-- > (0.01 secs, 1,381,696 bytes) +-- > > binomial' !! 1000 !! 1000 :: Integer +-- > 1 +-- > (0.01 secs, 391,152 bytes) + +{-# LANGUAGE CPP #-} + +module Math.NumberTheory.Recurrencies.Bilinear + ( binomial + , stirling1 + , stirling2 + , lah + , eulerian1 + , eulerian2 + , bernoulli + ) where + +import Data.List +import Data.Ratio +import Numeric.Natural + +#if MIN_VERSION_base(4,8,0) +#else +import Data.Word +#endif + +import Math.NumberTheory.Recurrencies.Linear (factorial) + +-- | Infinite zero-based table of binomial coefficients (also known as Pascal triangle): +-- @binomial !! n !! k == n! \/ k! \/ (n - k)!@. +-- +-- > > take 5 (map (take 5) binomial) +-- > [[1],[1,1],[1,2,1],[1,3,3,1],[1,4,6,4,1]] +-- +-- Complexity: @binomial !! n !! k@ is O(n) bits long, its computation +-- takes O(k n) time and forces thunks @binomial !! n !! i@ for @0 <= i <= k@. +-- Use the symmetry of Pascal triangle @binomial !! n !! k == binomial !! n !! (n - k)@ to speed up computations. +-- +-- One could also consider 'Math.Combinat.Numbers.binomial' to compute stand-alone values. +binomial :: Integral a => [[a]] +binomial = map f [0..] + where + f n = scanl (\x k -> x * (n - k + 1) `div` k) 1 [1..n] +{-# SPECIALIZE binomial :: [[Int]] #-} +{-# SPECIALIZE binomial :: [[Word]] #-} +{-# SPECIALIZE binomial :: [[Integer]] #-} +{-# SPECIALIZE binomial :: [[Natural]] #-} + +-- | Infinite zero-based table of . +-- +-- > > take 5 (map (take 5) stirling1) +-- > [[1],[0,1],[0,1,1],[0,2,3,1],[0,6,11,6,1]] +-- +-- Complexity: @stirling1 !! n !! k@ is O(n ln n) bits long, its computation +-- takes O(k n^2 ln n) time and forces thunks @stirling1 !! i !! j@ for @0 <= i <= n@ and @max(0, k - n + i) <= j <= k@. +-- +-- One could also consider 'Math.Combinat.Numbers.unsignedStirling1st' to compute stand-alone values. +stirling1 :: (Num a, Enum a) => [[a]] +stirling1 = scanl f [1] [0..] + where + f xs n = 0 : zipIndexedListWithTail (\_ x y -> x + n * y) 1 xs 0 +{-# SPECIALIZE stirling1 :: [[Int]] #-} +{-# SPECIALIZE stirling1 :: [[Word]] #-} +{-# SPECIALIZE stirling1 :: [[Integer]] #-} +{-# SPECIALIZE stirling1 :: [[Natural]] #-} + +-- | Infinite zero-based table of . +-- +-- > > take 5 (map (take 5) stirling2) +-- > [[1],[0,1],[0,1,1],[0,1,3,1],[0,1,7,6,1]] +-- +-- Complexity: @stirling2 !! n !! k@ is O(n ln n) bits long, its computation +-- takes O(k n^2 ln n) time and forces thunks @stirling2 !! i !! j@ for @0 <= i <= n@ and @max(0, k - n + i) <= j <= k@. +-- +-- One could also consider 'Math.Combinat.Numbers.stirling2nd' to compute stand-alone values. +stirling2 :: (Num a, Enum a) => [[a]] +stirling2 = iterate f [1] + where + f xs = 0 : zipIndexedListWithTail (\k x y -> x + k * y) 1 xs 0 +{-# SPECIALIZE stirling2 :: [[Int]] #-} +{-# SPECIALIZE stirling2 :: [[Word]] #-} +{-# SPECIALIZE stirling2 :: [[Integer]] #-} +{-# SPECIALIZE stirling2 :: [[Natural]] #-} + +-- | Infinite one-based table of . +-- @lah !! n !! k@ equals to lah(n + 1, k + 1). +-- +-- > > take 5 (map (take 5) lah) +-- > [[1],[2,1],[6,6,1],[24,36,12,1],[120,240,120,20,1]] +-- +-- Complexity: @lah !! n !! k@ is O(n ln n) bits long, its computation +-- takes O(k n ln n) time and forces thunks @lah !! n !! i@ for @0 <= i <= k@. +lah :: Integral a => [[a]] +-- Implementation was derived from code by https://github.com/grandpascorpion +lah = zipWith f (tail factorial) [1..] + where + f nf n = scanl (\x k -> x * (n - k) `div` (k * (k + 1))) nf [1..n-1] +{-# SPECIALIZE lah :: [[Int]] #-} +{-# SPECIALIZE lah :: [[Word]] #-} +{-# SPECIALIZE lah :: [[Integer]] #-} +{-# SPECIALIZE lah :: [[Natural]] #-} + +-- | Infinite zero-based table of . +-- +-- > > take 5 (map (take 5) eulerian1) +-- > [[],[1],[1,1],[1,4,1],[1,11,11,1]] +-- +-- Complexity: @eulerian1 !! n !! k@ is O(n ln n) bits long, its computation +-- takes O(k n^2 ln n) time and forces thunks @eulerian1 !! i !! j@ for @0 <= i <= n@ and @max(0, k - n + i) <= j <= k@. +-- +eulerian1 :: (Num a, Enum a) => [[a]] +eulerian1 = scanl f [] [1..] + where + f xs n = 1 : zipIndexedListWithTail (\k x y -> (n - k) * x + (k + 1) * y) 1 xs 0 +{-# SPECIALIZE eulerian1 :: [[Int]] #-} +{-# SPECIALIZE eulerian1 :: [[Word]] #-} +{-# SPECIALIZE eulerian1 :: [[Integer]] #-} +{-# SPECIALIZE eulerian1 :: [[Natural]] #-} + +-- | Infinite zero-based table of . +-- +-- > > take 5 (map (take 5) eulerian2) +-- > [[],[1],[1,2],[1,8,6],[1,22,58,24]] +-- +-- Complexity: @eulerian2 !! n !! k@ is O(n ln n) bits long, its computation +-- takes O(k n^2 ln n) time and forces thunks @eulerian2 !! i !! j@ for @0 <= i <= n@ and @max(0, k - n + i) <= j <= k@. +-- +eulerian2 :: (Num a, Enum a) => [[a]] +eulerian2 = scanl f [] [1..] + where + f xs n = 1 : zipIndexedListWithTail (\k x y -> (2 * n - k - 1) * x + (k + 1) * y) 1 xs 0 +{-# SPECIALIZE eulerian2 :: [[Int]] #-} +{-# SPECIALIZE eulerian2 :: [[Word]] #-} +{-# SPECIALIZE eulerian2 :: [[Integer]] #-} +{-# SPECIALIZE eulerian2 :: [[Natural]] #-} + +-- | Infinite zero-based sequence of , +-- computed via +-- with 'stirling2'. +-- +-- > > take 5 bernoulli +-- > [1 % 1,(-1) % 2,1 % 6,0 % 1,(-1) % 30] +-- +-- Complexity: @bernoulli !! n@ is O(n ln n) bits long, its computation +-- takes O(n^3 ln n) time and forces thunks @stirling2 !! i !! j@ for @0 <= i <= n@ and @0 <= j <= i@. +-- +-- One could also consider 'Math.Combinat.Numbers.bernoulli' to compute stand-alone values. +bernoulli :: Integral a => [Ratio a] +bernoulli = map f stirling2 + where + f = sum . zipWith4 (\sgn denom fact stir -> sgn * fact * stir % denom) (cycle [1, -1]) [1..] factorial +{-# SPECIALIZE bernoulli :: [Ratio Int] #-} +{-# SPECIALIZE bernoulli :: [Rational] #-} + +------------------------------------------------------------------------------- +-- Utils + +-- zipIndexedListWithTail f n as a == zipWith3 f [n..] as (tail as ++ [a]) +-- but inlines much better and avoids checks for distinct sizes of lists. +zipIndexedListWithTail :: Enum b => (b -> a -> a -> b) -> b -> [a] -> a -> [b] +zipIndexedListWithTail f n as a = case as of + [] -> [] + (x : xs) -> go n x xs + where + go m y ys = case ys of + [] -> let v = f m y a in [v] + (z : zs) -> let v = f m y z in (v : go (succ m) z zs) +{-# INLINE zipIndexedListWithTail #-} diff --git a/Math/NumberTheory/Recurrencies/Linear.hs b/Math/NumberTheory/Recurrencies/Linear.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Recurrencies/Linear.hs @@ -0,0 +1,121 @@ +-- | +-- Module: Math.NumberTheory.Recurrencies.Linear +-- Copyright: (c) 2011 Daniel Fischer +-- Licence: MIT +-- Maintainer: Daniel Fischer +-- Stability: Provisional +-- Portability: Non-portable (GHC extensions) +-- +-- Efficient calculation of linear recurrent sequences, including Fibonacci and Lucas sequences. + +{-# LANGUAGE CPP #-} +module Math.NumberTheory.Recurrencies.Linear + ( factorial + , fibonacci + , fibonacciPair + , lucas + , lucasPair + , generalLucas + ) where + +#include "MachDeps.h" + +import Data.Bits +import Numeric.Natural + +#if MIN_VERSION_base(4,8,0) +#else +import Data.Word +#endif + +-- | Infinite zero-based table of factorials. +-- +-- > > take 5 factorial +-- > [1,1,2,6,24] +-- +-- The time-and-space behaviour of 'factorial' is similar to described in +-- "Math.NumberTheory.Recurrencies.Bilinear#memory". +factorial :: (Num a, Enum a) => [a] +factorial = scanl (*) 1 [1..] +{-# SPECIALIZE factorial :: [Int] #-} +{-# SPECIALIZE factorial :: [Word] #-} +{-# SPECIALIZE factorial :: [Integer] #-} +{-# SPECIALIZE factorial :: [Natural] #-} + +-- | @'fibonacci' k@ calculates the @k@-th Fibonacci number in +-- /O/(@log (abs k)@) steps. The index may be negative. This +-- is efficient for calculating single Fibonacci numbers (with +-- large index), but for computing many Fibonacci numbers in +-- close proximity, it is better to use the simple addition +-- formula starting from an appropriate pair of successive +-- Fibonacci numbers. +fibonacci :: Int -> Integer +fibonacci = fst . fibonacciPair + +-- | @'fibonacciPair' k@ returns the pair @(F(k), F(k+1))@ of the @k@-th +-- Fibonacci number and its successor, thus it can be used to calculate +-- the Fibonacci numbers from some index on without needing to compute +-- the previous. The pair is efficiently calculated +-- in /O/(@log (abs k)@) steps. The index may be negative. +fibonacciPair :: Int -> (Integer, Integer) +fibonacciPair n + | n < 0 = let (f,g) = fibonacciPair (-(n+1)) in if testBit n 0 then (g, -f) else (-g, f) + | n == 0 = (0, 1) + | otherwise = look (WORD_SIZE_IN_BITS - 2) + where + look k + | testBit n k = go (k-1) 0 1 + | otherwise = look (k-1) + go k g f + | k < 0 = (f, f+g) + | testBit n k = go (k-1) (f*(f+shiftL g 1)) ((f+g)*shiftL f 1 + g*g) + | otherwise = go (k-1) (f*f+g*g) (f*(f+shiftL g 1)) + +-- | @'lucas' k@ computes the @k@-th Lucas number. Very similar +-- to @'fibonacci'@. +lucas :: Int -> Integer +lucas = fst . lucasPair + +-- | @'lucasPair' k@ computes the pair @(L(k), L(k+1))@ of the @k@-th +-- Lucas number and its successor. Very similar to @'fibonacciPair'@. +lucasPair :: Int -> (Integer, Integer) +lucasPair n + | n < 0 = let (f,g) = lucasPair (-(n+1)) in if testBit n 0 then (-g, f) else (g, -f) + | n == 0 = (2, 1) + | otherwise = look (WORD_SIZE_IN_BITS - 2) + where + look k + | testBit n k = go (k-1) 0 1 + | otherwise = look (k-1) + go k g f + | k < 0 = (shiftL g 1 + f,g+3*f) + | otherwise = go (k-1) g' f' + where + (f',g') + | testBit n k = (shiftL (f*(f+g)) 1 + g*g,f*(shiftL g 1 + f)) + | otherwise = (f*(shiftL g 1 + f),f*f+g*g) + + +-- | @'generalLucas' p q k@ calculates the quadruple @(U(k), U(k+1), V(k), V(k+1))@ +-- where @U(i)@ is the Lucas sequence of the first kind and @V(i)@ the Lucas +-- sequence of the second kind for the parameters @p@ and @q@, where @p^2-4q /= 0@. +-- Both sequences satisfy the recurrence relation @A(j+2) = p*A(j+1) - q*A(j)@, +-- the starting values are @U(0) = 0, U(1) = 1@ and @V(0) = 2, V(1) = p@. +-- The Fibonacci numbers form the Lucas sequence of the first kind for the +-- parameters @p = 1, q = -1@ and the Lucas numbers form the Lucas sequence of +-- the second kind for these parameters. +-- Here, the index must be non-negative, since the terms of the sequence for +-- negative indices are in general not integers. +generalLucas :: Integer -> Integer -> Int -> (Integer, Integer, Integer, Integer) +generalLucas p q k + | k < 0 = error "generalLucas: negative index" + | k == 0 = (0,1,2,p) + | otherwise = look (WORD_SIZE_IN_BITS - 2) + where + look i + | testBit k i = go (i-1) 1 p p q + | otherwise = look (i-1) + go i un un1 vn qn + | i < 0 = (un, un1, vn, p*un1 - shiftL (q*un) 1) + | testBit k i = go (i-1) (un1*vn-qn) ((p*un1-q*un)*vn - p*qn) ((p*un1 - (2*q)*un)*vn - p*qn) (qn*qn*q) + | otherwise = go (i-1) (un*vn) (un1*vn-qn) (vn*vn - 2*qn) (qn*qn) diff --git a/Math/NumberTheory/UniqueFactorisation.hs b/Math/NumberTheory/UniqueFactorisation.hs --- a/Math/NumberTheory/UniqueFactorisation.hs +++ b/Math/NumberTheory/UniqueFactorisation.hs @@ -18,6 +18,7 @@ ) where import Control.Arrow +import Data.Coerce #if MIN_VERSION_base(4,8,0) #else @@ -28,15 +29,6 @@ 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) diff --git a/Math/NumberTheory/Utils.hs b/Math/NumberTheory/Utils.hs --- a/Math/NumberTheory/Utils.hs +++ b/Math/NumberTheory/Utils.hs @@ -20,9 +20,6 @@ , bitCountWord# , uncheckedShiftR , splitOff -#if __GLASGOW_HASKELL__ < 707 - , isTrue# -#endif ) where #include "MachDeps.h" @@ -205,9 +202,3 @@ go !k m = case m `quotRem` p of (q,r) | r == 0 -> go (k+1) q | otherwise -> (k,m) - -#if __GLASGOW_HASKELL__ < 707 --- The times they are a-changing. The types of primops too :( -isTrue# :: Bool -> Bool -isTrue# = id -#endif diff --git a/Math/NumberTheory/Zeta.hs b/Math/NumberTheory/Zeta.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Zeta.hs @@ -0,0 +1,115 @@ +-- | +-- Module: Math.NumberTheory.Zeta +-- Copyright: (c) 2016 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- Stability: Provisional +-- Portability: Non-portable (GHC extensions) +-- +-- Riemann zeta-function. + +{-# LANGUAGE ScopedTypeVariables #-} + +module Math.NumberTheory.Zeta + ( zetas + , zetasEven + , approximateValue + ) where + +import Data.ExactPi +import Data.Ratio + +import Math.NumberTheory.Recurrencies.Bilinear (bernoulli) +import Math.NumberTheory.Recurrencies.Linear (factorial) + +-- | Infinite sequence of exact values of Riemann zeta-function at even arguments, starting with @ζ(0)@. +-- Note that due to numerical errors convertation to 'Double' may return values below 1: +-- +-- > > approximateValue (zetasEven !! 25) :: Double +-- > 0.9999999999999996 +-- +-- Use your favorite type for long-precision arithmetic. For instance, 'Data.Number.Fixed.Fixed' works fine: +-- +-- > > approximateValue (zetasEven !! 25) :: Fixed Prec50 +-- > 1.00000000000000088817842111574532859293035196051773 +-- +zetasEven :: [ExactPi] +zetasEven = zipWith Exact [0, 2 ..] $ zipWith (*) (skipOdds bernoulli) cs + where + cs = (- 1 % 2) : zipWith (\i f -> i * (-4) / fromInteger (2 * f * (2 * f - 1))) cs [1..] + +skipOdds :: [a] -> [a] +skipOdds (x : _ : xs) = x : skipOdds xs +skipOdds xs = xs + +zetasEven' :: Floating a => [a] +zetasEven' = map approximateValue zetasEven + +zetasOdd :: forall a. (Floating a, Ord a) => a -> [a] +zetasOdd eps = (1 / 0) : zets + where + zets :: [a] -- [zeta(3), zeta(5), zeta(7)...] + zets = zipWith (*) zs (tail (iterate (* (- pi * pi)) 1)) + + zs :: [a] -- [zeta(3) / (-pi^2), zeta(5) / pi^4, zeta(7) / (-pi^6)...] + zs = zipWith (\w f -> negate (w / (1 + f))) ws fourth + + ys :: [a] -- [(1 - 1/4) * zeta(3) / (-pi^2), (1 - 1/4^2) * zeta(5) / pi^4...] + ys = zipWith (*) zs fourth + yss :: [[a]] -- [[], [ys !! 0], [ys !! 1, ys !! 0], [ys !! 2, ys !! 1, ys !! 0]...] + yss = scanl (flip (:)) [] ys + + xs :: [a] -- first summand of RHS in (57) for m=[1..] + xs = map (sum . zipWith (flip (/)) factorial2) yss + + ws :: [a] -- RHS in (57) for m=[1..] + ws = zipWith (+) xs cs + + rs :: [a] -- [1, 1/2, 1/3, 1/4...] + rs = map (\n -> recip (fromInteger n)) [1..] + rss :: [[a]] -- [[1, 1/2, 1/3...], [1/2, 1/3, 1/4...], [1/3, 1/4...]] + rss = iterate tail rs + + factorial2 :: [a] -- [2!, 4!, 6!..] + factorial2 = map fromInteger $ tail $ skipOdds factorial + + fourth :: [a] -- [1 - 1/4, 1 - 1/4^2, 1 - 1/4^3...] + fourth = tail $ map (1 -) $ iterate (/ 4) 1 + + as :: [a] -- [zeta(0), zeta(2)/4, zeta(2*2)/4^2, zeta(2*3)/4^3...] + as = zipWith (/) zetasEven' (iterate (* 4) 1) + + bs :: [a] -- map (+ log 2) [b(1), b(2), b(3)...], + -- where b(m) = \sum_{n=0}^\infty zeta(2n) / 4^n / (n + m) + bs = map ((+ log 2) . suminf eps . zipWith (*) as) rss + + cs :: [a] -- second summand of RHS in (57) for m = [1..] + cs = zipWith (\b f -> b / f) bs factorial2 + +suminf :: (Floating a, Ord a) => a -> [a] -> a +suminf eps = sum . takeWhile ((>= eps / 111) . abs) + +-- | Infinite sequence of approximate (up to given precision) +-- values of Riemann zeta-function at integer arguments, starting with @ζ(0)@. +-- Computations for odd arguments are performed in accordance to +-- +-- by J. M. Borwein, D. M. Bradley, R. E. Crandall, formula (57). +-- +-- > > take 5 (zetas 1e-14) :: [Double] +-- > [-0.5,Infinity,1.6449340668482262,1.2020569031595942,1.0823232337111381] +-- +-- Beware to force evaluation of @zetas !! 1@, if the type @a@ does not support infinite values +-- (for instance, 'Data.Number.Fixed.Fixed'). +-- +zetas :: (Floating a, Ord a) => a -> [a] +zetas eps = e : o : scanl1 f (intertwine es os) + where + e : es = zetasEven' + o : os = zetasOdd eps + + intertwine (x : xs) (y : ys) = x : y : intertwine xs ys + intertwine xs ys = xs ++ ys + + -- Cap-and-floor to improve numerical stability: + -- 0 < zeta(n + 1) - 1 < (zeta(n) - 1) / 2 + f x y = 1 `max` (y `min` (1 + (x - 1) / 2)) diff --git a/TODO b/TODO deleted file mode 100644 --- a/TODO +++ /dev/null @@ -1,7 +0,0 @@ -- Atkin sieve -- General number field sieve -- Portability -- Check whether bit twiddling can be as fast as the lookup table for leading and trailing zeros - Using bit twiddling already, faster on my x86_64, not benchmarked on x86 recently, - but it used to be only a marginal difference anyway. -- More Certificates? diff --git a/arithmoi.cabal b/arithmoi.cabal --- a/arithmoi.cabal +++ b/arithmoi.cabal @@ -1,31 +1,30 @@ name : arithmoi -version : 0.4.3.0 +version : 0.5.0.0 cabal-version : >= 1.10 author : Daniel Fischer -copyright : (c) 2011 Daniel Fischer +copyright : (c) 2011 Daniel Fischer, 2016-2017 Andrew Lelechenko, Carter Schonwald license : MIT license-file : LICENSE -maintainer : Carter Schonwald carter at wellposed dot com +maintainer : Carter Schonwald carter at wellposed dot com, + Andrew Lelechenko andrew dot lelechenko at gmail dot com build-type : Simple stability : Provisional homepage : https://github.com/cartazio/arithmoi bug-reports : https://github.com/cartazio/arithmoi/issues synopsis : Efficient basic number-theoretic functions. - Primes, powers, integer logarithms. description : A library of basic functionality needed for number-theoretic calculations. The aim of this library is to provide efficient implementations of the functions. Primes and related things (totients, factorisation), - powers (integer roots and tests, modular exponentiation), - integer logarithms. + powers (integer roots and tests, modular exponentiation). category : Math, Algorithms, Number Theory -tested-with : GHC==7.6.3, GHC==7.8.4, GHC==7.10.3, GHC==8.0.1 +tested-with : GHC==7.8.4, GHC==7.10.3, GHC==8.0.2 -extra-source-files : Changes, TODO +extra-source-files : Changes flag check-bounds description : Replace unsafe array operations with safe ones @@ -34,13 +33,15 @@ library default-language: Haskell2010 - build-depends : base >= 4.6 && < 5 + build-depends : base >= 4.7 && < 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 + , exact-pi >= 0.4.1.1 + , integer-logarithms >= 1.0 if impl(ghc < 7.10) build-depends : nats >= 1 && <1.2 if impl(ghc < 8.0) @@ -49,11 +50,11 @@ 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 - Math.NumberTheory.Lucas + Math.NumberTheory.Recurrencies.Bilinear + Math.NumberTheory.Recurrencies.Linear Math.NumberTheory.GaussianIntegers Math.NumberTheory.GCD Math.NumberTheory.GCD.LowLevel @@ -73,12 +74,12 @@ Math.NumberTheory.Primes.Testing.Certificates Math.NumberTheory.Primes.Heap Math.NumberTheory.UniqueFactorisation + Math.NumberTheory.Zeta other-modules : Math.NumberTheory.Utils Math.NumberTheory.Unsafe Math.NumberTheory.Primes.Counting.Impl Math.NumberTheory.Primes.Counting.Approximate Math.NumberTheory.Primes.Factorisation.Montgomery - Math.NumberTheory.Primes.Factorisation.Utils Math.NumberTheory.Primes.Factorisation.TrialDivision Math.NumberTheory.Primes.Sieve.Eratosthenes Math.NumberTheory.Primes.Sieve.Indexing @@ -104,8 +105,13 @@ , criterion , containers , random + , integer-logarithms + if impl(ghc < 7.10) + build-depends : nats >= 1 && <1.2 other-modules: Math.NumberTheory.ArithmeticFunctionsBench , Math.NumberTheory.PowersBench + , Math.NumberTheory.PrimesBench + , Math.NumberTheory.RecurrenciesBench hs-source-dirs: benchmark main-is: Bench.hs type: exitcode-stdio-1.0 @@ -119,15 +125,15 @@ default-language: Haskell2010 build-depends: base >= 4.6 && < 5 , containers >= 0.5 && < 0.6 - , arithmoi >= 0.4 && < 0.5 + , arithmoi >= 0.5 && < 0.6 , 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 + , transformers >= 0.5 + , integer-gmp < 1.1 if impl(ghc < 7.10) build-depends : nats >= 1 && <1.2 @@ -135,21 +141,23 @@ , Math.NumberTheory.GaussianIntegersTests , Math.NumberTheory.GCDTests , Math.NumberTheory.GCD.LowLevelTests - , Math.NumberTheory.LogarithmsTests - , Math.NumberTheory.LucasTests + , Math.NumberTheory.Recurrencies.LinearTests + , Math.NumberTheory.Recurrencies.BilinearTests , Math.NumberTheory.ModuliTests , Math.NumberTheory.Powers.CubesTests , Math.NumberTheory.MoebiusInversionTests , Math.NumberTheory.MoebiusInversion.IntTests , Math.NumberTheory.Powers.FourthTests , Math.NumberTheory.Powers.GeneralTests - , Math.NumberTheory.Powers.IntegerTests , Math.NumberTheory.Powers.SquaresTests , Math.NumberTheory.PrimesTests , Math.NumberTheory.Primes.CountingTests + , Math.NumberTheory.Primes.FactorisationTests , Math.NumberTheory.Primes.HeapTests , Math.NumberTheory.Primes.SieveTests + , Math.NumberTheory.Primes.TestingTests , Math.NumberTheory.TestUtils , Math.NumberTheory.TestUtils.Wrappers , Math.NumberTheory.TestUtils.Compose , Math.NumberTheory.UniqueFactorisationTests + , Math.NumberTheory.ZetaTests diff --git a/benchmark/Bench.hs b/benchmark/Bench.hs --- a/benchmark/Bench.hs +++ b/benchmark/Bench.hs @@ -1,13 +1,15 @@ -{-# OPTIONS_GHC -fno-warn-warnings-deprecations #-} - module Main where import Criterion.Main import Math.NumberTheory.ArithmeticFunctionsBench as ArithmeticFunctions import Math.NumberTheory.PowersBench as Powers +import Math.NumberTheory.PrimesBench as Primes +import Math.NumberTheory.RecurrenciesBench as Recurrencies main = defaultMain [ ArithmeticFunctions.benchSuite , Powers.benchSuite + , Primes.benchSuite + , Recurrencies.benchSuite ] diff --git a/benchmark/Math/NumberTheory/ArithmeticFunctionsBench.hs b/benchmark/Math/NumberTheory/ArithmeticFunctionsBench.hs --- a/benchmark/Math/NumberTheory/ArithmeticFunctionsBench.hs +++ b/benchmark/Math/NumberTheory/ArithmeticFunctionsBench.hs @@ -1,5 +1,3 @@ -{-# OPTIONS_GHC -fno-warn-warnings-deprecations #-} - module Math.NumberTheory.ArithmeticFunctionsBench ( benchSuite ) where @@ -8,27 +6,20 @@ 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] - ] +compareFunctions :: String -> (Integer -> Integer) -> Benchmark +compareFunctions name new = bench name $ 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] - ] +compareSetFunctions :: String -> (Integer -> Set Integer) -> Benchmark +compareSetFunctions name new = bench name $ nf (map new) [1..100000] benchSuite = bgroup "ArithmeticFunctions" - [ compareSetFunctions "divisors" F.divisors A.divisors + [ compareSetFunctions "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) + , compareFunctions "totient" A.totient + , compareFunctions "carmichael" A.carmichael + , compareFunctions "moebius" A.moebius + , compareFunctions "tau" A.tau + , compareFunctions "sigma 1" (A.sigma 1) + , compareFunctions "sigma 2" (A.sigma 2) ] diff --git a/benchmark/Math/NumberTheory/PowersBench.hs b/benchmark/Math/NumberTheory/PowersBench.hs --- a/benchmark/Math/NumberTheory/PowersBench.hs +++ b/benchmark/Math/NumberTheory/PowersBench.hs @@ -1,5 +1,3 @@ -{-# OPTIONS_GHC -fno-warn-warnings-deprecations #-} - module Math.NumberTheory.PowersBench ( benchSuite ) where diff --git a/benchmark/Math/NumberTheory/PrimesBench.hs b/benchmark/Math/NumberTheory/PrimesBench.hs new file mode 100644 --- /dev/null +++ b/benchmark/Math/NumberTheory/PrimesBench.hs @@ -0,0 +1,32 @@ +{-# OPTIONS_GHC -fno-warn-warnings-deprecations #-} + +module Math.NumberTheory.PrimesBench + ( benchSuite + ) where + +import Criterion.Main +import System.Random + +import Math.NumberTheory.Logarithms (integerLog2) +import Math.NumberTheory.Primes + +genInteger :: Int -> Int -> Integer +genInteger salt bits + = head + . dropWhile ((< bits) . integerLog2) + . scanl (\a r -> a * 2^31 + abs r) 1 + . randoms + . mkStdGen + $ salt + bits + +comparePrimalityTests :: Int -> Benchmark +comparePrimalityTests bits = bgroup ("primality" ++ show bits) + [ bench "isPrime" $ nf (map isPrime) ns + , bench "millerRabinV 0" $ nf (map $ millerRabinV 0) ns + , bench "millerRabinV 10" $ nf (map $ millerRabinV 10) ns + , bench "millerRabinV 50" $ nf (map $ millerRabinV 50) ns + ] + where + ns = take bits [genInteger 0 bits ..] + +benchSuite = bgroup "Primes" $ map comparePrimalityTests [50, 100, 200, 500, 1000, 2000] diff --git a/benchmark/Math/NumberTheory/RecurrenciesBench.hs b/benchmark/Math/NumberTheory/RecurrenciesBench.hs new file mode 100644 --- /dev/null +++ b/benchmark/Math/NumberTheory/RecurrenciesBench.hs @@ -0,0 +1,30 @@ +{-# LANGUAGE RankNTypes #-} + +module Math.NumberTheory.RecurrenciesBench + ( benchSuite + ) where + +import Criterion.Main +import Numeric.Natural +import System.Random + +import Math.NumberTheory.Recurrencies.Bilinear + +benchTriangle :: String -> (forall a. (Integral a) => [[a]]) -> Int -> Benchmark +benchTriangle name triangle n = bgroup name + [ benchAt (10 * n) (1 * n) + , benchAt (10 * n) (2 * n) + , benchAt (10 * n) (5 * n) + , benchAt (10 * n) (9 * n) + ] + where + benchAt i j = bench ("!! " ++ show i ++ " !! " ++ show j) + $ nf (\(x, y) -> triangle !! x !! y :: Integer) (i, j) + +benchSuite = bgroup "Bilinear" + [ benchTriangle "binomial" binomial 1000 + , benchTriangle "stirling1" stirling1 100 + , benchTriangle "stirling2" stirling2 100 + , benchTriangle "eulerian1" eulerian1 100 + , benchTriangle "eulerian2" eulerian2 100 + ] diff --git a/test-suite/Math/NumberTheory/LogarithmsTests.hs b/test-suite/Math/NumberTheory/LogarithmsTests.hs deleted file mode 100644 --- a/test-suite/Math/NumberTheory/LogarithmsTests.hs +++ /dev/null @@ -1,112 +0,0 @@ --- | --- Module: Math.NumberTheory.LogarithmsTests --- Copyright: (c) 2016 Andrew Lelechenko --- Licence: MIT --- Maintainer: Andrew Lelechenko --- Stability: Provisional --- --- Tests for Math.NumberTheory.Logarithms --- - -{-# LANGUAGE CPP #-} - -{-# OPTIONS_GHC -fno-warn-type-defaults #-} - -module Math.NumberTheory.LogarithmsTests - ( testSuite - ) where - -import Test.Tasty - -#if MIN_VERSION_base(4,8,0) -#else -import Data.Word -#endif - -import Math.NumberTheory.Logarithms -import Math.NumberTheory.TestUtils - --- | Check that 'integerLogBase' returns the largest integer @l@ such that @b@ ^ @l@ <= @n@ and @b@ ^ (@l@+1) > @n@. -integerLogBaseProperty :: Positive Integer -> Positive Integer -> Bool -integerLogBaseProperty (Positive b) (Positive n) = b == 1 || b ^ l <= n && b ^ (l + 1) > n - where - l = toInteger $ integerLogBase b n - --- | Check that 'integerLog2' returns the largest integer @l@ such that 2 ^ @l@ <= @n@ and 2 ^ (@l@+1) > @n@. -integerLog2Property :: Positive Integer -> Bool -integerLog2Property (Positive n) = 2 ^ l <= n && 2 ^ (l + 1) > n - where - l = toInteger $ integerLog2 n - --- | Check that 'integerLog10' returns the largest integer @l@ such that 10 ^ @l@ <= @n@ and 10 ^ (@l@+1) > @n@. -integerLog10Property :: Positive Integer -> Bool -integerLog10Property (Positive n) = 10 ^ l <= n && 10 ^ (l + 1) > n - where - l = toInteger $ integerLog10 n - --- | Check that 'intLog2' returns the largest integer @l@ such that 2 ^ @l@ <= @n@ and 2 ^ (@l@+1) > @n@. -intLog2Property :: Positive Int -> Bool -intLog2Property (Positive n) = 2 ^ l <= n && (2 ^ (l + 1) > n || n > maxBound `div` 2) - where - l = intLog2 n - --- | Check that 'wordLog2' returns the largest integer @l@ such that 2 ^ @l@ <= @n@ and 2 ^ (@l@+1) > @n@. -wordLog2Property :: Positive Word -> Bool -wordLog2Property (Positive n) = 2 ^ l <= n && (2 ^ (l + 1) > n || n > maxBound `div` 2) - where - l = wordLog2 n - --- | Check that 'integerLogBase'' returns the largest integer @l@ such that @b@ ^ @l@ <= @n@ and @b@ ^ (@l@+1) > @n@. -integerLogBase'Property :: Positive Integer -> Positive Integer -> Bool -integerLogBase'Property (Positive b) (Positive n) = b == 1 || b ^ l <= n && b ^ (l + 1) > n - where - l = toInteger $ integerLogBase' b n - --- | Check that 'integerLogBase'' returns the largest integer @l@ such that @b@ ^ @l@ <= @n@ and @b@ ^ (@l@+1) > @n@ for @b@ > 32 and @n@ >= @b@ ^ 2. -integerLogBase'Property2 :: Positive Integer -> Positive Integer -> Bool -integerLogBase'Property2 (Positive b') (Positive n') = b ^ l <= n && b ^ (l + 1) > n - where - b = b' + 32 - n = n' + b ^ 2 - 1 - l = toInteger $ integerLogBase' b n - --- | Check that 'integerLog2'' returns the largest integer @l@ such that 2 ^ @l@ <= @n@ and 2 ^ (@l@+1) > @n@. -integerLog2'Property :: Positive Integer -> Bool -integerLog2'Property (Positive n) = 2 ^ l <= n && 2 ^ (l + 1) > n - where - l = toInteger $ integerLog2' n - --- | Check that 'integerLog10'' returns the largest integer @l@ such that 10 ^ @l@ <= @n@ and 10 ^ (@l@+1) > @n@. -integerLog10'Property :: Positive Integer -> Bool -integerLog10'Property (Positive n) = 10 ^ l <= n && 10 ^ (l + 1) > n - where - l = toInteger $ integerLog10' n - --- | Check that 'intLog2'' returns the largest integer @l@ such that 2 ^ @l@ <= @n@ and 2 ^ (@l@+1) > @n@. -intLog2'Property :: Positive Int -> Bool -intLog2'Property (Positive n) = 2 ^ l <= n && (2 ^ (l + 1) > n || n > maxBound `div` 2) - where - l = intLog2' n - --- | Check that 'wordLog2'' returns the largest integer @l@ such that 2 ^ @l@ <= @n@ and 2 ^ (@l@+1) > @n@. -wordLog2'Property :: Positive Word -> Bool -wordLog2'Property (Positive n) = 2 ^ l <= n && (2 ^ (l + 1) > n || n > maxBound `div` 2) - where - l = wordLog2' n - -testSuite :: TestTree -testSuite = testGroup "Logarithms" - [ testSmallAndQuick "integerLogBase" integerLogBaseProperty - , testSmallAndQuick "integerLog2" integerLog2Property - , testSmallAndQuick "integerLog10" integerLog10Property - , testSmallAndQuick "intLog2" intLog2Property - , testSmallAndQuick "wordLog2" wordLog2Property - - , testSmallAndQuick "integerLogBase'" integerLogBase'Property - , testSmallAndQuick "integerLogBase' with base > 32 and n >= base ^ 2" - integerLogBase'Property2 - , testSmallAndQuick "integerLog2'" integerLog2'Property - , testSmallAndQuick "integerLog10'" integerLog10'Property - , testSmallAndQuick "intLog2'" intLog2'Property - , testSmallAndQuick "wordLog2'" wordLog2'Property - ] diff --git a/test-suite/Math/NumberTheory/LucasTests.hs b/test-suite/Math/NumberTheory/LucasTests.hs deleted file mode 100644 --- a/test-suite/Math/NumberTheory/LucasTests.hs +++ /dev/null @@ -1,104 +0,0 @@ --- | --- Module: Math.NumberTheory.LucasTests --- Copyright: (c) 2016 Andrew Lelechenko --- Licence: MIT --- Maintainer: Andrew Lelechenko --- Stability: Provisional --- --- Tests for Math.NumberTheory.Lucas --- - -{-# LANGUAGE CPP #-} - -{-# OPTIONS_GHC -fno-warn-type-defaults #-} - -module Math.NumberTheory.LucasTests - ( testSuite - ) where - -import Test.Tasty -import Test.Tasty.HUnit - -import Math.NumberTheory.Lucas -import Math.NumberTheory.TestUtils - --- | Check that 'fibonacci' matches the definition of Fibonacci sequence. -fibonacciProperty1 :: AnySign Int -> Bool -fibonacciProperty1 (AnySign n) = fibonacci n + fibonacci (n + 1) == fibonacci (n +2) - --- | Check that 'fibonacci' for negative indices is correctly defined. -fibonacciProperty2 :: NonNegative Int -> Bool -fibonacciProperty2 (NonNegative n) = fibonacci n == (if even n then negate else id) (fibonacci (- n)) - --- | Check that 'fibonacciPair' is a pair of consequent 'fibonacci'. -fibonacciPairProperty :: AnySign Int -> Bool -fibonacciPairProperty (AnySign n) = fibonacciPair n == (fibonacci n, fibonacci (n + 1)) - --- | Check that 'fibonacci 0' is 0. -fibonacciSpecialCase0 :: Assertion -fibonacciSpecialCase0 = assertEqual "fibonacci" (fibonacci 0) 0 - --- | Check that 'fibonacci 1' is 1. -fibonacciSpecialCase1 :: Assertion -fibonacciSpecialCase1 = assertEqual "fibonacci" (fibonacci 1) 1 - - --- | Check that 'lucas' matches the definition of Lucas sequence. -lucasProperty1 :: AnySign Int -> Bool -lucasProperty1 (AnySign n) = lucas n + lucas (n + 1) == lucas (n +2) - --- | Check that 'lucas' for negative indices is correctly defined. -lucasProperty2 :: NonNegative Int -> Bool -lucasProperty2 (NonNegative n) = lucas n == (if odd n then negate else id) (lucas (- n)) - --- | Check that 'lucasPair' is a pair of consequent 'lucas'. -lucasPairProperty :: AnySign Int -> Bool -lucasPairProperty (AnySign n) = lucasPair n == (lucas n, lucas (n + 1)) - --- | Check that 'lucas 0' is 2. -lucasSpecialCase0 :: Assertion -lucasSpecialCase0 = assertEqual "lucas" (lucas 0) 2 - --- | Check that 'lucas 1' is 1. -lucasSpecialCase1 :: Assertion -lucasSpecialCase1 = assertEqual "lucas" (lucas 1) 1 - --- | Check that 'generalLucas' matches its definition. -generalLucasProperty1 :: AnySign Integer -> AnySign Integer -> NonNegative Int -> Bool -generalLucasProperty1 (AnySign p) (AnySign q) (NonNegative n) = un1 == un1' && vn1 == vn1' && un2 == p * un1 - q * un && vn2 == p * vn1 - q * vn - where - (un, un1, vn, vn1) = generalLucas p q n - (un1', un2, vn1', vn2) = generalLucas p q (n + 1) - --- | Check that 'generalLucas' 1 (-1) is 'fibonacciPair' plus 'lucasPair'. -generalLucasProperty2 :: NonNegative Int -> Bool -generalLucasProperty2 (NonNegative n) = (un, un1) == fibonacciPair n && (vn, vn1) == lucasPair n - where - (un, un1, vn, vn1) = generalLucas 1 (-1) n - --- | Check that 'generalLucas' p _ 0 is (0, 1, 2, p). -generalLucasProperty3 :: AnySign Integer -> AnySign Integer -> Bool -generalLucasProperty3 (AnySign p) (AnySign q) = generalLucas p q 0 == (0, 1, 2, p) - -testSuite :: TestTree -testSuite = testGroup "Lucas" - [ testGroup "fibonacci" - [ testSmallAndQuick "matches definition" fibonacciProperty1 - , testSmallAndQuick "negative indices" fibonacciProperty2 - , testSmallAndQuick "pair" fibonacciPairProperty - , testCase "fibonacci 0" fibonacciSpecialCase0 - , testCase "fibonacci 1" fibonacciSpecialCase1 - ] - , testGroup "lucas" - [ testSmallAndQuick "matches definition" lucasProperty1 - , testSmallAndQuick "negative indices" lucasProperty2 - , testSmallAndQuick "pair" lucasPairProperty - , testCase "lucas 0" lucasSpecialCase0 - , testCase "lucas 1" lucasSpecialCase1 - ] - , testGroup "generalLucas" - [ testSmallAndQuick "matches definition" generalLucasProperty1 - , testSmallAndQuick "generalLucas 1 (-1)" generalLucasProperty2 - , testSmallAndQuick "generalLucas _ _ 0" generalLucasProperty3 - ] - ] diff --git a/test-suite/Math/NumberTheory/Powers/CubesTests.hs b/test-suite/Math/NumberTheory/Powers/CubesTests.hs --- a/test-suite/Math/NumberTheory/Powers/CubesTests.hs +++ b/test-suite/Math/NumberTheory/Powers/CubesTests.hs @@ -142,9 +142,11 @@ , testSmallAndQuick "almost cube Int" integerCubeRootProperty2_Int , testSmallAndQuick "almost cube Word" integerCubeRootProperty2_Word +#if WORD_SIZE_IN_BITS == 64 , testCase "maxBound :: Int" integerCubeRootSpecialCase1_Int , testCase "maxBound / 2 :: Word" integerCubeRootSpecialCase1_Word , testCase "maxBound :: Word" integerCubeRootSpecialCase2 +#endif ] , testIntegralProperty "integerCubeRoot'" integerCubeRoot'Property , testIntegralProperty "isCube" isCubeProperty diff --git a/test-suite/Math/NumberTheory/Powers/FourthTests.hs b/test-suite/Math/NumberTheory/Powers/FourthTests.hs --- a/test-suite/Math/NumberTheory/Powers/FourthTests.hs +++ b/test-suite/Math/NumberTheory/Powers/FourthTests.hs @@ -133,9 +133,11 @@ , testSmallAndQuick "almost Fourth Int" integerFourthRootProperty2_Int , testSmallAndQuick "almost Fourth Word" integerFourthRootProperty2_Word +#if WORD_SIZE_IN_BITS == 64 , testCase "maxBound / 8 :: Int" integerFourthRootSpecialCase1_Int , testCase "maxBound / 16 :: Word" integerFourthRootSpecialCase1_Word , testCase "maxBound :: Word" integerFourthRootSpecialCase2 +#endif ] , testIntegralProperty "integerFourthRoot'" integerFourthRoot'Property , testIntegralProperty "isFourthPower" isFourthPowerProperty diff --git a/test-suite/Math/NumberTheory/Powers/IntegerTests.hs b/test-suite/Math/NumberTheory/Powers/IntegerTests.hs deleted file mode 100644 --- a/test-suite/Math/NumberTheory/Powers/IntegerTests.hs +++ /dev/null @@ -1,41 +0,0 @@ --- | --- Module: Math.NumberTheory.Powers.IntegerTests --- Copyright: (c) 2016 Andrew Lelechenko --- Licence: MIT --- Maintainer: Andrew Lelechenko --- Stability: Provisional --- --- Tests for Math.NumberTheory.Powers.Integer --- - -{-# LANGUAGE CPP #-} - -{-# OPTIONS_GHC -fno-warn-type-defaults #-} - -module Math.NumberTheory.Powers.IntegerTests - ( testSuite - ) where - -import Test.Tasty - -#if MIN_VERSION_base(4,8,0) -#else -import Data.Word -#endif - -import Math.NumberTheory.Powers.Integer -import Math.NumberTheory.TestUtils - --- | Check that 'integerPower' == '^'. -integerPowerProperty :: Integer -> Power Int -> Bool -integerPowerProperty a (Power b) = integerPower a b == a ^ b - --- | Check that 'integerWordPower' == '^'. -integerWordPowerProperty :: Integer -> Power Word -> Bool -integerWordPowerProperty a (Power b) = integerWordPower a b == a ^ b - -testSuite :: TestTree -testSuite = testGroup "Integer" - [ testSmallAndQuick "integerPower" integerPowerProperty - , testSmallAndQuick "integerWordPower" integerWordPowerProperty - ] diff --git a/test-suite/Math/NumberTheory/Primes/FactorisationTests.hs b/test-suite/Math/NumberTheory/Primes/FactorisationTests.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/NumberTheory/Primes/FactorisationTests.hs @@ -0,0 +1,44 @@ +-- | +-- Module: Math.NumberTheory.Primes.FactorisationTests +-- Copyright: (c) 2017 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- Stability: Provisional +-- +-- Tests for Math.NumberTheory.Primes.Factorisation +-- + +{-# OPTIONS_GHC -fno-warn-type-defaults #-} + +module Math.NumberTheory.Primes.FactorisationTests + ( testSuite + ) where + +import Test.Tasty +import Test.Tasty.HUnit + +import Math.NumberTheory.Primes.Factorisation +import Math.NumberTheory.Primes.Testing +import Math.NumberTheory.TestUtils + +factoriseProperty1 :: Assertion +factoriseProperty1 = assertEqual "0" [] (factorise 1) + +factoriseProperty2 :: Positive Integer -> Bool +factoriseProperty2 (Positive n) = (-1, 1) : factorise n == factorise (negate n) + +factoriseProperty3 :: Positive Integer -> Bool +factoriseProperty3 (Positive n) = all (isPrime . fst) (factorise n) + +factoriseProperty4 :: Positive Integer -> Bool +factoriseProperty4 (Positive n) = product (map (uncurry (^)) (factorise n)) == n + +testSuite :: TestTree +testSuite = testGroup "Factorisation" + [ testGroup "factorise" + [ testCase "0" factoriseProperty1 + , testSmallAndQuick "negate" factoriseProperty2 + , testSmallAndQuick "bases are prime" factoriseProperty3 + , testSmallAndQuick "factorback" factoriseProperty4 + ] + ] diff --git a/test-suite/Math/NumberTheory/Primes/TestingTests.hs b/test-suite/Math/NumberTheory/Primes/TestingTests.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/NumberTheory/Primes/TestingTests.hs @@ -0,0 +1,73 @@ +-- | +-- Module: Math.NumberTheory.Primes.TestingTests +-- Copyright: (c) 2017 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- Stability: Provisional +-- +-- Tests for Math.NumberTheory.Primes.Testing +-- + +{-# OPTIONS_GHC -fno-warn-type-defaults #-} + +module Math.NumberTheory.Primes.TestingTests + ( testSuite + ) where + +import Test.Tasty +import Test.Tasty.HUnit + +import GHC.Integer.GMP.Internals (nextPrimeInteger) + +import Math.NumberTheory.Primes.Testing +import Math.NumberTheory.TestUtils + +isPrimeProperty1 :: Assertion +isPrimeProperty1 = assertEqual "[0..100]" expected actual + where + expected = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97] + actual = filter isPrime [0..100] + +isPrimeProperty2 :: Integer -> Bool +isPrimeProperty2 n = isPrime n == isPrime (negate n) + +isPrimeProperty3 :: Assertion +isPrimeProperty3 = assertBool "Carmichael pseudoprimes" $ all (not . isPrime) pseudoprimes + where + -- OEIS A002997 + pseudoprimes = [561, 1105, 1729, 2465, 2821, 6601, 8911, 10585, 15841, 29341, 41041, 46657, 52633, 62745, 63973, 75361, 101101, 115921, 126217, 162401, 172081, 188461, 252601, 278545, 294409, 314821, 334153, 340561, 399001, 410041, 449065, 488881, 512461] + +isPrimeProperty4 :: Assertion +isPrimeProperty4 = assertBool "strong pseudoprimes to base 2" $ all (not . isPrime) pseudoprimes + where + -- OEIS A001262 + pseudoprimes = [2047, 3277, 4033, 4681, 8321, 15841, 29341, 42799, 49141, 52633, 65281, 74665, 80581, 85489, 88357, 90751, 104653, 130561, 196093, 220729, 233017, 252601, 253241, 256999, 271951, 280601, 314821, 357761, 390937, 458989, 476971, 486737] + +isPrimeProperty5 :: Assertion +isPrimeProperty5 = assertBool "strong Lucas pseudoprimes" $ all (not . isPrime) pseudoprimes + where + -- OEIS A217255 + pseudoprimes = [5459, 5777, 10877, 16109, 18971, 22499, 24569, 25199, 40309, 58519, 75077, 97439, 100127, 113573, 115639, 130139, 155819, 158399, 161027, 162133, 176399, 176471, 189419, 192509, 197801, 224369, 230691, 231703, 243629, 253259, 268349, 288919, 313499, 324899] + +isPrimeProperty6 :: NonNegative Integer -> Bool +isPrimeProperty6 (NonNegative n) = if isPrime n + then nextPrimeInteger (n - 1) == n + else isPrime (nextPrimeInteger n) + +isStrongFermatPPProperty :: NonNegative Integer -> Integer -> Bool +isStrongFermatPPProperty (NonNegative n) b = not (isPrime n) || isStrongFermatPP n b + +testSuite :: TestTree +testSuite = testGroup "Testing" + [ testGroup "isPrime" + [ testCase "[0..100]" isPrimeProperty1 + , testSmallAndQuick "negate" isPrimeProperty2 + , testCase "Carmichael pseudoprimes" isPrimeProperty3 + , testCase "strong pseudoprimes base 2" isPrimeProperty4 + , testCase "strong Lucas pseudoprimes" isPrimeProperty5 + , testSmallAndQuick "matches GMP" isPrimeProperty6 + ] + , testGroup "isStrongFermatPP" + [ testSmallAndQuick "matches isPrime" isStrongFermatPPProperty + ] + ] diff --git a/test-suite/Math/NumberTheory/Recurrencies/BilinearTests.hs b/test-suite/Math/NumberTheory/Recurrencies/BilinearTests.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/NumberTheory/Recurrencies/BilinearTests.hs @@ -0,0 +1,196 @@ +-- | +-- Module: Math.NumberTheory.Recurrencies.BilinearTests +-- Copyright: (c) 2016 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- Stability: Provisional +-- +-- Tests for Math.NumberTheory.Recurrencies.Bilinear +-- + +{-# OPTIONS_GHC -fno-warn-type-defaults #-} + +module Math.NumberTheory.Recurrencies.BilinearTests + ( testSuite + ) where + +import Test.Tasty +import Test.Tasty.HUnit + +import Data.Ratio + +import Math.NumberTheory.Recurrencies.Bilinear +import Math.NumberTheory.TestUtils + +binomialProperty1 :: NonNegative Int -> Bool +binomialProperty1 (NonNegative i) = length (binomial !! i) == i + 1 + +binomialProperty2 :: NonNegative Int -> Bool +binomialProperty2 (NonNegative i) = binomial !! i !! 0 == 1 + +binomialProperty3 :: NonNegative Int -> Bool +binomialProperty3 (NonNegative i) = binomial !! i !! i == 1 + +binomialProperty4 :: Positive Int -> Positive Int -> Bool +binomialProperty4 (Positive i) (Positive j) + = j >= i + || binomial !! i !! j + == binomial !! (i - 1) !! (j - 1) + + binomial !! (i - 1) !! j + +stirling1Property1 :: NonNegative Int -> Bool +stirling1Property1 (NonNegative i) = length (stirling1 !! i) == i + 1 + +stirling1Property2 :: NonNegative Int -> Bool +stirling1Property2 (NonNegative i) + = stirling1 !! i !! 0 + == if i == 0 then 1 else 0 + +stirling1Property3 :: NonNegative Int -> Bool +stirling1Property3 (NonNegative i) = stirling1 !! i !! i == 1 + +stirling1Property4 :: Positive Int -> Positive Int -> Bool +stirling1Property4 (Positive i) (Positive j) + = j >= i + || stirling1 !! i !! j + == stirling1 !! (i - 1) !! (j - 1) + + (toInteger i - 1) * stirling1 !! (i - 1) !! j + +stirling2Property1 :: NonNegative Int -> Bool +stirling2Property1 (NonNegative i) = length (stirling2 !! i) == i + 1 + +stirling2Property2 :: NonNegative Int -> Bool +stirling2Property2 (NonNegative i) + = stirling2 !! i !! 0 + == if i == 0 then 1 else 0 + +stirling2Property3 :: NonNegative Int -> Bool +stirling2Property3 (NonNegative i) = stirling2 !! i !! i == 1 + +stirling2Property4 :: Positive Int -> Positive Int -> Bool +stirling2Property4 (Positive i) (Positive j) + = j >= i + || stirling2 !! i !! j + == stirling2 !! (i - 1) !! (j - 1) + + toInteger j * stirling2 !! (i - 1) !! j + +lahProperty1 :: NonNegative Int -> Bool +lahProperty1 (NonNegative i) = length (lah !! i) == i + 1 + +lahProperty2 :: NonNegative Int -> Bool +lahProperty2 (NonNegative i) + = lah !! i !! 0 + == product [1 .. i+1] + +lahProperty3 :: NonNegative Int -> Bool +lahProperty3 (NonNegative i) = lah !! i !! i == 1 + +lahProperty4 :: Positive Int -> Positive Int -> Bool +lahProperty4 (Positive i) (Positive j) + = j >= i + || lah !! i !! j + == sum [ stirling1 !! (i + 1) !! k * stirling2 !! k !! (j + 1) | k <- [j + 1 .. i + 1] ] + +eulerian1Property1 :: NonNegative Int -> Bool +eulerian1Property1 (NonNegative i) = length (eulerian1 !! i) == i + +eulerian1Property2 :: Positive Int -> Bool +eulerian1Property2 (Positive i) = eulerian1 !! i !! 0 == 1 + +eulerian1Property3 :: Positive Int -> Bool +eulerian1Property3 (Positive i) = eulerian1 !! i !! (i - 1) == 1 + +eulerian1Property4 :: Positive Int -> Positive Int -> Bool +eulerian1Property4 (Positive i) (Positive j) + = j >= i - 1 + || eulerian1 !! i !! j + == (toInteger $ i - j) * eulerian1 !! (i - 1) !! (j - 1) + + (toInteger j + 1) * eulerian1 !! (i - 1) !! j + +eulerian2Property1 :: NonNegative Int -> Bool +eulerian2Property1 (NonNegative i) = length (eulerian2 !! i) == i + +eulerian2Property2 :: Positive Int -> Bool +eulerian2Property2 (Positive i) + = eulerian2 !! i !! 0 == 1 + +eulerian2Property3 :: Positive Int -> Bool +eulerian2Property3 (Positive i) + = eulerian2 !! i !! (i - 1) + == product [1 .. toInteger i] + +eulerian2Property4 :: Positive Int -> Positive Int -> Bool +eulerian2Property4 (Positive i) (Positive j) + = j >= i - 1 + || eulerian2 !! i !! j + == (toInteger $ 2 * i - j - 1) * eulerian2 !! (i - 1) !! (j - 1) + + (toInteger j + 1) * eulerian2 !! (i - 1) !! j + +bernoulliSpecialCase1 :: Assertion +bernoulliSpecialCase1 = assertEqual "B_0 = 1" (bernoulli !! 0) 1 + +bernoulliSpecialCase2 :: Assertion +bernoulliSpecialCase2 = assertEqual "B_1 = -1/2" (bernoulli !! 1) (- 1 % 2) + +bernoulliProperty1 :: NonNegative Int -> Bool +bernoulliProperty1 (NonNegative m) + = case signum (bernoulli !! m) of + 1 -> m == 0 || m `mod` 4 == 2 + 0 -> m /= 1 && odd m + -1 -> m == 1 || (m /= 0 && m `mod` 4 == 0) + _ -> False + +bernoulliProperty2 :: NonNegative Int -> Bool +bernoulliProperty2 (NonNegative m) + = bernoulli !! m + == (if m == 0 then 1 else 0) + - sum [ bernoulli !! k + * (binomial !! m !! k % (toInteger $ m - k + 1)) + | k <- [0 .. m - 1] + ] + +testSuite :: TestTree +testSuite = testGroup "Bilinear" + [ testGroup "binomial" + [ testSmallAndQuick "shape" binomialProperty1 + , testSmallAndQuick "left side" binomialProperty2 + , testSmallAndQuick "right side" binomialProperty3 + , testSmallAndQuick "recurrency" binomialProperty4 + ] + , testGroup "stirling1" + [ testSmallAndQuick "shape" stirling1Property1 + , testSmallAndQuick "left side" stirling1Property2 + , testSmallAndQuick "right side" stirling1Property3 + , testSmallAndQuick "recurrency" stirling1Property4 + ] + , testGroup "stirling2" + [ testSmallAndQuick "shape" stirling2Property1 + , testSmallAndQuick "left side" stirling2Property2 + , testSmallAndQuick "right side" stirling2Property3 + , testSmallAndQuick "recurrency" stirling2Property4 + ] + , testGroup "lah" + [ testSmallAndQuick "shape" lahProperty1 + , testSmallAndQuick "left side" lahProperty2 + , testSmallAndQuick "right side" lahProperty3 + , testSmallAndQuick "zip stirlings" lahProperty4 + ] + , testGroup "eulerian1" + [ testSmallAndQuick "shape" eulerian1Property1 + , testSmallAndQuick "left side" eulerian1Property2 + , testSmallAndQuick "right side" eulerian1Property3 + , testSmallAndQuick "recurrency" eulerian1Property4 + ] + , testGroup "eulerian2" + [ testSmallAndQuick "shape" eulerian2Property1 + , testSmallAndQuick "left side" eulerian2Property2 + , testSmallAndQuick "right side" eulerian2Property3 + , testSmallAndQuick "recurrency" eulerian2Property4 + ] + , testGroup "bernoulli" + [ testCase "B_0" bernoulliSpecialCase1 + , testCase "B_1" bernoulliSpecialCase2 + , testSmallAndQuick "sign" bernoulliProperty1 + , testSmallAndQuick "recursive definition" bernoulliProperty2 + ] + ] diff --git a/test-suite/Math/NumberTheory/Recurrencies/LinearTests.hs b/test-suite/Math/NumberTheory/Recurrencies/LinearTests.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/NumberTheory/Recurrencies/LinearTests.hs @@ -0,0 +1,104 @@ +-- | +-- Module: Math.NumberTheory.Recurrencies.LinearTests +-- Copyright: (c) 2016 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- Stability: Provisional +-- +-- Tests for Math.NumberTheory.Recurrencies.Linear +-- + +{-# LANGUAGE CPP #-} + +{-# OPTIONS_GHC -fno-warn-type-defaults #-} + +module Math.NumberTheory.Recurrencies.LinearTests + ( testSuite + ) where + +import Test.Tasty +import Test.Tasty.HUnit + +import Math.NumberTheory.Recurrencies.Linear +import Math.NumberTheory.TestUtils + +-- | Check that 'fibonacci' matches the definition of Fibonacci sequence. +fibonacciProperty1 :: AnySign Int -> Bool +fibonacciProperty1 (AnySign n) = fibonacci n + fibonacci (n + 1) == fibonacci (n +2) + +-- | Check that 'fibonacci' for negative indices is correctly defined. +fibonacciProperty2 :: NonNegative Int -> Bool +fibonacciProperty2 (NonNegative n) = fibonacci n == (if even n then negate else id) (fibonacci (- n)) + +-- | Check that 'fibonacciPair' is a pair of consequent 'fibonacci'. +fibonacciPairProperty :: AnySign Int -> Bool +fibonacciPairProperty (AnySign n) = fibonacciPair n == (fibonacci n, fibonacci (n + 1)) + +-- | Check that 'fibonacci 0' is 0. +fibonacciSpecialCase0 :: Assertion +fibonacciSpecialCase0 = assertEqual "fibonacci" (fibonacci 0) 0 + +-- | Check that 'fibonacci 1' is 1. +fibonacciSpecialCase1 :: Assertion +fibonacciSpecialCase1 = assertEqual "fibonacci" (fibonacci 1) 1 + + +-- | Check that 'lucas' matches the definition of Lucas sequence. +lucasProperty1 :: AnySign Int -> Bool +lucasProperty1 (AnySign n) = lucas n + lucas (n + 1) == lucas (n +2) + +-- | Check that 'lucas' for negative indices is correctly defined. +lucasProperty2 :: NonNegative Int -> Bool +lucasProperty2 (NonNegative n) = lucas n == (if odd n then negate else id) (lucas (- n)) + +-- | Check that 'lucasPair' is a pair of consequent 'lucas'. +lucasPairProperty :: AnySign Int -> Bool +lucasPairProperty (AnySign n) = lucasPair n == (lucas n, lucas (n + 1)) + +-- | Check that 'lucas 0' is 2. +lucasSpecialCase0 :: Assertion +lucasSpecialCase0 = assertEqual "lucas" (lucas 0) 2 + +-- | Check that 'lucas 1' is 1. +lucasSpecialCase1 :: Assertion +lucasSpecialCase1 = assertEqual "lucas" (lucas 1) 1 + +-- | Check that 'generalLucas' matches its definition. +generalLucasProperty1 :: AnySign Integer -> AnySign Integer -> NonNegative Int -> Bool +generalLucasProperty1 (AnySign p) (AnySign q) (NonNegative n) = un1 == un1' && vn1 == vn1' && un2 == p * un1 - q * un && vn2 == p * vn1 - q * vn + where + (un, un1, vn, vn1) = generalLucas p q n + (un1', un2, vn1', vn2) = generalLucas p q (n + 1) + +-- | Check that 'generalLucas' 1 (-1) is 'fibonacciPair' plus 'lucasPair'. +generalLucasProperty2 :: NonNegative Int -> Bool +generalLucasProperty2 (NonNegative n) = (un, un1) == fibonacciPair n && (vn, vn1) == lucasPair n + where + (un, un1, vn, vn1) = generalLucas 1 (-1) n + +-- | Check that 'generalLucas' p _ 0 is (0, 1, 2, p). +generalLucasProperty3 :: AnySign Integer -> AnySign Integer -> Bool +generalLucasProperty3 (AnySign p) (AnySign q) = generalLucas p q 0 == (0, 1, 2, p) + +testSuite :: TestTree +testSuite = testGroup "Linear" + [ testGroup "fibonacci" + [ testSmallAndQuick "matches definition" fibonacciProperty1 + , testSmallAndQuick "negative indices" fibonacciProperty2 + , testSmallAndQuick "pair" fibonacciPairProperty + , testCase "fibonacci 0" fibonacciSpecialCase0 + , testCase "fibonacci 1" fibonacciSpecialCase1 + ] + , testGroup "lucas" + [ testSmallAndQuick "matches definition" lucasProperty1 + , testSmallAndQuick "negative indices" lucasProperty2 + , testSmallAndQuick "pair" lucasPairProperty + , testCase "lucas 0" lucasSpecialCase0 + , testCase "lucas 1" lucasSpecialCase1 + ] + , testGroup "generalLucas" + [ testSmallAndQuick "matches definition" generalLucasProperty1 + , testSmallAndQuick "generalLucas 1 (-1)" generalLucasProperty2 + , testSmallAndQuick "generalLucas _ _ 0" generalLucasProperty3 + ] + ] 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 @@ -92,24 +92,14 @@ instance (f (g x)) => (f `Compose` g) x type family ConcatMap (w :: * -> Constraint) (cs :: [*]) :: Constraint -#if __GLASGOW_HASKELL__ >= 708 where ConcatMap w '[] = () ConcatMap w (c ': cs) = (w c, ConcatMap w cs) -#else -type instance ConcatMap w '[] = () -type instance ConcatMap w (c ': cs) = (w c, ConcatMap w cs) -#endif type family Matrix (as :: [* -> Constraint]) (w :: * -> *) (bs :: [*]) :: Constraint -#if __GLASGOW_HASKELL__ >= 708 where Matrix '[] w bs = () Matrix (a ': as) w bs = (ConcatMap (a `Compose` w) bs, Matrix as w bs) -#else -type instance Matrix '[] w bs = () -type instance Matrix (a ': as) w bs = (ConcatMap (a `Compose` w) bs, Matrix as w bs) -#endif type TestableIntegral wrapper = ( Matrix '[Arbitrary, Show, Serial IO] wrapper '[Int, Word, Integer] diff --git a/test-suite/Math/NumberTheory/TestUtils/Compose.hs b/test-suite/Math/NumberTheory/TestUtils/Compose.hs --- a/test-suite/Math/NumberTheory/TestUtils/Compose.hs +++ b/test-suite/Math/NumberTheory/TestUtils/Compose.hs @@ -9,7 +9,6 @@ -- Utils to test Math.NumberTheory -- -{-# LANGUAGE CPP #-} {-# LANGUAGE DeriveGeneric #-} {-# LANGUAGE FlexibleContexts #-} {-# LANGUAGE FlexibleInstances #-} @@ -22,11 +21,8 @@ module Math.NumberTheory.TestUtils.Compose where +import Data.Functor.Classes 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) @@ -35,13 +31,9 @@ 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 (Ord1 f, Ord1 g, Ord a, Real (f (g a))) => Real (Compose f g a) +deriving instance (Ord1 f, Ord1 g, Ord 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 --- a/test-suite/Math/NumberTheory/TestUtils/Wrappers.hs +++ b/test-suite/Math/NumberTheory/TestUtils/Wrappers.hs @@ -47,25 +47,13 @@ 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 @@ -81,25 +69,13 @@ 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 @@ -115,25 +91,13 @@ 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 @@ -148,25 +112,13 @@ 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 @@ -182,25 +134,13 @@ 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 @@ -216,25 +156,13 @@ 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 diff --git a/test-suite/Math/NumberTheory/ZetaTests.hs b/test-suite/Math/NumberTheory/ZetaTests.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/NumberTheory/ZetaTests.hs @@ -0,0 +1,116 @@ +-- | +-- Module: Math.NumberTheory.ZetaTests +-- Copyright: (c) 2016 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- Stability: Provisional +-- +-- Tests for Math.NumberTheory.Zeta +-- + +{-# OPTIONS_GHC -fno-warn-type-defaults #-} + +module Math.NumberTheory.ZetaTests + ( testSuite + ) where + +import Test.Tasty +import Test.Tasty.HUnit + +import Math.NumberTheory.Zeta +import Math.NumberTheory.TestUtils + +assertEqualUpToEps :: String -> Double -> Double -> Double -> Assertion +assertEqualUpToEps msg eps expected actual + = assertBool msg (abs (expected - actual) < eps) + +epsilon :: Double +epsilon = 1e-14 + +zetasEvenSpecialCase1 :: Assertion +zetasEvenSpecialCase1 + = assertEqual "zeta(0) = -1/2" + (approximateValue $ zetasEven !! 0) + (-1 / 2) + +zetasEvenSpecialCase2 :: Assertion +zetasEvenSpecialCase2 + = assertEqualUpToEps "zeta(2) = pi^2/6" epsilon + (approximateValue $ zetasEven !! 1) + (pi * pi / 6) + +zetasEvenSpecialCase3 :: Assertion +zetasEvenSpecialCase3 + = assertEqualUpToEps "zeta(4) = pi^4/90" epsilon + (approximateValue $ zetasEven !! 2) + (pi ^ 4 / 90) + +zetasEvenProperty1 :: Positive Int -> Bool +zetasEvenProperty1 (Positive m) + = zetaM < 1 + || zetaM > zetaM1 + where + zetaM = approximateValue (zetasEven !! m) + zetaM1 = approximateValue (zetasEven !! (m + 1)) + +zetasEvenProperty2 :: Positive Int -> Bool +zetasEvenProperty2 (Positive m) + = abs (zetaM - zetaM') < epsilon + where + zetaM = approximateValue (zetasEven !! m) + zetaM' = zetas' !! (2 * m) + +zetas' :: [Double] +zetas' = zetas epsilon + +zetasSpecialCase1 :: Assertion +zetasSpecialCase1 + = assertEqual "zeta(1) = Infinity" + (zetas' !! 1) + (1 / 0) + +zetasSpecialCase2 :: Assertion +zetasSpecialCase2 + = assertEqualUpToEps "zeta(3) = 1.2020569" epsilon + (zetas' !! 3) + 1.2020569031595942853997381615114499908 + +zetasSpecialCase3 :: Assertion +zetasSpecialCase3 + = assertEqualUpToEps "zeta(5) = 1.0369277" epsilon + (zetas' !! 5) + 1.0369277551433699263313654864570341681 + +zetasProperty1 :: Positive Int -> Bool +zetasProperty1 (Positive m) + = zetaM >= zetaM1 + && zetaM1 >= 1 + where + zetaM = zetas' !! m + zetaM1 = zetas' !! (m + 1) + +zetasProperty2 :: NonNegative Int -> NonNegative Int -> Bool +zetasProperty2 (NonNegative e1) (NonNegative e2) + = maximum (take 25 $ drop 2 $ zipWith ((abs .) . (-)) (zetas eps1) (zetas eps2)) < eps1 + eps2 + where + eps1, eps2 :: Double + eps1 = 1.0 / 2 ^ e1 + eps2 = 1.0 / 2 ^ e2 + +testSuite :: TestTree +testSuite = testGroup "Zeta" + [ testGroup "zetasEven" + [ testCase "zeta(0)" zetasEvenSpecialCase1 + , testCase "zeta(2)" zetasEvenSpecialCase2 + , testCase "zeta(4)" zetasEvenSpecialCase3 + , testSmallAndQuick "zeta(2n) > zeta(2n+2)" zetasEvenProperty1 + , testSmallAndQuick "zetasEven matches zetas" zetasEvenProperty2 + ] + , testGroup "zetas" + [ testCase "zeta(1)" zetasSpecialCase1 + , testCase "zeta(3)" zetasSpecialCase2 + , testCase "zeta(5)" zetasSpecialCase3 + , testSmallAndQuick "zeta(n) > zeta(n+1)" zetasProperty1 + , testSmallAndQuick "precision" zetasProperty2 + ] + ] diff --git a/test-suite/Test.hs b/test-suite/Test.hs --- a/test-suite/Test.hs +++ b/test-suite/Test.hs @@ -3,9 +3,8 @@ import qualified Math.NumberTheory.GCDTests as GCD import qualified Math.NumberTheory.GCD.LowLevelTests as GCDLowLevel -import qualified Math.NumberTheory.LogarithmsTests as Logarithms - -import qualified Math.NumberTheory.LucasTests as Lucas +import qualified Math.NumberTheory.Recurrencies.BilinearTests as RecurrenciesBilinear +import qualified Math.NumberTheory.Recurrencies.LinearTests as RecurrenciesLinear import qualified Math.NumberTheory.ModuliTests as Moduli @@ -15,18 +14,20 @@ import qualified Math.NumberTheory.Powers.CubesTests as Cubes import qualified Math.NumberTheory.Powers.FourthTests as Fourth import qualified Math.NumberTheory.Powers.GeneralTests as General -import qualified Math.NumberTheory.Powers.IntegerTests as Integer import qualified Math.NumberTheory.Powers.SquaresTests as Squares import qualified Math.NumberTheory.PrimesTests as Primes import qualified Math.NumberTheory.Primes.CountingTests as Counting +import qualified Math.NumberTheory.Primes.FactorisationTests as Factorisation import qualified Math.NumberTheory.Primes.HeapTests as Heap import qualified Math.NumberTheory.Primes.SieveTests as Sieve +import qualified Math.NumberTheory.Primes.TestingTests as Testing import qualified Math.NumberTheory.GaussianIntegersTests as Gaussian import qualified Math.NumberTheory.ArithmeticFunctionsTests as ArithmeticFunctions import qualified Math.NumberTheory.UniqueFactorisationTests as UniqueFactorisation +import qualified Math.NumberTheory.ZetaTests as Zeta main :: IO () main = defaultMain tests @@ -37,18 +38,15 @@ [ Cubes.testSuite , Fourth.testSuite , General.testSuite - , Integer.testSuite , Squares.testSuite ] , testGroup "GCD" [ GCD.testSuite , GCDLowLevel.testSuite ] - , testGroup "Logarithms" - [ Logarithms.testSuite - ] - , testGroup "Lucas" - [ Lucas.testSuite + , testGroup "Recurrencies" + [ RecurrenciesLinear.testSuite + , RecurrenciesBilinear.testSuite ] , testGroup "Moduli" [ Moduli.testSuite @@ -60,8 +58,10 @@ , testGroup "Primes" [ Primes.testSuite , Counting.testSuite + , Factorisation.testSuite , Heap.testSuite , Sieve.testSuite + , Testing.testSuite ] , testGroup "Gaussian" [ Gaussian.testSuite @@ -71,5 +71,8 @@ ] , testGroup "UniqueFactorisation" [ UniqueFactorisation.testSuite + ] + , testGroup "Zeta" + [ Zeta.testSuite ] ]