arithmoi 0.4.3.0 → 0.5.0.0

raw patch · 40 files changed

+1204/−1044 lines, 40 filesdep +exact-pidep +integer-logarithmsdep −transformers-compatdep ~arithmoidep ~basedep ~transformersPVP ok

version bump matches the API change (PVP)

Dependencies added: exact-pi, integer-logarithms

Dependencies removed: transformers-compat

Dependency ranges changed: arithmoi, base, transformers

API changes (from Hackage documentation)

Files

• Changes +37/−0
• LICENSE +1/−1
• Math/NumberTheory/ArithmeticFunctions/Standard.hs +1/−9
• Math/NumberTheory/Logarithms.hs +0/−195
• Math/NumberTheory/Lucas.hs +0/−99
• Math/NumberTheory/Powers/Cubes.hs +0/−3
• Math/NumberTheory/Powers/Fourth.hs +0/−3
• Math/NumberTheory/Powers/General.hs +0/−3
• Math/NumberTheory/Powers/Integer.hs +9/−23
• Math/NumberTheory/Powers/Squares/Internal.hs +0/−3
• Math/NumberTheory/Primes/Factorisation.hs +0/−122
• Math/NumberTheory/Primes/Factorisation/Utils.hs +0/−89
• Math/NumberTheory/Primes/Sieve/Misc.hs +23/−2
• Math/NumberTheory/Primes/Testing/Probabilistic.hs +43/−65
• Math/NumberTheory/Recurrencies/Bilinear.hs +197/−0
• Math/NumberTheory/Recurrencies/Linear.hs +121/−0
• Math/NumberTheory/UniqueFactorisation.hs +1/−9
• Math/NumberTheory/Utils.hs +0/−9
• Math/NumberTheory/Zeta.hs +115/−0
• TODO +0/−7
• arithmoi.cabal +26/−18
• benchmark/Bench.hs +4/−2
• benchmark/Math/NumberTheory/ArithmeticFunctionsBench.hs +11/−20
• benchmark/Math/NumberTheory/PowersBench.hs +0/−2
• benchmark/Math/NumberTheory/PrimesBench.hs +32/−0
• benchmark/Math/NumberTheory/RecurrenciesBench.hs +30/−0
• test-suite/Math/NumberTheory/LogarithmsTests.hs +0/−112
• test-suite/Math/NumberTheory/LucasTests.hs +0/−104
• test-suite/Math/NumberTheory/Powers/CubesTests.hs +2/−0
• test-suite/Math/NumberTheory/Powers/FourthTests.hs +2/−0
• test-suite/Math/NumberTheory/Powers/IntegerTests.hs +0/−41
• test-suite/Math/NumberTheory/Primes/FactorisationTests.hs +44/−0
• test-suite/Math/NumberTheory/Primes/TestingTests.hs +73/−0
• test-suite/Math/NumberTheory/Recurrencies/BilinearTests.hs +196/−0
• test-suite/Math/NumberTheory/Recurrencies/LinearTests.hs +104/−0
• test-suite/Math/NumberTheory/TestUtils.hs +0/−10
• test-suite/Math/NumberTheory/TestUtils/Compose.hs +3/−11
• test-suite/Math/NumberTheory/TestUtils/Wrappers.hs +0/−72
• test-suite/Math/NumberTheory/ZetaTests.hs +116/−0
• test-suite/Test.hs +13/−10

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

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

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

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

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

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

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

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

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

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

@@ -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 @&#8484;/(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 @&#8484;/(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

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

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

@@ -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 <daniel.is.fischer@googlemail.com> -- 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 <https://miller-rabin.appspot.com Deterministic variants of the Miller-Rabin primality test>+--   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 --   <http://mathworld.wolfram.com/Baillie-PSWPrimalityTest.html> 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

@@ -0,0 +1,197 @@+-- |+-- Module:      Math.NumberTheory.Recurrencies.Bilinear+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- 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 <https://en.wikipedia.org/wiki/Stirling_numbers_of_the_first_kind Stirling numbers of the first kind>.+--+-- > > 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 <https://en.wikipedia.org/wiki/Stirling_numbers_of_the_second_kind Stirling numbers of the second kind>.+--+-- > > 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 <https://en.wikipedia.org/wiki/Lah_number Lah numbers>.+-- @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 <https://en.wikipedia.org/wiki/Eulerian_number Eulerian numbers of the first kind>.+--+-- > > 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 <https://en.wikipedia.org/wiki/Eulerian_number#Eulerian_numbers_of_the_second_kind Eulerian numbers of the second kind>.+--+-- > > 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 <https://en.wikipedia.org/wiki/Bernoulli_number Bernoulli numbers>,+-- computed via <https://en.wikipedia.org/wiki/Bernoulli_number#Connection_with_Stirling_numbers_of_the_second_kind connection>+-- 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 #-}

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

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

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

@@ -0,0 +1,115 @@+-- |+-- Module:      Math.NumberTheory.Zeta+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- 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+-- <https://cr.yp.to/bib/2000/borwein.pdf Computational strategies for the Riemann zeta function>+-- 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))

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

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

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

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

@@ -1,5 +1,3 @@-{-# OPTIONS_GHC -fno-warn-warnings-deprecations #-}- module Math.NumberTheory.PowersBench   ( benchSuite   ) where

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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