arithmoi 0.4.1.3 → 0.4.2.0

raw patch · 44 files changed

+2526/−122 lines, 44 filesdep +QuickCheckdep +smallcheckdep +tastydep −hspecdep ~arithmoidep ~basedep ~ghc-primnew-uploaderPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: QuickCheck, smallcheck, tasty, tasty-hunit, tasty-quickcheck, tasty-smallcheck

Dependencies removed: hspec

Dependency ranges changed: arithmoi, base, ghc-prim

API changes (from Hackage documentation)

Files

• Changes +15/−0
• Math/NumberTheory/GCD.hs +15/−7
• Math/NumberTheory/GaussianIntegers.hs +238/−0
• Math/NumberTheory/Logarithms.hs +1/−1
• Math/NumberTheory/Moduli.hs +50/−48
• Math/NumberTheory/MoebiusInversion.hs +5/−3
• Math/NumberTheory/MoebiusInversion/Int.hs +6/−4
• Math/NumberTheory/Powers/Cubes.hs +4/−4
• Math/NumberTheory/Powers/Fourth.hs +8/−4
• Math/NumberTheory/Powers/General.hs +14/−8
• Math/NumberTheory/Powers/Squares.hs +1/−1
• Math/NumberTheory/Primes/Counting.hs +4/−0
• Math/NumberTheory/Primes/Counting/Approximate.hs +23/−10
• Math/NumberTheory/Primes/Counting/Impl.hs +19/−6
• Math/NumberTheory/Primes/Factorisation.hs +16/−1
• Math/NumberTheory/Primes/Factorisation/Montgomery.hs +1/−1
• Math/NumberTheory/Primes/Factorisation/Utils.hs +10/−0
• Math/NumberTheory/Primes/Heap.hs +5/−2
• Math/NumberTheory/Primes/Sieve/Eratosthenes.hs +3/−3
• Math/NumberTheory/Primes/Sieve/Indexing.hs +5/−2
• Math/NumberTheory/Primes/Sieve/Misc.hs +2/−2
• Math/NumberTheory/Unsafe.hs +73/−0
• Math/NumberTheory/Utils.hs +2/−2
• arithmoi.cabal +36/−12
• test-suite/Math/NumberTheory/GCD/LowLevelTests.hs +74/−0
• test-suite/Math/NumberTheory/GCDTests.hs +50/−0
• test-suite/Math/NumberTheory/GaussianIntegersTests.hs +79/−0
• test-suite/Math/NumberTheory/LogarithmsTests.hs +112/−0
• test-suite/Math/NumberTheory/LucasTests.hs +104/−0
• test-suite/Math/NumberTheory/ModuliTests.hs +222/−0
• test-suite/Math/NumberTheory/MoebiusInversion/IntTests.hs +45/−0
• test-suite/Math/NumberTheory/MoebiusInversionTests.hs +45/−0
• test-suite/Math/NumberTheory/Powers/CubesTests.hs +154/−0
• test-suite/Math/NumberTheory/Powers/FourthTests.hs +145/−0
• test-suite/Math/NumberTheory/Powers/GeneralTests.hs +128/−0
• test-suite/Math/NumberTheory/Powers/IntegerTests.hs +41/−0
• test-suite/Math/NumberTheory/Powers/SquaresTests.hs +155/−0
• test-suite/Math/NumberTheory/Primes/CountingTests.hs +148/−0
• test-suite/Math/NumberTheory/Primes/HeapTests.hs +67/−0
• test-suite/Math/NumberTheory/Primes/SieveTests.hs +69/−0
• test-suite/Math/NumberTheory/PrimesTests.hs +40/−0
• test-suite/Math/NumberTheory/TestUtils.hs +226/−0
• test-suite/Spec.hs +0/−1
• test-suite/Test.hs +66/−0

@@ -1,3 +1,18 @@+0.4.2.0:+    This release supports GHC 7.6, 7.8 and 8.0.++    Add new cabal flag check-bounds, which replaces all unsafe array functions with safe ones.++    Add basic functions on Gaussian integers.+    Add Moebius mu-function.++    Forbid non-positive moduli in Math.NumberTheory.Moduli.++    Fix out-of-bounds error in Math.NumberTheory.Primes.Heap, Math.NumberTheory.Primes.Sieve and Math.NumberTheory.MoebiusInversion.+    Fix 32-bit build.+    Fix binaryGCD on negative numbers.+    Fix highestPower (various issues).+ 0.4.1.0:     Add integerLog10 variants at Bas van Dijk's request and expose     Math.NumberTheory.Powers.Integer, with an added integerWordPower.

@@ -83,10 +83,13 @@ binaryGCDImpl a 0 = abs a binaryGCDImpl 0 b = abs b binaryGCDImpl a b =-    case shiftToOddCount a of+    case shiftToOddCount a' of       (!za, !oa) ->-        case shiftToOddCount b of+        case shiftToOddCount b' of           (!zb, !ob) -> gcdOdd (abs oa) (abs ob) `shiftL` min za zb+    where+      a' = abs a+      b' = abs b  {-# SPECIALISE extendedGCD :: Int -> Int -> (Int, Int, Int),                               Word -> Word -> (Word, Word, Word),@@ -95,16 +98,21 @@ -- | Calculate the greatest common divisor of two numbers and coefficients --   for the linear combination. -----   Satisfies:+--   For signed types satisfies: -- -- > case extendedGCD a b of -- >   (d, u, v) -> u*a + v*b == d--- >--- > d == gcd a b+-- >                && d == gcd a b -----   and, for signed types,+--   For unsigned and bounded types the property above holds, but since @u@ and @v@ must also be unsigned,+--   the result may look weird. E. g., on 64-bit architecture ----- >+-- > extendedGCD (2 :: Word) (3 :: Word) == (1, 2^64-1, 1)+--+--   For unsigned and unbounded types (like 'Numeric.Natural.Natural') the result is undefined.+--+--   For signed types we also have+-- -- > abs u < abs b || abs b <= 1 -- > -- > abs v < abs a || abs a <= 1

@@ -0,0 +1,238 @@+-- |+-- Module:      Math.NumberTheory.GaussianIntegers+-- Copyright:   (c) 2016 Chris Fredrickson+-- Licence:     MIT+-- Maintainer:  Chris Fredrickson <chris.p.fredrickson@gmail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- This module exports functions for manipulating Gaussian integers, including+-- computing their prime factorisations.+--++{-# LANGUAGE BangPatterns #-}+module Math.NumberTheory.GaussianIntegers (+    GaussianInteger((:+)),+    ι,+    real,+    imag,+    conjugate,+    norm,+    divModG,+    divG,+    modG,+    quotRemG,+    quotG,+    remG,+    (.^),+    isPrime,+    primes,+    gcdG,+    gcdG',+    findPrime,+    findPrime',+    factorise,+) where++import qualified Math.NumberTheory.Moduli as Moduli+import qualified Math.NumberTheory.Powers as Powers+import qualified Math.NumberTheory.Primes.Factorisation as Factorisation+import qualified Math.NumberTheory.Primes.Sieve as Sieve+import qualified Math.NumberTheory.Primes.Testing as Testing++infix 6 :++infixr 8 .^+-- |A Gaussian integer is a+bi, where a and b are both integers.+data GaussianInteger = (:+) { real :: !Integer, imag :: !Integer } deriving (Eq)++-- |The imaginary unit, where+--+-- > ι .^ 2 == -1+ι :: GaussianInteger+ι = 0 :+ 1++instance Show GaussianInteger where+    show (a :+ b)+        | b == 0     = show a+        | a == 0     = s ++ b'+        | otherwise  = show a ++ op ++ b'+        where+            b' = if abs b == 1 then "ι" else show (abs b) ++ "*ι"+            op = if b > 0 then "+" else "-"+            s  = if b > 0 then "" else "-"++instance Num GaussianInteger where+    (+) (a :+ b) (c :+ d) = (a + c) :+ (b + d)+    (*) (a :+ b) (c :+ d) = (a * c - b * d) :+ (a * d + b * c)+    abs z@(a :+ b)+        | a == 0 && b == 0 =   z             -- origin+        | a >  0 && b >= 0 =   z             -- first quadrant: (0, inf) x [0, inf)i+        | a <= 0 && b >  0 =   b  :+ (-a)    -- second quadrant: (-inf, 0] x (0, inf)i+        | a <  0 && b <= 0 = (-a) :+ (-b)    -- third quadrant: (-inf, 0) x (-inf, 0]i+        | otherwise        = (-b) :+   a     -- fourth quadrant: [0, inf) x (-inf, 0)i+    negate (a :+ b) = (-a) :+ (-b)+    fromInteger n = n :+ 0+    signum z@(a :+ b)+        | a == 0 && b == 0 = z               -- hole at origin+        | otherwise        = z `divG` abs z++-- |Simultaneous 'quot' and 'rem'.+quotRemG :: GaussianInteger -> GaussianInteger -> (GaussianInteger, GaussianInteger)+quotRemG = divHelper quot++-- |Gaussian integer division, truncating toward zero.+quotG :: GaussianInteger -> GaussianInteger -> GaussianInteger+n `quotG` d = q where (q,_) = quotRemG n d++-- |Gaussian integer remainder, satisfying+--+-- > (x `quotG` y)*y + (x `remG` y) == x+remG :: GaussianInteger -> GaussianInteger -> GaussianInteger+n `remG`  d = r where (_,r) = quotRemG n d++-- |Simultaneous 'div' and 'mod'.+divModG :: GaussianInteger -> GaussianInteger -> (GaussianInteger, GaussianInteger)+divModG = divHelper div++-- |Gaussian integer division, truncating toward negative infinity.+divG :: GaussianInteger -> GaussianInteger -> GaussianInteger+n `divG` d = q where (q,_) = divModG n d++-- |Gaussian integer remainder, satisfying+--+-- > (x `divG` y)*y + (x `modG` y) == x+modG :: GaussianInteger -> GaussianInteger -> GaussianInteger+n `modG` d = r where (_,r) = divModG n d++divHelper :: (Integer -> Integer -> Integer) -> GaussianInteger -> GaussianInteger -> (GaussianInteger, GaussianInteger)+divHelper divide g h =+    let nr :+ ni = g * conjugate h+        denom = norm h+        q = divide nr denom :+ divide ni denom+        p = h * q+    in (q, g - p)++-- |Conjugate a Gaussian integer.+conjugate :: GaussianInteger -> GaussianInteger+conjugate (r :+ i) = r :+ (-i)++-- |The square of the magnitude of a Gaussian integer.+norm :: GaussianInteger -> Integer+norm (x :+ y) = x * x + y * y++-- |Compute whether a given Gaussian integer is prime.+isPrime :: GaussianInteger -> Bool+isPrime g@(x :+ y)+    | x == 0 && y /= 0 = abs y `mod` 4 == 3 && Testing.isPrime y+    | y == 0 && x /= 0 = abs x `mod` 4 == 3 && Testing.isPrime x+    | otherwise        = Testing.isPrime $ norm g++-- |An infinite list of the Gaussian primes. Uses primes in Z to exhaustively+-- generate all Gaussian primes, but not quite in order of ascending magnitude.+primes :: [GaussianInteger]+primes = [ g+         | p <- Sieve.primes+         , g <- if p `mod` 4 == 3+                then [p :+ 0]+                else+                    if p == 2+                    then [1 :+ 1]+                    else let x :+ y = findPrime' p+                         in [x :+ y, y :+ x]+         ]++-- |Compute the GCD of two Gaussian integers. Enforces the precondition that each+-- integer must be in the first quadrant (or zero).+gcdG :: GaussianInteger -> GaussianInteger -> GaussianInteger+gcdG g h = gcdG' (abs g) (abs h)++-- |Compute the GCD of two Gauss integers. Does not check the precondition.+gcdG' :: GaussianInteger -> GaussianInteger -> GaussianInteger+gcdG' g h+    | h == 0    = g --done recursing+    | otherwise = gcdG' h (abs (g `modG` h))++-- |Find a Gaussian integer whose norm is the given prime number.+-- Checks the precondition that p is prime and that p `mod` 4 /= 3.+findPrime :: Integer -> GaussianInteger+findPrime p+    | p == 2 || (p `mod` 4 == 1 && Testing.isPrime p) = findPrime' p+    | otherwise = error "p must be prime, and not congruent to 3 (mod 4)"++-- |Find a Gaussian integer whose norm is the given prime number. Does not+-- check the precondition.+findPrime' :: Integer -> GaussianInteger+findPrime' p =+    let (Just c) = Moduli.sqrtModP (-1) p+        k  = Powers.integerSquareRoot p+        bs = [1 .. k]+        asbs = map (\b' -> ((b' * c) `mod` p, b')) bs+        (a, b) = head [ (a', b') | (a', b') <- asbs, a' <= k]+    in a :+ b++-- |Raise a Gaussian integer to a given power.+(.^) :: (Integral a) => GaussianInteger -> a -> GaussianInteger+a .^ e+    | e < 0 && norm a == 1 =+        case a of+            1    :+ 0 -> 1+            (-1) :+ 0 -> if even e then 1 else (-1)+            0    :+ 1 -> (0 :+ (-1)) .^ (abs e `mod` 4)+            _         -> (0 :+ 1) .^ (abs e `mod` 4)+    | e < 0     = error "Cannot exponentiate non-unit Gaussian Int to negative power"+    | a == 0    = 0+    | e == 0    = 1+    | even e    = s * s+    | otherwise = a * a .^ (e - 1)+    where+    s = a .^ div e 2++-- |Compute the prime factorization of a Gaussian integer. This is unique up to units (+/- 1, +/- i).+factorise :: GaussianInteger -> [(GaussianInteger, Int)]+factorise g+    | g == 0    = error "0 has no prime factorisation"+    | g == 1    = []+    | otherwise =+        let helper :: [(Integer, Int)] -> GaussianInteger -> [(GaussianInteger, Int)] -> [(GaussianInteger, Int)]+            helper [] g' fs = (if g' == 1 then [] else [(g', 1)]) ++ fs    -- include the unit, if it isn't 1+            helper ((!p, !e) : pt) g' fs+                | p `mod` 4 == 3 =+                    -- prime factors congruent to 3 mod 4 are simple.+                    let pow = div e 2+                        gp = fromInteger p+                    in helper pt (g' `divG` (gp .^ pow)) ((gp, pow) : fs)+                | otherwise      =+                    -- general case: for every prime factor of the magnitude+                    -- squared, find a Gaussian prime whose magnitude squared+                    -- is that prime. Then find out how many times the original+                    -- number is divisible by that Gaussian prime and its+                    -- conjugate. The order that the prime factors are+                    -- processed doesn't really matter, but it is reversed so+                    -- that the Gaussian primes will be in order of increasing+                    -- magnitude.+                    let gp = findPrime' p+                        (!gNext, !facs) = trialDivide g' [gp, abs $ conjugate gp] []+                    in helper pt gNext (facs ++ fs)+        in helper (reverse . Factorisation.factorise $ norm g) g []++-- Divide a Gaussian integer by a set of (relatively prime) Gaussian integers,+-- as many times as possible, and return the final quotient as well as a count+-- of how many times each factor divided the original.+trialDivide :: GaussianInteger -> [GaussianInteger] -> [(GaussianInteger, Int)] -> (GaussianInteger, [(GaussianInteger, Int)])+trialDivide g [] fs = (g, fs)+trialDivide g (pf : pft) fs+    | g `modG` pf == 0 =+        let (cnt, g') = countEvenDivisions g pf+        in trialDivide g' pft ((pf, cnt) : fs)+    | otherwise    = trialDivide g pft fs++-- Divide a Gaussian integer by a possible factor, and return how many times+-- the factor divided it evenly, as well as the result of dividing the original+-- that many times.+countEvenDivisions :: GaussianInteger -> GaussianInteger -> (Int, GaussianInteger)+countEvenDivisions g pf = helper g 0+    where+    helper :: GaussianInteger -> Int -> (Int, GaussianInteger)+    helper g' acc+        | g' `modG` pf == 0 = helper (g' `divG` pf) (1 + acc)+        | otherwise     = (acc, g')

@@ -36,10 +36,10 @@  import Data.Bits import Data.Array.Unboxed-import Data.Array.Base (unsafeAt)  import Math.NumberTheory.Logarithms.Internal import Math.NumberTheory.Powers.Integer+import Math.NumberTheory.Unsafe #if __GLASGOW_HASKELL__ < 707 import Math.NumberTheory.Utils  (isTrue#) #endif

@@ -34,52 +34,50 @@  #include "MachDeps.h" -#if __GLASGOW_HASKELL__ < 709+#if __GLASGOW_HASKELL__ < 709 || WORD_SIZE_IN_BITS == 32 import Data.Word #endif import Data.Bits import Data.Array.Unboxed-import Data.Array.Base (unsafeAt)-import Data.Maybe (fromJust) import Data.List (nub) import Control.Monad (foldM, liftM2)  import Math.NumberTheory.Utils (shiftToOddCount, splitOff) import Math.NumberTheory.GCD (extendedGCD) import Math.NumberTheory.Primes.Heap (sieveFrom)+import Math.NumberTheory.Unsafe+ -- Guesstimated startup time for the Heap algorithm is lower than -- the cost to sieve an entire chunk. --- | Invert a number relative to a modulus.+-- | Invert a number relative to a positive modulus. --   If @number@ and @modulus@ are coprime, the result is --   @Just inverse@ where ----- >    (number * inverse) `mod` (abs modulus) == 1--- >    0 <= inverse < abs modulus+-- >    (number * inverse) `mod` modulus == 1+-- >    0 <= inverse < modulus -----   unless @modulus == 0@ and @abs number == 1@, in which case the---   result is @Just number@.---   If @gcd number modulus > 1@, the result is @Nothing@.+--   If @number `mod` modulus == 0@ or @gcd number modulus > 1@, the result is @Nothing@. invertMod :: Integer -> Integer -> Maybe Integer-invertMod k 0 = if k == 1 || k == (-1) then Just k else Nothing-invertMod k m = wrap $ go False 1 0 m' k'+invertMod k m+  | m <= 0 = error "Math.NumberTheory.Moduli.invertMod: non-positive modulus"+  | otherwise = wrap $ go False 1 0 m k'   where-    m' = abs m-    k' | r < 0     = r+m'+    k' | r < 0     = r+m        | otherwise = r          where-           r = k `rem` m'-    wrap x = case (x*k') `rem` m' of+           r = k `rem` m+    wrap x = case (x*k') `rem` m of                1 -> Just x                _ -> Nothing-    -- Calculate modular inverse of k' modulo m' by continued fraction expansion-    -- of m'/k', say [a_0,a_1,...,a_s]. Let the convergents be p_j/q_j.+    -- Calculate modular inverse of k' modulo m by continued fraction expansion+    -- of m/k', say [a_0,a_1,...,a_s]. Let the convergents be p_j/q_j.     -- Starting from j = -2, the arguments of go are-    -- (p_j/q_j) > m'/k', p_{j+1}, p_j, and n, d with n/d = [a_{j+2},...,a_s].-    -- Since m'/k' = p_s/q_s, and p_j*q_{j+1} - p_{j+1}*q_j = (-1)^(j+1), we have-    -- p_{s-1}*k' - q_{s-1}*m' = (-1)^s * gcd m' k', so if the inverse exists,+    -- (p_j/q_j) > m/k', p_{j+1}, p_j, and n, d with n/d = [a_{j+2},...,a_s].+    -- Since m/k' = p_s/q_s, and p_j*q_{j+1} - p_{j+1}*q_j = (-1)^(j+1), we have+    -- p_{s-1}*k' - q_{s-1}*m = (-1)^s * gcd m k', so if the inverse exists,     -- it is either p_{s-1} or -p_{s-1}, depending on whether s is even or odd.-    go !b _ po _ 0 = if b then po else (m'-po)+    go !b _ po _ 0 = if b then po else (m-po)     go b !pn po n d = case n `quotRem` d of                         (q,r) -> go (not b) (q*pn+po) pn d r @@ -161,7 +159,7 @@ -- -- > powerMod base exponent modulus -----   calculates @(base ^ exponent) \`mod\` modulus@ by repeated squaring and reduction.+--   calculates @(base ^ exponent) \`mod\` modulus@ by repeated squaring and reduction. Modulus must be positive. --   If @exponent < 0@ and @base@ is invertible modulo @modulus@, @(inverse ^ |exponent|) \`mod\` modulus@ --   is calculated. This function does some input checking and sanitation before calling the unsafe worker. {-# RULES@@ -176,18 +174,17 @@   #-} powerModImpl :: (Integral a, Bits a) => Integer -> a -> Integer -> Integer powerModImpl base expo md-  | md == 0     = base ^ expo-  | md' == 1    = 0+  | md <= 0     = error "Math.NumberTheory.Moduli.powerMod: non-positive modulus"+  | md == 1     = 0   | expo == 0   = 1   | bse' == 1   = 1-  | expo < 0    = case invertMod bse' md' of-                    Just i  -> powerMod'Impl i (negate expo) md'+  | expo < 0    = case invertMod bse' md of+                    Just i  -> powerMod'Impl i (negate expo) md                     Nothing -> error "Math.NumberTheory.Moduli.powerMod: Base isn't invertible with respect to modulus"   | bse' == 0   = 0-  | otherwise   = powerMod'Impl bse' expo md'+  | otherwise   = powerMod'Impl bse' expo md     where-      md' = abs md-      bse' = if base < 0 || md' <= base then base `mod` md' else base+      bse' = if base < 0 || md <= base then base `mod` md else base  -- | Modular power worker without input checking. --   Assumes all arguments strictly positive and modulus greater than 1.@@ -213,21 +210,20 @@ -- | Specialised version of 'powerMod' for 'Integer' exponents. --   Reduces the number of shifts of the exponent since shifting --   large 'Integer's is expensive. Call this function directly---   if you don't want or can't rely on rewrite rules.+--   if you don't want or can't rely on rewrite rules. Modulus must be positive. powerModInteger :: Integer -> Integer -> Integer -> Integer powerModInteger base ex mdl-  | mdl == 0    = base ^ ex-  | mdl' == 1   = 0+  | mdl <= 0     = error "Math.NumberTheory.Moduli.powerModInteger: non-positive modulus"+  | mdl == 1    = 0   | ex == 0     = 1-  | ex < 0      = case invertMod bse' mdl' of-                    Just i  -> powerModInteger' i (negate ex) mdl'+  | ex < 0      = case invertMod bse' mdl of+                    Just i  -> powerModInteger' i (negate ex) mdl                     Nothing -> error "Math.NumberTheory.Moduli.powerMod: Base isn't invertible with respect to modulus"   | bse' == 0   = 0   | bse' == 1   = 1-  | otherwise   = powerModInteger' bse' ex mdl'+  | otherwise   = powerModInteger' bse' ex mdl     where-      mdl' = abs mdl-      bse' = if base < 0 || mdl' <= base then base `mod` mdl' else base+      bse' = if base < 0 || mdl <= base then base `mod` mdl else base  -- | Specialised worker without input checks. Makes the same assumptions --   as the general version 'powerMod''.@@ -322,7 +318,7 @@                         _      -> []  -- | @sqrtModP' square prime@ finds a square root of @square@ modulo---   prime. @prime@ /must/ be a (positive) prime, and @sqaure@ /must/ be a+--   prime. @prime@ /must/ be a (positive) prime, and @square@ /must/ be a positive --   quadratic residue modulo @prime@, i.e. @'jacobi square prime == 1@. --   The precondition is /not/ checked. sqrtModP' :: Integer -> Integer -> Integer@@ -333,7 +329,7 @@  -- | @tonelliShanks square prime@ calculates a square root of @square@ --   modulo @prime@, where @prime@ is a prime of the form @4*k + 1@ and---   @square@ is a quadratic residue modulo @prime@, using the+--   @square@ is a positive quadratic residue modulo @prime@, using the --   Tonelli-Shanks algorithm. --   No checks on the input are performed. tonelliShanks :: Integer -> Integer -> Integer@@ -366,15 +362,15 @@ sqrtModPP :: Integer -> (Integer,Int) -> Maybe Integer sqrtModPP n (2,e) = sqM2P n e sqrtModPP n (prime,expo) = case sqrtModP n prime of-                             Just r -> Just $ fixup r+                             Just r -> fixup r                              _      -> Nothing   where     fixup r = let diff' = r*r-n               in if diff' == 0-                   then r+                   then Just r                    else case splitOff prime diff' of-                          (e,q) | expo <= e -> r-                                | otherwise -> hoist (fromJust $ invertMod (2*r) prime) r (q `mod` prime) (prime^e)+                          (e,q) | expo <= e -> Just r+                                | otherwise -> fmap (\inv -> hoist inv r (q `mod` prime) (prime^e)) (invertMod (2*r) prime)                       --     hoist inv root elim pp         | diff' == 0    = root'@@ -418,13 +414,17 @@ -- | @sqrtModF n primePowers@ calculates a square root of @n@ modulo --   @product [p^k | (p,k) <- primePowers]@ if one exists and all primes --   are distinct.+--   The list must be non-empty, @n@ must be coprime with all primes. sqrtModF :: Integer -> [(Integer,Int)] -> Maybe Integer+sqrtModF _ []  = Nothing sqrtModF n pps = do roots <- mapM (sqrtModPP n) pps                     chineseRemainder $ zip roots (map (uncurry (^)) pps)  -- | @sqrtModFList n primePowers@ calculates all square roots of @n@ modulo --   @product [p^k | (p,k) <- primePowers]@ if all primes are distinct.+--   The list must be non-empty, @n@ must be coprime with all primes. sqrtModFList :: Integer -> [(Integer,Int)] -> [Integer]+sqrtModFList _ []  = [] sqrtModFList n pps = map fst $ foldl1 (liftM2 comb) cs   where     ms :: [Integer]@@ -439,6 +439,7 @@ --   square roots of @n@ modulo @prime^expo@. The same restriction --   as in 'sqrtModPP' applies to the arguments. sqrtModPPList :: Integer -> (Integer,Int) -> [Integer]+sqrtModPPList n (2,1) = [n `mod` 2] sqrtModPPList n (2,expo)     = case sqM2P n expo of         Just r -> let m = 1 `shiftL` (expo-1)@@ -458,13 +459,14 @@ -- > r ≡ r_k (mod m_k) -- > -----   if all moduli are pairwise coprime. If not all moduli are---   pairwise coprime, the result is @Nothing@ regardless of whether+--   if all moduli are positive and pairwise coprime. Otherwise+--   the result is @Nothing@ regardless of whether --   a solution exists. chineseRemainder :: [(Integer,Integer)] -> Maybe Integer chineseRemainder remainders = foldM addRem 0 remainders   where     !modulus = product (map snd remainders)+    addRem acc (_,1) = Just acc     addRem acc (r,m) = do         let cf = modulus `quot` m         inv <- invertMod cf m@@ -489,21 +491,21 @@                         Int -> Bool,                         Word -> Bool   #-}-evenI :: (Integral a, Bits a) => a -> Bool+evenI :: Integral a => a -> Bool evenI n = fromIntegral n .&. 1 == (0 :: Int)  {-# SPECIALISE rem4 :: Integer -> Int,                        Int -> Int,                        Word -> Int   #-}-rem4 :: (Integral a, Bits a) => a -> Int+rem4 :: Integral a => a -> Int rem4 n = fromIntegral n .&. 3  {-# SPECIALISE rem8 :: Integer -> Int,                        Int -> Int,                        Word -> Int   #-}-rem8 :: (Integral a, Bits a) => a -> Int+rem8 :: Integral a => a -> Int rem8 n = fromIntegral n .&. 7  jac2 :: UArray Int Int

@@ -15,10 +15,11 @@     ) where  import Data.Array.ST-import Data.Array.Base+import Control.Monad import Control.Monad.ST  import Math.NumberTheory.Powers.Squares+import Math.NumberTheory.Unsafe  -- | @totientSum n@ is, for @n > 0@, the sum of @[totient k | k <- [1 .. n]]@, --   computed via generalised Moebius inversion.@@ -69,7 +70,7 @@ --   Since the function arguments are used as array indices, the domain of --   @f@ is restricted to 'Int'. -----   The value @f n@ is then computed by @generalInversion g n@). Note that when+--   The value @f n@ is then computed by @generalInversion g n@. Note that when --   many values of @f@ are needed, there are far more efficient methods, this --   method is only appropriate to compute isolated values of @f@. generalInversion :: (Int -> Integer) -> Int -> Integer@@ -90,7 +91,8 @@         small <- newArray_ (0,mk0) :: ST s (STArray s Int Integer)         unsafeWrite small 0 0         unsafeWrite small 1 $! (fun 1)-        unsafeWrite small 2 $! (fun 2 - fun 1)+        when (mk0 >= 2) $+            unsafeWrite small 2 $! (fun 2 - fun 1)         let calcit switch change i                 | mk0 < i   = return (switch,change)                 | i == change = calcit (switch+1) (change + 4*switch+6) i

@@ -16,11 +16,11 @@     ) where  import Data.Array.ST-import Data.Array.Base+import Control.Monad import Control.Monad.ST  import Math.NumberTheory.Powers.Squares-+import Math.NumberTheory.Unsafe  -- | @totientSum n@ is, for @n > 0@, the sum of @[totient k | k <- [1 .. n]]@, --   computed via generalised Moebius inversion.@@ -71,11 +71,12 @@ --   That bears the risk of overflow, however, so be sure to use it only when it's --   safe. -----   The value @f n@ is then computed by @generalInversion g n@). Note that when+--   The value @f n@ is then computed by @generalInversion g n@. Note that when --   many values of @f@ are needed, there are far more efficient methods, this --   method is only appropriate to compute isolated values of @f@. generalInversion :: (Int -> Int) -> Int -> Int generalInversion fun n+    | n < 1     = error "Moebius inversion only defined on positive domain"     | n == 1    = fun 1     | n == 2    = fun 2 - fun 1     | n == 3    = fun 3 - 2*fun 1@@ -91,7 +92,8 @@         small <- newArray_ (0,mk0) :: ST s (STUArray s Int Int)         unsafeWrite small 0 0         unsafeWrite small 1 (fun 1)-        unsafeWrite small 2 (fun 2 - fun 1)+        when (mk0 >= 2) $+            unsafeWrite small 2 (fun 2 - fun 1)         let calcit switch change i                 | mk0 < i   = return (switch,change)                 | i == change = calcit (switch+1) (change + 4*switch+6) i

@@ -21,7 +21,6 @@ #include "MachDeps.h"  import Data.Array.Unboxed-import Data.Array.Base import Data.Array.ST  import Data.Bits@@ -34,6 +33,7 @@ import GHC.Integer.GMP.Internals  import Math.NumberTheory.Logarithms.Internal (integerLog2#)+import Math.NumberTheory.Unsafe #if __GLASGOW_HASKELL__ < 707 import Math.NumberTheory.Utils (isTrue#) #endif@@ -61,9 +61,9 @@ --   that is, the largest integer @r@ such that @r^3 <= n@. --   The precondition @n >= 0@ is not checked. {-# RULES-"integerCubeRoot'/Int"  integerCubeRoot' = cubeRootInt'-"integerCubeRoot'/Word" integerCubeRoot' = cubeRootWord-"integerCubeRoot'/Igr"  integerCubeRoot' = cubeRootIgr+"integerCubeRoot'/Int"     integerCubeRoot' = cubeRootInt'+"integerCubeRoot'/Word"    integerCubeRoot' = cubeRootWord+"integerCubeRoot'/Integer" integerCubeRoot' = cubeRootIgr   #-} {-# INLINE [1] integerCubeRoot' #-} integerCubeRoot' :: Integral a => a -> a

@@ -26,7 +26,6 @@  import Data.Array.Unboxed import Data.Array.ST-import Data.Array.Base (unsafeAt, unsafeWrite)  import Data.Bits #if __GLASGOW_HASKELL__ < 705@@ -34,6 +33,7 @@ #endif  import Math.NumberTheory.Logarithms.Internal (integerLog2#)+import Math.NumberTheory.Unsafe #if __GLASGOW_HASKELL__ < 707 import Math.NumberTheory.Utils (isTrue#) #endif@@ -41,6 +41,10 @@ -- | Calculate the integer fourth root of a nonnegative number, --   that is, the largest integer @r@ with @r^4 <= n@. --   Throws an error on negaitve input.+{-# SPECIALISE integerFourthRoot :: Int -> Int,+                                    Integer -> Integer,+                                    Word -> Word+  #-} integerFourthRoot :: Integral a => a -> a integerFourthRoot n     | n < 0     = error "integerFourthRoot: negative argument"@@ -50,9 +54,9 @@ --   that is, the largest integer @r@ with @r^4 <= n@. --   The condition is /not/ checked. {-# RULES-"integerFourthRoot'/Int"  integerFourthRoot' = biSqrtInt-"integerFourthRoot'/Word" integerFourthRoot' = biSqrtWord-"integerFourthRoot'/Igr"  integerFourthRoot' = biSqrtIgr+"integerFourthRoot'/Int"     integerFourthRoot' = biSqrtInt+"integerFourthRoot'/Word"    integerFourthRoot' = biSqrtWord+"integerFourthRoot'/Integer" integerFourthRoot' = biSqrtIgr   #-} {-# INLINE [1] integerFourthRoot' #-} integerFourthRoot' :: Integral a => a -> a

@@ -162,14 +162,17 @@ --   of the prime exponents if some have been found, otherwise by trying --   prime exponents recursively. highestPower :: Integral a => a -> (a, Int)-highestPower n-  | abs n <= 1  = (n,3)-  | n < 0       = case integerHighPower (toInteger $ negate n) of+highestPower n'+  | abs n <= 1  = (n', 3)+  | n < 0       = case integerHighPower (negate n) of                     (r,e) -> case shiftToOddCount e of                                (k, o) -> (negate $ fromInteger (sqr k r), o)-  | otherwise   = case integerHighPower (toInteger n) of+  | otherwise   = case integerHighPower n of                     (r,e) -> (fromInteger r, e)     where+      n :: Integer+      n = toInteger n'+       sqr :: Int -> Integer -> Integer       sqr 0 m = m       sqr k m = sqr (k-1) (m*m)@@ -199,7 +202,8 @@ newtonK :: Integral a => a -> a -> a -> a newtonK k n a = go (step a)   where-    step m = ((k-1)*m + n `quot` (m^(k-1))) `quot` k+    -- Beware integer overflow in m^(k-1)+    step m = ((k-1)*m + fromInteger (toInteger n `quot` (toInteger m^(k-1)))) `quot` k     go m       | l < m     = go l       | otherwise = m@@ -321,13 +325,14 @@                           | otherwise -> go (Set.fromDistinctAscList $ takeWhile (<= mx) $ take (k+1) (iterate (*2) 1)) o iops   where     mset k st = fst (Set.split (mx+1) (Set.mapMonotonic (*k) st))+    -- unP p m = (k, m / p ^ k), where k is as large as possible such that p ^ k still divides m     unP :: Int -> Int -> (Int,Int)     unP p m = goP 0 m       where         goP :: Int -> Int -> (Int,Int)-        goP !i j = case m `quotRem` p of+        goP !i j = case j `quotRem` p of                      (q,r) | r == 0 -> goP (i+1) q-                           | otherwise -> (0,j)+                           | otherwise -> (i,j)     iops :: [Int]     iops = 3:5:prs     prs :: [Int]@@ -343,5 +348,6 @@       | otherwise =         case unP p m of           (0,_) -> go st m ps-          (k,r) -> go (Set.unions (st:take k (iterate (mset p) st))) r ps+          -- iterate f x = [x, f x, f (f x)...]+          (k,r) -> go (Set.unions (take (k + 1) (iterate (mset p) st))) r ps     go st m [] = go st m [m+1]

@@ -29,7 +29,6 @@  import Data.Array.Unboxed import Data.Array.ST-import Data.Array.Base (unsafeAt, unsafeWrite)  import Data.Bits #if __GLASGOW_HASKELL__ < 705@@ -37,6 +36,7 @@ #endif  import Math.NumberTheory.Logarithms.Internal (integerLog2#)+import Math.NumberTheory.Unsafe #if __GLASGOW_HASKELL__ < 707 import Math.NumberTheory.Utils (isTrue#) #endif

@@ -11,10 +11,14 @@ module Math.NumberTheory.Primes.Counting     ( -- * Exact functions       primeCount+    , primeCountMaxArg     , nthPrime+    , nthPrimeMaxArg       -- * Approximations     , approxPrimeCount+    , approxPrimeCountOverestimateLimit     , nthPrimeApprox+    , nthPrimeApproxUnderestimateLimit     ) where  import Math.NumberTheory.Primes.Counting.Impl

@@ -12,24 +12,37 @@ {-# OPTIONS_HADDOCK hide #-} module Math.NumberTheory.Primes.Counting.Approximate     ( approxPrimeCount+    , approxPrimeCountOverestimateLimit     , nthPrimeApprox+    , nthPrimeApproxUnderestimateLimit     ) where --- | @'approxPrimeCount' n@ gives (for @n > 0@) an+-- For prime p = 3742914359 we have+--   approxPrimeCount p = 178317879+--         primeCount p = 178317880++-- | Following property holds:+--+-- > approxPrimeCount n >= primeCount n || n >= approxPrimeCountOverestimateLimit+approxPrimeCountOverestimateLimit :: Integral a => a+approxPrimeCountOverestimateLimit = 3742914359++-- | @'approxPrimeCount' n@ gives an --   approximation of the number of primes not exceeding --   @n@. The approximation is fairly good for @n@ large enough.---   The number of primes should be slightly overestimated---   (so it is suitable for allocation of storage) and is---   never underestimated for @n <= 10^12@. approxPrimeCount :: Integral a => a -> a-approxPrimeCount = truncate . appi . fromIntegral+approxPrimeCount = truncate . max 0 . appi . fromIntegral --- | @'nthPrimeApprox' n@ gives (for @n > 0@) an+-- | Following property holds:+--+-- > nthPrimeApprox n <= nthPrime n || n >= nthPrimeApproxUnderestimateLimit+nthPrimeApproxUnderestimateLimit :: Integer+nthPrimeApproxUnderestimateLimit = 1000000000000++-- | @'nthPrimeApprox' n@ gives an --   approximation to the n-th prime. The approximation---   is fairly good for @n@ large enough. Dual to---   @'approxPrimeCount'@, this estimate should err---   on the low side (and does for @n < 10^12@).-nthPrimeApprox :: Integral a => a -> a+--   is fairly good for @n@ large enough.+nthPrimeApprox :: Integer -> Integer nthPrimeApprox = max 1 . truncate . nthApp . fromIntegral . max 3  -- Basically the approximation of the prime count by Li(x),

@@ -13,7 +13,12 @@ {-# OPTIONS_GHC -fspec-constr-count=24 #-} #endif {-# OPTIONS_HADDOCK hide #-}-module Math.NumberTheory.Primes.Counting.Impl (primeCount, nthPrime) where+module Math.NumberTheory.Primes.Counting.Impl+    ( primeCount+    , primeCountMaxArg+    , nthPrime+    , nthPrimeMaxArg+    ) where  #include "MachDeps.h" @@ -23,8 +28,8 @@ import Math.NumberTheory.Powers.Squares import Math.NumberTheory.Powers.Cubes import Math.NumberTheory.Logarithms+import Math.NumberTheory.Unsafe -import Data.Array.Base import Data.Array.ST #if !MIN_VERSION_array(0,5,0)     hiding (unsafeThaw)@@ -39,10 +44,14 @@ #define COUNT_T Int #endif +-- | Maximal allowed argument of 'primeCount'. Currently 8e18.+primeCountMaxArg :: Integer+primeCountMaxArg = 8000000000000000000+ -- | @'primeCount' n == &#960;(n)@ is the number of (positive) primes not exceeding @n@. -- --   For efficiency, the calculations are done on 64-bit signed integers, therefore @n@ must---   not exceed @8 * 10^18@.+--   not exceed 'primeCountMaxArg'. -- --   Requires @/O/(n^0.5)@ space, the time complexity is roughly @/O/(n^0.7)@. --   For small bounds, @'primeCount' n@ simply counts the primes not exceeding @n@,@@ -51,7 +60,7 @@ --   <http://en.wikipedia.org/wiki/Prime_counting_function#Algorithms_for_evaluating_.CF.80.28x.29>. primeCount :: Integer -> Integer primeCount n-    | n > 8000000000000000000   = error $ "primeCount: can't handle bound " ++ show n+    | n > primeCountMaxArg = error $ "primeCount: can't handle bound " ++ show n     | n < 2     = 0     | n < 1000  = fromIntegral . length . takeWhile (<= n) . primeList . primeSieve $ max 242 n     | n < 30000 = runST $ do@@ -69,15 +78,19 @@             !pdf = sieveCount ub cs sr         in phn1 - pdf +-- | Maximal allowed argument of 'nthPrime'. Currently 1.5e17.+nthPrimeMaxArg :: Integer+nthPrimeMaxArg = 150000000000000000+ -- | @'nthPrime' n@ calculates the @n@-th prime. Numbering of primes is --   @1@-based, so @'nthPrime' 1 == 2@. -- --   Requires @/O/((n*log n)^0.5)@ space, the time complexity is roughly @/O/((n*log n)^0.7@.---   The argument must be strictly positive, and must not exceed @1.5 * 10^17@.+--   The argument must be strictly positive, and must not exceed 'nthPrimeMaxArg'. nthPrime :: Integer -> Integer nthPrime n     | n < 1         = error "Prime indexing starts at 1"-    | n > 150000000000000000    = error $ "nthPrime: can't handle index " ++ show n+    | n > nthPrimeMaxArg = error $ "nthPrime: can't handle index " ++ show n     | n < 200000    = nthPrimeCt n     | ct0 < n       = tooLow n p0 (n-ct0) approxGap     | otherwise     = tooHigh n p0 (ct0-n) approxGap

@@ -48,6 +48,10 @@     , carmichaelSieve     , sieveCarmichael     , carmichaelFromCanonical+      -- * Moebius function+    , moebius+    , μ+    , moebiusFromCanonical       -- * Divisors     , divisors     , tau@@ -101,7 +105,7 @@ φ = totient  -- | Calculates the Carmichael function for a positive integer, that is,---   the (smallest) exponent of the group of units in @&8484;/(n)@.+--   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"@@ -111,6 +115,17 @@ -- | 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.

@@ -46,7 +46,6 @@ #if __GLASGOW_HASKELL__ < 705 import GHC.Word     -- Moved to GHC.Types #endif-import Data.Array.Base  import System.Random import Control.Monad.State.Strict@@ -63,6 +62,7 @@ import Math.NumberTheory.Primes.Sieve.Eratosthenes import Math.NumberTheory.Primes.Sieve.Indexing import Math.NumberTheory.Primes.Testing.Probabilistic+import Math.NumberTheory.Unsafe import Math.NumberTheory.Utils  -- | @'factorise' n@ produces the prime factorisation of @n@, including

@@ -13,6 +13,7 @@     ( ppTotient     , totientFromCanonical     , carmichaelFromCanonical+    , moebiusFromCanonical     , divisorsFromCanonical     , tauFromCanonical     , divisorSumFromCanonical@@ -49,6 +50,15 @@     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

@@ -22,7 +22,6 @@ module Math.NumberTheory.Primes.Heap (primes, sieveFrom) where  import Data.Array.Unboxed-import Data.Array.Base (unsafeAt, unsafeRead, unsafeWrite) import Data.Array.ST import Control.Monad.ST import Data.List (foldl')@@ -30,6 +29,8 @@ import Data.Word #endif +import Math.NumberTheory.Unsafe+ #ifndef SH_SIZE #define SH_SIZE 31 #endif@@ -331,12 +332,14 @@     | r < m     = down (a-1)     | otherwise = up a       where-        a = min 5758 ((192*r) `quot` 1001 - 1)+        a = max 0 (min 5758 ((192*r) `quot` 1001 - 1))         m = remainders `unsafeAt` a         down k+            | k < 0                         = 0             | r < (remainders `unsafeAt` k) = down (k-1)             | otherwise                     = k         up k+            | k+1 > 5759                        = 5759             | r < (remainders `unsafeAt` (k+1)) = k             | otherwise                         = up (k+1) 

@@ -35,18 +35,18 @@ #include "MachDeps.h"  import Control.Monad.ST-import Data.Array.Base import Data.Array.ST #if !MIN_VERSION_array(0,5,0)                      hiding (unsafeFreeze, unsafeThaw, castSTUArray) #endif import Control.Monad (when) import Data.Bits-#if __GLASGOW_HASKELL__ < 709+#if __GLASGOW_HASKELL__ < 709 || WORD_SIZE_IN_BITS == 32 import Data.Word #endif  import Math.NumberTheory.Powers.Squares (integerSquareRoot)+import Math.NumberTheory.Unsafe import Math.NumberTheory.Utils import Math.NumberTheory.Primes.Counting.Approximate import Math.NumberTheory.Primes.Sieve.Indexing@@ -97,7 +97,7 @@ -- | Compact store of primality flags. data PrimeSieve = PS !Integer {-# UNPACK #-} !(UArray Int Bool) --- | Sieve primes up to (and including) a bound.+-- | Sieve primes up to (and including) a bound (or 7, if bound is smaller). --   For small enough bounds, this is more efficient than --   using the segmented sieve. --

@@ -21,9 +21,10 @@     ) where  import Data.Array.Unboxed-import Data.Array.Base (unsafeAt) import Data.Bits +import Math.NumberTheory.Unsafe+ -- Auxiliary stuff, conversion between number and index, -- remainders modulo 30 and related things. @@ -33,7 +34,9 @@ --   #-} {-# INLINE idxPr #-} idxPr :: Integral a => a -> (Int,Int)-idxPr n0 = (fromIntegral bytes0, rm3)+idxPr n0+    | n0 < 7    = (0, 0)+    | otherwise = (fromIntegral bytes0, rm3)   where     n = if (fromIntegral n0 .&. 1 == (1 :: Int))             then n0 else (n0-1)

@@ -31,7 +31,6 @@     ) where  import Control.Monad.ST-import Data.Array.Base (unsafeRead, unsafeWrite, unsafeAt) import Data.Array.ST import Data.Array.Unboxed import Control.Monad (when)@@ -44,6 +43,7 @@ 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  -- | A compact store of smallest prime factors.@@ -149,7 +149,7 @@                                                  | otherwise -> tdLoop j (integerSquareRoot' j) (ix+1)           where             p = toPrim ix-            pix = unsafeAt sve ix+            pix = unsafeAt sve $ fromIntegral p     curve n = stdGenFactorisation (Just (bound*(bound+2))) (mkStdGen $ fromIntegral n `xor` 0xdecaf00d) Nothing n  -- | @'totientSieve' n@ creates a store of the totients of the numbers not exceeding @n@.

@@ -0,0 +1,73 @@+-- |+-- Module:      Math.NumberTheory.Unsafe+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Layer to switch between safe and unsafe arrays.+--++{-# LANGUAGE CPP #-}++module Math.NumberTheory.Unsafe+  ( UArray+  , bounds+  , castSTUArray+  , unsafeAt+  , unsafeFreeze+  , unsafeNewArray_+  , unsafeRead+  , unsafeThaw+  , unsafeWrite+  ) where++#ifdef CheckBounds++import Data.Array.Base+  ( UArray+  , castSTUArray+  )+import Data.Array.IArray+  ( IArray+  , bounds+  , (!)+  )+import Data.Array.MArray+#if !MIN_VERSION_array(0,5,0)+  hiding (unsafeFreeze, unsafeThaw)+#endif++unsafeAt :: (IArray a e, Ix i) => a i e -> i -> e+unsafeAt = (!)++unsafeFreeze :: (Ix i, MArray a e m, IArray b e) => a i e -> m (b i e)+unsafeFreeze = freeze++unsafeNewArray_ :: (Ix i, MArray a e m) => (i, i) -> m (a i e)+unsafeNewArray_ = newArray_++unsafeRead :: (MArray a e m, Ix i) => a i e -> i -> m e+unsafeRead = readArray++unsafeThaw :: (Ix i, IArray a e, MArray b e m) => a i e -> m (b i e)+unsafeThaw = thaw++unsafeWrite :: (MArray a e m, Ix i) => a i e -> i -> e -> m ()+unsafeWrite = writeArray++#else++import Data.Array.Base+  ( UArray+  , bounds+  , castSTUArray+  , unsafeAt+  , unsafeFreeze+  , unsafeNewArray_+  , unsafeRead+  , unsafeThaw+  , unsafeWrite+  )++#endif

@@ -63,7 +63,7 @@ "shiftToOddCount/Integer"   shiftToOddCount = shiftOCInteger   #-} {-# INLINE [1] shiftToOddCount #-}-shiftToOddCount :: (Integral a, Bits a) => a -> (Int, a)+shiftToOddCount :: Integral a => a -> (Int, a) shiftToOddCount n = case shiftOCInteger (fromIntegral n) of                       (z, o) -> (z, fromInteger o) @@ -125,7 +125,7 @@ "shiftToOdd/Integer"   shiftToOdd = shiftOInteger   #-} {-# INLINE [1] shiftToOdd #-}-shiftToOdd :: (Integral a, Bits a) => a -> a+shiftToOdd :: Integral a => a -> a shiftToOdd n = fromInteger (shiftOInteger (fromIntegral n))  -- | Specialised version for @'Int'@.

@@ -1,5 +1,5 @@ name                : arithmoi-version             : 0.4.1.3+version             : 0.4.2.0 cabal-version       : >= 1.10 author              : Daniel Fischer copyright           : (c) 2011 Daniel Fischer@@ -27,12 +27,12 @@  category            : Math, Algorithms, Number Theory -tested-with         : GHC == 7.4.2, GHC==7.6.3, GHC==7.8.3+tested-with         : GHC==7.6.3, GHC==7.8.4, GHC==7.10.3, GHC==8.0.1  extra-source-files  : Changes, TODO -flag llvm-    description         : Compile the library with the LLVM backend+flag check-bounds+    description         : Replace unsafe array operations with safe ones     default             : False     manual              : True @@ -40,7 +40,7 @@     default-language: Haskell2010     build-depends       : base >= 4 && < 5                           , array >= 0.3 && < 0.6-                          , ghc-prim < 0.5+                          , ghc-prim < 0.6                           , integer-gmp < 1.1                           , containers >= 0.3 && < 0.6                           , random >= 1.0 && < 1.2@@ -51,6 +51,7 @@                           Math.NumberTheory.MoebiusInversion                           Math.NumberTheory.MoebiusInversion.Int                           Math.NumberTheory.Lucas+                          Math.NumberTheory.GaussianIntegers                           Math.NumberTheory.GCD                           Math.NumberTheory.GCD.LowLevel                           Math.NumberTheory.Powers@@ -68,6 +69,7 @@                           Math.NumberTheory.Primes.Testing.Certificates                           Math.NumberTheory.Primes.Heap     other-modules       : Math.NumberTheory.Utils+                          Math.NumberTheory.Unsafe                           Math.NumberTheory.Logarithms.Internal                           Math.NumberTheory.Primes.Counting.Impl                           Math.NumberTheory.Primes.Counting.Approximate@@ -83,9 +85,9 @@     other-extensions          : BangPatterns      ghc-options         : -O2 -Wall-    if flag(llvm)-        ghc-options     : -fllvm     ghc-prof-options    : -O2 -auto+    if flag(check-bounds)+        cpp-options     : -DCheckBounds  source-repository head   type:     git@@ -105,9 +107,31 @@   type:                 exitcode-stdio-1.0   hs-source-dirs:       test-suite   ghc-options:          -Wall-  main-is:              Spec.hs+  main-is:              Test.hs   default-language: Haskell2010-  build-depends:        base-                        ,hspec-                        ,arithmoi-+  build-depends:        base >= 4 && < 5+                        , arithmoi >= 0.4 && < 0.5+                        , tasty >= 0.10 && < 0.12+                        , tasty-smallcheck >= 0.8 && < 0.9+                        , tasty-quickcheck >= 0.8 && < 0.9+                        , tasty-hunit >= 0.9 && < 0.10+                        , QuickCheck >= 2.8 && < 2.9+                        , smallcheck >= 1.1 && < 1.2+  other-modules :   Math.NumberTheory.GaussianIntegersTests+                  , Math.NumberTheory.GCDTests+                  , Math.NumberTheory.GCD.LowLevelTests+                  , Math.NumberTheory.LogarithmsTests+                  , Math.NumberTheory.LucasTests+                  , 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.HeapTests+                  , Math.NumberTheory.Primes.SieveTests+                  , Math.NumberTheory.TestUtils

@@ -0,0 +1,74 @@+-- |+-- Module:      Math.NumberTheory.GCD.LowLevelTests+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.GCD.LowLevel+--++{-# LANGUAGE CPP       #-}+{-# LANGUAGE MagicHash #-}++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.GCD.LowLevelTests+  ( testSuite+  ) where++import Test.Tasty++#if MIN_VERSION_base(4,8,0)+#else+import Data.Word+#endif++import GHC.Exts++import Math.NumberTheory.GCD.LowLevel+import Math.NumberTheory.TestUtils++-- | Check that 'gcdInt' matches 'gcd'.+gcdIntProperty :: Int -> Int -> Bool+gcdIntProperty a b = gcdInt a b == gcd a b++-- | Check that 'gcdWord' matches 'gcd'.+gcdWordProperty :: Word -> Word -> Bool+gcdWordProperty a b = gcdWord a b == gcd a b++-- | Check that 'gcdInt#' matches 'gcd'.+gcdIntProperty# :: Int -> Int -> Bool+gcdIntProperty# a@(I# a') b@(I# b') = I# (gcdInt# a' b') == gcd a b++-- | Check that 'gcdWord#' matches 'gcd'.+gcdWordProperty# :: Word -> Word -> Bool+gcdWordProperty# a@(W# a') b@(W# b') = W# (gcdWord# a' b') == gcd a b++-- | Check that numbers are coprime iff their gcd equals to 1.+coprimeIntProperty :: Int -> Int -> Bool+coprimeIntProperty a b = coprimeInt a b == (gcd a b == 1)++-- | Check that numbers are coprime iff their gcd equals to 1.+coprimeWordProperty :: Word -> Word -> Bool+coprimeWordProperty a b = coprimeWord a b == (gcd a b == 1)++-- | Check that numbers are coprime iff their gcd equals to 1.+coprimeIntProperty# :: Int -> Int -> Bool+coprimeIntProperty# a@(I# a') b@(I# b') = coprimeInt# a' b' == (gcd a b == 1)++-- | Check that numbers are coprime iff their gcd equals to 1.+coprimeWordProperty# :: Word -> Word -> Bool+coprimeWordProperty# a@(W# a') b@(W# b') = coprimeWord# a' b' == (gcd a b == 1)++testSuite :: TestTree+testSuite = testGroup "LowLevel"+  [ testSmallAndQuick "gcdInt"       gcdIntProperty+  , testSmallAndQuick "gcdWord"      gcdWordProperty+  , testSmallAndQuick "gcdInt#"      gcdIntProperty#+  , testSmallAndQuick "gcdWord#"     gcdWordProperty#+  , testSmallAndQuick "coprimeInt"   coprimeIntProperty+  , testSmallAndQuick "coprimeWord"  coprimeWordProperty+  , testSmallAndQuick "coprimeInt#"  coprimeIntProperty#+  , testSmallAndQuick "coprimeWord#" coprimeWordProperty#+  ]

@@ -0,0 +1,50 @@+-- |+-- Module:      Math.NumberTheory.GCDTests+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.GCD+--++{-# LANGUAGE ScopedTypeVariables #-}++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.GCDTests+  ( testSuite+  ) where++import Test.Tasty++import Data.Bits++import Math.NumberTheory.GCD+import Math.NumberTheory.TestUtils++-- | Check that 'binaryGCD' matches 'gcd'.+binaryGCDProperty :: (Integral a, Bits a) => AnySign a -> AnySign a -> Bool+binaryGCDProperty (AnySign a) (AnySign b) = binaryGCD a b == gcd a b++-- | Check that 'extendedGCD' is consistent with documentation.+extendedGCDProperty :: forall a. Integral a => AnySign a -> AnySign a -> Bool+extendedGCDProperty (AnySign a) (AnySign b) =+  u * a + v * b == d+  && d == gcd a b+  -- (-1) >= 0 is true for unsigned types+  && (abs u < abs b || abs b <= 1 || (-1 :: a) >= 0)+  && (abs v < abs a || abs a <= 1 || (-1 :: a) >= 0)+  where+    (d, u, v) = extendedGCD a b++-- | Check that numbers are coprime iff their gcd equals to 1.+coprimeProperty :: (Integral a, Bits a) => AnySign a -> AnySign a -> Bool+coprimeProperty (AnySign a) (AnySign b) = coprime a b == (gcd a b == 1)++testSuite :: TestTree+testSuite = testGroup "GCD"+  [ testSameIntegralProperty "binaryGCD"   binaryGCDProperty+  , testSameIntegralProperty "extendedGCD" extendedGCDProperty+  , testSameIntegralProperty "coprime"     coprimeProperty+  ]

@@ -0,0 +1,79 @@+{-# OPTIONS_GHC -fno-warn-type-defaults #-}++-- |+-- Module:      Math.NumberTheory.GaussianIntegersTests+-- Copyright:   (c) 2016 Chris Fredrickson+-- Licence:     MIT+-- Maintainer:  Chris Fredrickson <chris.p.fredrickson@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.GaussianIntegers+--++module Math.NumberTheory.GaussianIntegersTests+  ( testSuite+  ) where++import Test.Tasty+import Test.Tasty.HUnit++import Math.NumberTheory.GaussianIntegers+import Math.NumberTheory.TestUtils++-- | Number is zero or is equal to the product of its factors.+factoriseProperty :: Integer -> Integer -> Bool+factoriseProperty x y+  =  x == 0 && y == 0+  || g == g'+  where+    g = x :+ y+    factors = factorise g+    g' = product $ map (uncurry (.^)) factors++-- | Number is prime iff it is non-zero+--   and has exactly one (non-unit) factor.+isPrimeProperty :: Integer -> Integer -> Bool+isPrimeProperty x y+  =  x == 0 && y == 0+  || isPrime g && n == 1+  || not (isPrime g) && n /= 1+  where+    g = x :+ y+    factors = factorise g+    nonUnitFactors = filter (\(p, _) -> norm p /= 1) factors+    -- Count factors taking into account multiplicity+    n = sum $ map snd nonUnitFactors++-- | The list of primes should include only primes.+primesGeneratesPrimesProperty :: NonNegative Int -> Bool+primesGeneratesPrimesProperty (NonNegative i) = isPrime (primes !! i)++-- | signum and abs should satisfy: z == signum z * abs z+signumAbsProperty :: Integer -> Integer -> Bool+signumAbsProperty x y = z == signum z * abs z+  where+    z = x :+ y++-- | abs maps a Gaussian integer to its associate in first quadrant.+absProperty :: Integer -> Integer -> Bool+absProperty x y = isOrigin || (inFirstQuadrant && isAssociate)+  where+    z = x :+ y+    z'@(x' :+ y') = abs z+    isOrigin = z' == 0 && z == 0+    inFirstQuadrant = x' > 0 && y' >= 0     -- first quadrant includes the positive real axis, but not the origin or the positive imaginary axis+    isAssociate = z' `elem` map (\e -> z * (0 :+ 1) .^ e) [0 .. 3]++-- | a special case that tests rounding/truncating in GCD.+gcdGSpecialCase1 :: Assertion+gcdGSpecialCase1 = assertEqual "gcdG" 1 $ gcdG (12 :+ 23) (23 :+ 34)++testSuite :: TestTree+testSuite = testGroup "GaussianIntegers"+  [ testSmallAndQuick "factorise"         factoriseProperty+  , testSmallAndQuick "isPrime"           isPrimeProperty+  , testSmallAndQuick "primes"            primesGeneratesPrimesProperty+  , testSmallAndQuick "signumAbsProperty" signumAbsProperty+  , testSmallAndQuick "absProperty"       absProperty+  , testCase          "gcdG (12 :+ 23) (23 :+ 34)" gcdGSpecialCase1+  ]

@@ -0,0 +1,112 @@+-- |+-- 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+  ]

@@ -0,0 +1,104 @@+-- |+-- 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+    ]+  ]

@@ -0,0 +1,222 @@+-- |+-- Module:      Math.NumberTheory.ModuliTests+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.Moduli+--++{-# LANGUAGE CPP          #-}+{-# LANGUAGE ViewPatterns #-}++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.ModuliTests+  ( testSuite+  ) where++import Test.Tasty++import Control.Arrow+import Data.Bits+import Data.List (tails, nub)+import Data.Maybe++import Math.NumberTheory.Moduli+import Math.NumberTheory.TestUtils++toOdd :: Num a => a -> a+toOdd n = n * 2 + 1++unwrapPP :: (Prime, Power Int) -> (Integer, Int)+unwrapPP (Prime p, Power e) = (p, e)++-- | Check that 'jacobi' matches 'jacobi''.+jacobiProperty1 :: (Integral a, Bits a) => AnySign a -> NonNegative a -> Bool+jacobiProperty1 (AnySign a) (NonNegative (toOdd -> n)) = n == 1 && j == 1 || n > 1 && j == j'+  where+    j = jacobi a n+    j' = jacobi' a n++-- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 2+jacobiProperty2 :: (Integral a, Bits a) => AnySign a -> NonNegative a -> Bool+jacobiProperty2 (AnySign a) (NonNegative (toOdd -> n)) = jacobi a n == jacobi (a + n) n++-- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 3+jacobiProperty3 :: (Integral a, Bits a) => AnySign a -> NonNegative a -> Bool+jacobiProperty3 (AnySign a) (NonNegative (toOdd -> n)) = j == 0 && g /= 1 || abs j == 1 && g == 1+  where+    j = jacobi a n+    g = gcd a n++-- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 4+jacobiProperty4 :: (Integral a, Bits a) => AnySign a -> AnySign a -> NonNegative a -> Bool+jacobiProperty4 (AnySign a) (AnySign b) (NonNegative (toOdd -> n)) = jacobi (a * b) n == jacobi a n * jacobi b n++jacobiProperty4_Integer :: AnySign Integer -> AnySign Integer -> NonNegative Integer -> Bool+jacobiProperty4_Integer = jacobiProperty4++-- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 5+jacobiProperty5 :: (Integral a, Bits a) => AnySign a -> NonNegative a -> NonNegative a -> Bool+jacobiProperty5 (AnySign a) (NonNegative (toOdd -> m)) (NonNegative (toOdd -> n)) = jacobi a (m * n) == jacobi a m * jacobi a n++jacobiProperty5_Integer :: AnySign Integer -> NonNegative Integer -> NonNegative Integer -> Bool+jacobiProperty5_Integer = jacobiProperty5++-- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 6+jacobiProperty6 :: (Integral a, Bits a) => NonNegative a -> NonNegative a -> Bool+jacobiProperty6 (NonNegative (toOdd -> m)) (NonNegative (toOdd -> n)) = gcd m n /= 1 || jacobi m n * jacobi n m == (if m `mod` 4 == 1 || n `mod` 4 == 1 then 1 else -1)++-- | Check that 'invertMod' inverts numbers modulo.+invertModProperty :: AnySign Integer -> Positive Integer -> Bool+invertModProperty (AnySign k) (Positive m) = case invertMod k m of+  Nothing  -> k `mod` m == 0 || gcd k m > 1+  Just inv -> gcd k m == 1+      && k * inv `mod` m == 1 && 0 <= inv && inv < m++-- | Check that the result of 'powerMod' is between 0 and modulo (non-inclusive).+powerModProperty1 :: (Integral a, Bits a) => AnySign a -> AnySign Integer -> Positive Integer -> Bool+powerModProperty1 (AnySign e) (AnySign b) (Positive m)+  =  e < 0 && isNothing (invertMod b m)+  || (0 <= pm && pm < m)+  where+    pm = powerMod b e m++-- | Check that 'powerMod' is multiplicative by first argument.+powerModProperty2 :: (Integral a, Bits a) => AnySign a -> AnySign Integer -> AnySign Integer -> Positive Integer -> Bool+powerModProperty2 (AnySign e) (AnySign b1) (AnySign b2) (Positive m)+  =  e < 0 && (isNothing (invertMod b1 m) || isNothing (invertMod b2 m))+  || pm1 * pm2 `mod` m == pm12+  where+    pm1  = powerMod b1  e m+    pm2  = powerMod b2  e m+    pm12 = powerMod (b1 * b2) e m++-- | Check that 'powerMod' is additive by second argument.+powerModProperty3 :: (Integral a, Bits a) => AnySign a -> AnySign a -> AnySign Integer -> Positive Integer -> Bool+powerModProperty3 (AnySign e1) (AnySign e2) (AnySign b) (Positive m)+  =  (e1 < 0 || e2 < 0) && isNothing (invertMod b m)+  || pm1 * pm2 `mod` m == pm12+  where+    pm1  = powerMod b e1 m+    pm2  = powerMod b e2 m+    pm12 = powerMod b (e1 + e2) m++-- | Specialized to trigger 'powerModInteger'.+powerModProperty1_Integer :: AnySign Integer -> AnySign Integer -> Positive Integer -> Bool+powerModProperty1_Integer = powerModProperty1++-- | Specialized to trigger 'powerModInteger'.+powerModProperty2_Integer :: AnySign Integer -> AnySign Integer -> AnySign Integer -> Positive Integer -> Bool+powerModProperty2_Integer = powerModProperty2++-- | Specialized to trigger 'powerModInteger'.+powerModProperty3_Integer :: AnySign Integer -> AnySign Integer -> AnySign Integer -> Positive Integer -> Bool+powerModProperty3_Integer = powerModProperty3++-- | Check that 'powerMod' matches 'powerMod''.+powerMod'Property :: (Integral a, Bits a) => Positive a -> Positive Integer -> Positive Integer -> Bool+powerMod'Property (Positive e) (Positive b) (Positive m) = m == 1 || powerMod' b e m == powerMod b e m++-- | Specialized to trigger 'powerModInteger''.+powerMod'Property_Integer :: Positive Integer -> Positive Integer -> Positive Integer -> Bool+powerMod'Property_Integer = powerMod'Property++-- | Check that 'chineseRemainder' is defined iff modulos are coprime.+--   Also check that the result is a solution of input modular equations.+chineseRemainderProperty :: [(Integer, Positive Integer)] -> Bool+chineseRemainderProperty rms' = case chineseRemainder rms of+  Nothing -> not areCoprime+  Just n  -> areCoprime && map (n `mod`) ms == zipWith mod rs ms+  where+    rms = map (second getPositive) rms'+    (rs, ms) = unzip rms+    areCoprime = all (== 1) [ gcd m1 m2 | (m1 : m2s) <- tails ms, m2 <- m2s ]++-- | Check that 'chineseRemainder' matches 'chineseRemainder2'.+chineseRemainder2Property :: Integer -> Positive Integer -> Integer -> Positive Integer -> Bool+chineseRemainder2Property r1 (Positive m1) r2 (Positive m2) = gcd m1 m2 /= 1+  || Just (chineseRemainder2 (r1, m1) (r2, m2)) == chineseRemainder [(r1, m1), (r2, m2)]++-- | Check that 'sqrtMod' is defined iff a quadratic residue exists.+--   Also check that the result is a solution of input modular equation.+sqrtModPProperty :: AnySign Integer -> Prime -> Bool+sqrtModPProperty (AnySign n) (Prime p) = case sqrtModP n p of+  Nothing -> jacobi n p == -1+  Just rt -> (p == 2 || jacobi n p /= -1) && rt ^ 2 `mod` p == n `mod` p++sqrtModPListProperty :: AnySign Integer -> Prime -> Bool+sqrtModPListProperty (AnySign n) (Prime p) = all (\rt -> rt ^ 2 `mod` p == n `mod` p) (sqrtModPList n p)++sqrtModP'Property :: Positive Integer -> Prime -> Bool+sqrtModP'Property (Positive n) (Prime p) = (p /= 2 && jacobi n p /= 1) || rt ^ 2 `mod` p == n `mod` p+  where+    rt = sqrtModP' n p++tonelliShanksProperty :: Positive Integer -> Prime -> Bool+tonelliShanksProperty (Positive n) (Prime p) = p `mod` 4 /= 1 || jacobi n p /= 1 || rt ^ 2 `mod` p == n `mod` p+  where+    rt = tonelliShanks n p++sqrtModPPProperty :: AnySign Integer -> (Prime, Power Int) -> Bool+sqrtModPPProperty (AnySign n) (Prime p, Power e) = gcd n p > 1 || case sqrtModPP n (p, e) of+  Nothing -> True+  Just rt -> rt ^ 2 `mod` (p ^ e) == n `mod` (p ^ e)++sqrtModPPListProperty :: AnySign Integer -> (Prime, Power Int) -> Bool+sqrtModPPListProperty (AnySign n) (Prime p, Power e) = gcd n p > 1+  || all (\rt -> rt ^ 2 `mod` (p ^ e) == n `mod` (p ^ e)) (sqrtModPPList n (p, e))++sqrtModFProperty :: AnySign Integer -> [(Prime, Power Int)] -> Bool+sqrtModFProperty (AnySign n) (map unwrapPP -> pes) = case sqrtModF n pes of+  Nothing -> True+  Just rt -> all (\(p, e) -> rt ^ 2 `mod` (p ^ e) == n `mod` (p ^ e)) pes++sqrtModFListProperty :: AnySign Integer -> [(Prime, Power Int)] -> Bool+sqrtModFListProperty (AnySign n) (map unwrapPP -> pes)+  = nub ps /= ps || all+    (\rt -> all (\(p, e) -> rt ^ 2 `mod` (p ^ e) == n `mod` (p ^ e)) pes)+    (sqrtModFList n pes)+  where+    ps = map fst pes++testSuite :: TestTree+testSuite = testGroup "Moduli"+  [ testGroup "jacobi"+    [ testSameIntegralProperty "matches jacobi'"              jacobiProperty1+    , testSameIntegralProperty "same modulo n"                jacobiProperty2+    , testSameIntegralProperty "consistent with gcd"          jacobiProperty3+    , testSmallAndQuick        "multiplicative 1"             jacobiProperty4_Integer+    , testSmallAndQuick        "multiplicative 2"             jacobiProperty5_Integer+    , testSameIntegralProperty "law of quadratic reciprocity" jacobiProperty6+    ]+  , testSmallAndQuick "invertMod" invertModProperty+  , testGroup "powerMod"+    [ testGroup "generic"+      [ testIntegralProperty "bounded between 0 and m"  powerModProperty1+      , testIntegralProperty "multiplicative by base"   powerModProperty2+      , testSameIntegralProperty "additive by exponent" powerModProperty3+      , testIntegralProperty "matches powerMod'"        powerMod'Property+      ]+    , testGroup "Integer"+      [ testSmallAndQuick "bounded between 0 and m" powerModProperty1_Integer+      , testSmallAndQuick "multiplicative by base"  powerModProperty2_Integer+      , testSmallAndQuick "additive by exponent"    powerModProperty3_Integer+      , testSmallAndQuick "matches powerMod'"       powerMod'Property_Integer+      ]+    ]+    , testSmallAndQuick "chineseRemainder"  chineseRemainderProperty+    , testSmallAndQuick "chineseRemainder2" chineseRemainder2Property+    , testGroup "sqrtMod"+      [ testSmallAndQuick "sqrtModP"      sqrtModPProperty+      , testSmallAndQuick "sqrtModPList"  sqrtModPListProperty+      , testSmallAndQuick "sqrtModP'"     sqrtModP'Property+      , testSmallAndQuick "tonelliShanks" tonelliShanksProperty+      , testSmallAndQuick "sqrtModPP"     sqrtModPPProperty+      , testSmallAndQuick "sqrtModPPList" sqrtModPPListProperty+      , testSmallAndQuick "sqrtModF"      sqrtModFProperty+      , testSmallAndQuick "sqrtModFList"  sqrtModFListProperty+      ]+  ]

@@ -0,0 +1,45 @@+-- |+-- Module:      Math.NumberTheory.MoebiusInversion.IntTests+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.MoebiusInversion.Int+--++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.MoebiusInversion.IntTests+  ( testSuite+  ) where++import Test.Tasty+import Test.Tasty.HUnit+import Test.Tasty.QuickCheck as QC hiding (Positive)++import Math.NumberTheory.MoebiusInversion.Int+import Math.NumberTheory.Primes.Factorisation+import Math.NumberTheory.TestUtils++totientSumProperty :: Positive Int -> Bool+totientSumProperty (Positive n) = toInteger (totientSum n) == sum (map totient [1 .. toInteger n])++totientSumSpecialCase1 :: Assertion+totientSumSpecialCase1 = assertEqual "totientSum" 4496 (totientSum 121)++generalInversionProperty :: (Int -> Int) -> Positive Int -> Bool+generalInversionProperty g (Positive n)+  =  g n == sum [f (n `quot` k) | k <- [1 .. n]]+  && f n == sum [fromInteger (moebius (toInteger k)) * g (n `quot` k) | k <- [1 .. n]]+  where+    f = generalInversion g++testSuite :: TestTree+testSuite = testGroup "Int"+  [ testGroup "totientSum"+    [ testSmallAndQuick "matches definitions" totientSumProperty+    , testCase          "special case 1"      totientSumSpecialCase1+    ]+  , QC.testProperty "generalInversion" generalInversionProperty+  ]

@@ -0,0 +1,45 @@+-- |+-- Module:      Math.NumberTheory.MoebiusInversionTests+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.MoebiusInversion+--++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.MoebiusInversionTests+  ( testSuite+  ) where++import Test.Tasty+import Test.Tasty.HUnit+import Test.Tasty.QuickCheck as QC hiding (Positive)++import Math.NumberTheory.MoebiusInversion+import Math.NumberTheory.Primes.Factorisation+import Math.NumberTheory.TestUtils++totientSumProperty :: Positive Int -> Bool+totientSumProperty (Positive n) = totientSum n == sum (map totient [1 .. toInteger n])++totientSumSpecialCase1 :: Assertion+totientSumSpecialCase1 = assertEqual "totientSum" 4496 (totientSum 121)++generalInversionProperty :: (Int -> Integer) -> Positive Int -> Bool+generalInversionProperty g (Positive n)+  =  g n == sum [f (n `quot` k) | k <- [1 .. n]]+  && f n == sum [moebius (toInteger k) * g (n `quot` k) | k <- [1 .. n]]+  where+    f = generalInversion g++testSuite :: TestTree+testSuite = testGroup "MoebiusInversion"+  [ testGroup "totientSum"+    [ testSmallAndQuick "matches definitions" totientSumProperty+    , testCase          "special case 1"      totientSumSpecialCase1+    ]+  , QC.testProperty "generalInversion" generalInversionProperty+  ]

@@ -0,0 +1,154 @@+-- |+-- Module:      Math.NumberTheory.Powers.CubesTests+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.Powers.Cubes+--++{-# LANGUAGE CPP #-}++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.Powers.CubesTests+  ( testSuite+  ) where++import Test.Tasty+import Test.Tasty.HUnit++import Data.Maybe+#if MIN_VERSION_base(4,8,0)+#else+import Data.Word+#endif+++import Math.NumberTheory.Powers.Cubes+import Math.NumberTheory.TestUtils++#include "MachDeps.h"++-- | Check that 'integerCubeRoot' returns the largest integer @m@ with @m^3 <= n@.+--+-- (m + 1) ^ 3 /= n && cond+-- means+-- (m + 1) ^ 3 > n+-- but without overflow for bounded types+integerCubeRootProperty :: Integral a => AnySign a -> Bool+integerCubeRootProperty (AnySign n) = m ^ 3 <= n && (m + 1) ^ 3 /= n && cond+  where+    m = integerCubeRoot n+    cond+      | m == -1   = n == -1+      | m < 0     = (m + 1) ^ 2 <= n `div` (m + 1)+      | otherwise = (m + 1) ^ 2 >= n `div` (m + 1)++-- | Specialized to trigger 'cubeRootInt''.+integerCubeRootProperty_Int :: AnySign Int -> Bool+integerCubeRootProperty_Int = integerCubeRootProperty++-- | Specialized to trigger 'cubeRootWord'.+integerCubeRootProperty_Word :: AnySign Word -> Bool+integerCubeRootProperty_Word = integerCubeRootProperty++-- | Specialized to trigger 'cubeRootIgr'.+integerCubeRootProperty_Integer :: AnySign Integer -> Bool+integerCubeRootProperty_Integer = integerCubeRootProperty++-- | Check that 'integerCubeRoot' returns the largest integer @m@ with @m^3 <= n@, , where @n@ has form @k@^3-1.+integerCubeRootProperty2 :: Integral a => AnySign a -> Bool+integerCubeRootProperty2 (AnySign k) = m ^ 3 <= n && (m + 1) ^ 3 /= n && cond+  where+    n = k ^ 3 - 1+    m = integerCubeRoot n+    cond+      | m == -1   = n == -1+      | m < 0     = (m + 1) ^ 2 <= n `div` (m + 1)+      | otherwise = (m + 1) ^ 2 >= n `div` (m + 1)++-- | Specialized to trigger 'cubeRootInt''.+integerCubeRootProperty2_Int :: AnySign Int -> Bool+integerCubeRootProperty2_Int = integerCubeRootProperty2++-- | Specialized to trigger 'cubeRootWord'.+integerCubeRootProperty2_Word :: AnySign Word -> Bool+integerCubeRootProperty2_Word = integerCubeRootProperty2++#if WORD_SIZE_IN_BITS == 64++-- | Check that 'integerCubeRoot' of 2^63-1 is 2^21-1, not 2^21.+integerCubeRootSpecialCase1_Int :: Assertion+integerCubeRootSpecialCase1_Int =+  assertEqual "integerCubeRoot" (integerCubeRoot (maxBound :: Int)) (2 ^ 21 - 1)++-- | Check that 'integerCubeRoot' of 2^63-1 is 2^21-1, not 2^21.+integerCubeRootSpecialCase1_Word :: Assertion+integerCubeRootSpecialCase1_Word =+  assertEqual "integerCubeRoot" (integerCubeRoot (maxBound `div` 2 :: Word)) (2 ^ 21 - 1)++-- | Check that 'integerCubeRoot' of 2^64-1 is 2642245.+integerCubeRootSpecialCase2 :: Assertion+integerCubeRootSpecialCase2 =+  assertEqual "integerCubeRoot" (integerCubeRoot (maxBound :: Word)) 2642245++#endif++-- | Check that 'integerCubeRoot'' returns the largest integer @m@ with @m^3 <= n@.+integerCubeRoot'Property :: Integral a => NonNegative a -> Bool+integerCubeRoot'Property (NonNegative n) = m ^ 3 <= n && (m + 1) ^ 3 /= n && (m + 1) ^ 2 >= n `div` (m + 1)+  where+    m = integerCubeRoot' n++-- | Check that the number 'isCube' iff its 'integerCubeRoot' is exact.+isCubeProperty :: Integral a => AnySign a -> Bool+isCubeProperty (AnySign n) = (n /= m ^ 3 && not t) || (n == m ^ 3 && t)+  where+    t = isCube n+    m = integerCubeRoot n++-- | Check that the number 'isCube'' iff its 'integerCubeRoot'' is exact.+isCube'Property :: Integral a => NonNegative a -> Bool+isCube'Property (NonNegative n) = (n /= m ^ 3 && not t) || (n == m ^ 3 && t)+  where+    t = isCube' n+    m = integerCubeRoot' n++-- | Check that 'exactCubeRoot' returns an exact integer cubic root+-- and is consistent with 'isCube'.+exactCubeRootProperty :: Integral a => AnySign a -> Bool+exactCubeRootProperty (AnySign n) = case exactCubeRoot n of+  Nothing -> not (isCube n)+  Just m  -> isCube n && n == m ^ 3++-- | Check that 'isPossibleCube' is consistent with 'exactCubeRoot'.+isPossibleCubeProperty :: Integral a => NonNegative a -> Bool+isPossibleCubeProperty (NonNegative n) = t || not t && isNothing m+  where+    t = isPossibleCube n+    m = exactCubeRoot n++testSuite :: TestTree+testSuite = testGroup "Cubes"+  [ testGroup "integerCubeRoot"+    [ testIntegralProperty "generic"         integerCubeRootProperty+    , testSmallAndQuick    "generic Int"     integerCubeRootProperty_Int+    , testSmallAndQuick    "generic Word"    integerCubeRootProperty_Word+    , testSmallAndQuick    "generic Integer" integerCubeRootProperty_Integer++    , testIntegralProperty "almost cube"      integerCubeRootProperty2+    , testSmallAndQuick    "almost cube Int"  integerCubeRootProperty2_Int+    , testSmallAndQuick    "almost cube Word" integerCubeRootProperty2_Word++    , testCase             "maxBound :: Int"      integerCubeRootSpecialCase1_Int+    , testCase             "maxBound / 2 :: Word" integerCubeRootSpecialCase1_Word+    , testCase             "maxBound :: Word"     integerCubeRootSpecialCase2+    ]+  , testIntegralProperty "integerCubeRoot'" integerCubeRoot'Property+  , testIntegralProperty "isCube"           isCubeProperty+  , testIntegralProperty "isCube'"          isCube'Property+  , testIntegralProperty "exactCubeRoot"    exactCubeRootProperty+  , testIntegralProperty "isPossibleCube"   isPossibleCubeProperty+  ]

@@ -0,0 +1,145 @@+-- |+-- Module:      Math.NumberTheory.Powers.FourthTests+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.Powers.Fourth+--++{-# LANGUAGE CPP #-}++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.Powers.FourthTests+  ( testSuite+  ) where++import Test.Tasty+import Test.Tasty.HUnit++import Data.Maybe+#if MIN_VERSION_base(4,8,0)+#else+import Data.Word+#endif++import Math.NumberTheory.Powers.Fourth+import Math.NumberTheory.TestUtils++#include "MachDeps.h"++-- | Check that 'integerFourthRoot' returns the largest integer @m@ with @m^4 <= n@.+--+-- (m + 1) ^ 4 /= n && (m + 1) ^ 3 >= n `div` (m + 1)+-- means+-- (m + 1) ^ 4 > n+-- but without overflow for bounded types+integerFourthRootProperty :: Integral a => NonNegative a -> Bool+integerFourthRootProperty (NonNegative n) = m >= 0 && m ^ 4 <= n && (m + 1) ^ 4 /= n && (m + 1) ^ 3 >= n `div` (m + 1)+  where+    m = integerFourthRoot n++-- | Specialized to trigger 'biSqrtInt'.+integerFourthRootProperty_Int :: NonNegative Int -> Bool+integerFourthRootProperty_Int = integerFourthRootProperty++-- | Specialized to trigger 'biSqrtWord'.+integerFourthRootProperty_Word :: NonNegative Word -> Bool+integerFourthRootProperty_Word = integerFourthRootProperty++-- | Specialized to trigger 'biSqrtIgr'.+integerFourthRootProperty_Integer :: NonNegative Integer -> Bool+integerFourthRootProperty_Integer = integerFourthRootProperty++-- | Check that 'integerFourthRoot' returns the largest integer @m@ with @m^4 <= n@, , where @n@ has form @k@^4-1.+integerFourthRootProperty2 :: Integral a => NonNegative a -> Bool+integerFourthRootProperty2 (NonNegative k) = n < 0 || m >= 0 && m ^ 4 <= n && (m + 1) ^ 4 /= n && (m + 1) ^ 3 >= n `div` (m + 1)+  where+    n = k ^ 4 - 1+    m = integerFourthRoot n++-- | Specialized to trigger 'biSqrtInt.+integerFourthRootProperty2_Int :: NonNegative Int -> Bool+integerFourthRootProperty2_Int = integerFourthRootProperty2++-- | Specialized to trigger 'biSqrtWord'.+integerFourthRootProperty2_Word :: NonNegative Word -> Bool+integerFourthRootProperty2_Word = integerFourthRootProperty2++#if WORD_SIZE_IN_BITS == 64++-- | Check that 'integerFourthRoot' of 2^60-1 is 2^15-1, not 2^15.+integerFourthRootSpecialCase1_Int :: Assertion+integerFourthRootSpecialCase1_Int =+  assertEqual "integerFourthRoot" (integerFourthRoot (maxBound `div` 8 :: Int)) (2 ^ 15 - 1)++-- | Check that 'integerFourthRoot' of 2^60-1 is 2^15-1, not 2^15.+integerFourthRootSpecialCase1_Word :: Assertion+integerFourthRootSpecialCase1_Word =+  assertEqual "integerFourthRoot" (integerFourthRoot (maxBound `div` 16 :: Word)) (2 ^ 15 - 1)++-- | Check that 'integerFourthRoot' of 2^64-1 is 2^16-1, not 2^16.+integerFourthRootSpecialCase2 :: Assertion+integerFourthRootSpecialCase2 =+  assertEqual "integerFourthRoot" (integerFourthRoot (maxBound :: Word)) (2 ^ 16 - 1)++#endif++-- | Check that 'integerFourthRoot'' returns the largest integer @m@ with @m^4 <= n@.+integerFourthRoot'Property :: Integral a => NonNegative a -> Bool+integerFourthRoot'Property (NonNegative n) = m >= 0 && m ^ 4 <= n && (m + 1) ^ 4 /= n && (m + 1) ^ 3 >= n `div` (m + 1)+  where+    m = integerFourthRoot' n++-- | Check that the number 'isFourthPower' iff its 'integerFourthRoot' is exact.+isFourthPowerProperty :: Integral a => AnySign a -> Bool+isFourthPowerProperty (AnySign n) = (n < 0 && not t) || (n /= m ^ 4 && not t) || (n == m ^ 4 && t)+  where+    t = isFourthPower n+    m = integerFourthRoot n++-- | Check that the number 'isFourthPower'' iff its 'integerFourthRoot'' is exact.+isFourthPower'Property :: Integral a => NonNegative a -> Bool+isFourthPower'Property (NonNegative n) = (n /= m ^ 4 && not t) || (n == m ^ 4 && t)+  where+    t = isFourthPower' n+    m = integerFourthRoot' n++-- | Check that 'exactFourthRoot' returns an exact integer root of fourth power+-- and is consistent with 'isFourthPower'.+exactFourthRootProperty :: Integral a => AnySign a -> Bool+exactFourthRootProperty (AnySign n) = case exactFourthRoot n of+  Nothing -> not (isFourthPower n)+  Just m  -> isFourthPower n && n == m ^ 4++-- | Check that 'isPossibleFourthPower' is consistent with 'exactFourthRoot'.+isPossibleFourthPowerProperty :: Integral a => NonNegative a -> Bool+isPossibleFourthPowerProperty (NonNegative n) = t || not t && isNothing m+  where+    t = isPossibleFourthPower n+    m = exactFourthRoot n++testSuite :: TestTree+testSuite = testGroup "Fourth"+  [ testGroup "integerFourthRoot"+    [ testIntegralProperty "generic"         integerFourthRootProperty+    , testSmallAndQuick    "generic Int"     integerFourthRootProperty_Int+    , testSmallAndQuick    "generic Word"    integerFourthRootProperty_Word+    , testSmallAndQuick    "generic Integer" integerFourthRootProperty_Integer++    , testIntegralProperty "almost Fourth"      integerFourthRootProperty2+    , testSmallAndQuick    "almost Fourth Int"  integerFourthRootProperty2_Int+    , testSmallAndQuick    "almost Fourth Word" integerFourthRootProperty2_Word++    , testCase             "maxBound / 8 :: Int"   integerFourthRootSpecialCase1_Int+    , testCase             "maxBound / 16 :: Word" integerFourthRootSpecialCase1_Word+    , testCase             "maxBound :: Word"      integerFourthRootSpecialCase2+    ]+  , testIntegralProperty "integerFourthRoot'"    integerFourthRoot'Property+  , testIntegralProperty "isFourthPower"         isFourthPowerProperty+  , testIntegralProperty "isFourthPower'"        isFourthPower'Property+  , testIntegralProperty "exactFourthRoot"       exactFourthRootProperty+  , testIntegralProperty "isPossibleFourthPower" isPossibleFourthPowerProperty+  ]

@@ -0,0 +1,128 @@+-- |+-- Module:      Math.NumberTheory.Powers.GeneralTests+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.Powers.General+--++{-# LANGUAGE CPP #-}++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.Powers.GeneralTests+  ( testSuite+  ) where++import Test.Tasty+import Test.Tasty.HUnit++import Math.NumberTheory.Powers.General+import Math.NumberTheory.TestUtils++-- | Check that 'integerRoot' @pow@ returns the largest integer @m@ with @m^pow <= n@.+integerRootProperty :: (Integral a, Integral b) => AnySign a -> Power b -> Bool+integerRootProperty (AnySign n) (Power pow) = (even pow && n < 0)+  || (toInteger root ^ pow <= toInteger n && toInteger n < toInteger (root + 1) ^ pow)+    where+      root = integerRoot pow n++-- | Check that the number 'isKthPower' iff its 'integerRoot' is exact.+isKthPowerProperty :: (Integral a, Integral b) => AnySign a -> Power b -> Bool+isKthPowerProperty (AnySign n) (Power pow) = (even pow && n < 0 && not t) || (n /= root ^ pow && not t) || (n == root ^ pow && t)+  where+    t = isKthPower pow n+    root = integerRoot pow n++-- | Check that 'exactRoot' returns an exact integer root+-- and is consistent with 'isKthPower'.+exactRootProperty :: (Integral a, Integral b) => AnySign a -> Power b -> Bool+exactRootProperty (AnySign n) (Power pow) = case exactRoot pow n of+  Nothing   -> not (isKthPower pow n)+  Just root -> isKthPower pow n && n == root ^ pow++-- | Check that 'isPerfectPower' is consistent with 'highestPower'.+isPerfectPowerProperty :: Integral a => AnySign a -> Bool+isPerfectPowerProperty (AnySign n) = (k > 1 && t) || (k == 1 && not t)+  where+    t = isPerfectPower n+    (_, k) = highestPower n++-- | Check that the first component of 'highestPower' is square-free.+highestPowerProperty :: Integral a => AnySign a -> Bool+highestPowerProperty (AnySign n) = (n `elem` [-1, 0, 1] && k == 3) || (b ^ k == n && b' == b && k' == 1)+  where+    (b, k) = highestPower n+    (b', k') = highestPower b++-- | Check that 'largePFPower' is consistent with documentation.+largePFPowerProperty :: Positive Integer -> Integer -> Bool+largePFPowerProperty (Positive bd) n = bd == 1 || b == 0 || d' /= 0 || n <= b * d * d || any (\p -> gcd n p > 1) [2..bd] || b ^ k == n+  where+    (b, k) = largePFPower bd n+    (d, d') = bd `quotRem` b++highestPowerSpecialCases :: [Assertion]+highestPowerSpecialCases =+  -- Freezes before d44a13b.+  [ a ( 1013582159576576+      , 1013582159576576+      , 1)+  -- Freezes before d44a13b.+  , a ( 1013582159576576 ^ 7+      , 1013582159576576+      , 7)++  , a ( -2 ^ 63 :: Int+      , -2 :: Int+      , 63)++  , a ( (2 ^ 63 - 1) ^ 21+      , 2 ^ 63 - 1+      , 21)++  , a ( 576116746989720969230211509779286598589421531472851155101032940901763389787901933902294777750323196846498573545522289802689311975294763847414975335235584+      , 576116746989720969230211509779286598589421531472851155101032940901763389787901933902294777750323196846498573545522289802689311975294763847414975335235584+      , 1)++  , a ( -340282366920938463500268095579187314689+      , -340282366920938463500268095579187314689+      , 1)++  , a ( 268398749 :: Int+      , 268398749 :: Int+      , 1)++  , a ( 118372752099 :: Int+      , 118372752099 :: Int+      , 1)++  , a ( 1409777209 :: Int+      , 37547 :: Int+      , 2)++  , a ( -6277101735386680764856636523970481806547819498980467802113+      , -18446744073709551617+      , 3)++  , a ( -18446744073709551619 ^ 5+      , -18446744073709551619+      , 5)+  ]+  where+    a (n, b, k) = assertEqual "highestPower" (b, k) (highestPower n)++testSuite :: TestTree+testSuite = testGroup "General"+  [ testIntegral2Property "integerRoot"    integerRootProperty+  , testIntegral2Property "isKthPower"     isKthPowerProperty+  , testIntegral2Property "exactRoot"      exactRootProperty+  , testIntegralProperty  "isPerfectPower" isPerfectPowerProperty+  , testGroup "highestPower"+    ( testIntegralProperty  "highestPower"   highestPowerProperty+    : zipWith (\i a -> testCase ("special case " ++ show i) a) [1..] highestPowerSpecialCases+    )+  , testSmallAndQuick     "largePFPower"   largePFPowerProperty+  ]

@@ -0,0 +1,41 @@+-- |+-- 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,155 @@+-- |+-- Module:      Math.NumberTheory.Powers.SquaresTests+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.Powers.Squares+--++{-# LANGUAGE CPP #-}++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.Powers.SquaresTests+  ( testSuite+  ) where++import Test.Tasty+import Test.Tasty.HUnit++import Data.Maybe+#if MIN_VERSION_base(4,8,0)+#else+import Data.Word+#endif++import Math.NumberTheory.Powers.Squares+import Math.NumberTheory.TestUtils++#include "MachDeps.h"++-- | Check that 'integerSquareRoot' returns the largest integer @m@ with @m*m <= n@.+--+-- (m + 1) ^ 2 /= n && m + 1 >= n `div` (m + 1)+-- means+-- (m + 1) ^ 2 > n+-- but without overflow for bounded types+integerSquareRootProperty :: Integral a => NonNegative a -> Bool+integerSquareRootProperty (NonNegative n) = m >=0 && m * m <= n && (m + 1) ^ 2 /= n && m + 1 >= n `div` (m + 1)+  where+    m = integerSquareRoot n++-- | Specialized to trigger 'isqrtInt''.+integerSquareRootProperty_Int :: NonNegative Int -> Bool+integerSquareRootProperty_Int = integerSquareRootProperty++-- | Specialized to trigger 'isqrtWord'.+integerSquareRootProperty_Word :: NonNegative Word -> Bool+integerSquareRootProperty_Word = integerSquareRootProperty++-- | Check that 'integerSquareRoot' returns the largest integer @m@ with @m*m <= n@, where @n@ has form @k@^2-1.+integerSquareRootProperty2 :: Integral a => NonNegative a -> Bool+integerSquareRootProperty2 (NonNegative k) = n < 0+  || m >=0 && m * m <= n && (m + 1) ^ 2 /= n && m + 1 >= n `div` (m + 1)+  where+    n = k ^ 2 - 1+    m = integerSquareRoot n++-- | Specialized to trigger 'isqrtInt''.+integerSquareRootProperty2_Int :: NonNegative Int -> Bool+integerSquareRootProperty2_Int = integerSquareRootProperty2++-- | Specialized to trigger 'isqrtWord'.+integerSquareRootProperty2_Word :: NonNegative Word -> Bool+integerSquareRootProperty2_Word = integerSquareRootProperty2++#if WORD_SIZE_IN_BITS == 64++-- | Check that 'integerSquareRoot' of 2^62-1 is 2^31-1, not 2^31.+integerSquareRootSpecialCase1_Int :: Assertion+integerSquareRootSpecialCase1_Int =+  assertEqual "integerSquareRoot" (integerSquareRoot (maxBound `div` 2 :: Int)) (2 ^ 31 - 1)++-- | Check that 'integerSquareRoot' of 2^62-1 is 2^31-1, not 2^31.+integerSquareRootSpecialCase1_Word :: Assertion+integerSquareRootSpecialCase1_Word =+  assertEqual "integerSquareRoot" (integerSquareRoot (maxBound `div` 4 :: Word)) (2 ^ 31 - 1)++-- | Check that 'integerSquareRoot' of 2^64-1 is 2^32-1, not 2^32.+integerSquareRootSpecialCase2 :: Assertion+integerSquareRootSpecialCase2 =+  assertEqual "integerSquareRoot" (integerSquareRoot (maxBound :: Word)) (2 ^ 32 - 1)++#endif++-- | Check that 'integerSquareRoot'' returns the largest integer @r@ with @r*r <= n@.+integerSquareRoot'Property :: Integral a => NonNegative a -> Bool+integerSquareRoot'Property (NonNegative n) = m >=0 && m * m <= n && (m + 1) ^ 2 /= n && m + 1 >= n `div` (m + 1)+  where+    m = integerSquareRoot' n++-- | Check that the number 'isSquare' iff its 'integerSquareRoot' is exact.+isSquareProperty :: Integral a => AnySign a -> Bool+isSquareProperty (AnySign n) = (n < 0 && not t) || (n /= m * m && not t) || (n == m * m && t)+  where+    t = isSquare n+    m = integerSquareRoot n++-- | Check that the number 'isSquare'' iff its 'integerSquareRoot'' is exact.+isSquare'Property :: Integral a => NonNegative a -> Bool+isSquare'Property (NonNegative n) = (n /= m * m && not t) || (n == m * m && t)+  where+    t = isSquare' n+    m = integerSquareRoot' n++-- | Check that 'exactSquareRoot' returns an exact integer square root+-- and is consistent with 'isSquare'.+exactSquareRootProperty :: Integral a => AnySign a -> Bool+exactSquareRootProperty (AnySign n) = case exactSquareRoot n of+  Nothing -> not (isSquare n)+  Just m  -> isSquare n && n == m * m++-- | Check that 'isPossibleSquare' is consistent with 'exactSquareRoot'+-- and that 'isPossibleSquare2' is a refinement of 'isPossibleSquare'.+isPossibleSquareProperty :: Integral a => NonNegative a -> Bool+isPossibleSquareProperty (NonNegative n) = t || not t && not t2 && isNothing m+  where+    t = isPossibleSquare n+    t2 = isPossibleSquare2 n+    m = exactSquareRoot n++-- | Check that 'isPossibleSquare2'' is consistent with 'exactSquareRoot'.+isPossibleSquare2Property :: Integral a => NonNegative a -> Bool+isPossibleSquare2Property (NonNegative n) = t || not t && isNothing m+  where+    t = isPossibleSquare2 n+    m = exactSquareRoot n+++testSuite :: TestTree+testSuite = testGroup "Squares"+  [ testGroup "integerSquareRoor"+    [ testIntegralProperty "generic"          integerSquareRootProperty+    , testSmallAndQuick    "generic Int"      integerSquareRootProperty_Int+    , testSmallAndQuick    "generic Word"     integerSquareRootProperty_Word++    , testIntegralProperty "almost square"      integerSquareRootProperty2+    , testSmallAndQuick    "almost square Int"  integerSquareRootProperty2_Int+    , testSmallAndQuick    "almost square Word" integerSquareRootProperty2_Word++#if WORD_SIZE_IN_BITS == 64+    , testCase             "maxBound / 2 :: Int"  integerSquareRootSpecialCase1_Int+    , testCase             "maxBound / 4 :: Word" integerSquareRootSpecialCase1_Word+    , testCase             "maxBound :: Word"     integerSquareRootSpecialCase2+#endif+    ]++  , testIntegralProperty "integerSquareRoot'" integerSquareRoot'Property+  , testIntegralProperty "isSquare"           isSquareProperty+  , testIntegralProperty "isSquare'"          isSquare'Property+  , testIntegralProperty "exactSquareRoot"    exactSquareRootProperty+  , testIntegralProperty "isPossibleSquare"   isPossibleSquareProperty+  , testIntegralProperty "isPossibleSquare2"  isPossibleSquare2Property+  ]

@@ -0,0 +1,148 @@+-- |+-- Module:      Math.NumberTheory.Primes.CountingTests+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.Primes.Counting+--++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.Primes.CountingTests+  ( testSuite+  ) where++import Test.Tasty+import Test.Tasty.HUnit++import Math.NumberTheory.Primes.Counting+import Math.NumberTheory.Primes.Testing+import Math.NumberTheory.TestUtils++-- | https://en.wikipedia.org/wiki/Prime-counting_function#Table_of_.CF.80.28x.29.2C_x_.2F_ln_x.2C_and_li.28x.29+table :: [(Integer, Integer)]+table =+  [ (10^1,   4)+  , (10^2,   25)+  , (10^3,   168)+  , (10^4,   1229)+  , (10^5,   9592)+  , (10^6,   78498)+  , (10^7,   664579)+  , (10^8,   5761455)+  , (10^9,   50847534)+  , (10^10,  455052511)+  , (10^11,  4118054813)+  , (10^12,  37607912018)+  , (10^13,  346065536839)+  -- , (10^14,  3204941750802)+  -- , (10^15,  29844570422669)+  -- , (10^16,  279238341033925)+  -- , (10^17,  2623557157654233)+  -- , (10^18,  24739954287740860)+  -- , (10^19,  234057667276344607)+  -- , (10^20,  2220819602560918840)+  ]++-- | Check that values of 'primeCount' are non-negative.+primeCountProperty1 :: Integer -> Bool+primeCountProperty1 n = n > primeCountMaxArg+  || n >  0 && primeCount n >= 0+  || n <= 0 && primeCount n == 0++-- | Check that 'primeCount' is monotonically increasing function.+primeCountProperty2 :: Positive Integer -> Positive Integer -> Bool+primeCountProperty2 (Positive n1) (Positive n2)+  =  n1 > primeCountMaxArg+  || n2 > primeCountMaxArg+  || n1 <= n2 && p1 <= p2+  || n1 >  n2 && p1 >= p2+  where+    p1 = primeCount n1+    p2 = primeCount n2++-- | Check that 'primeCount' is strictly increasing iff an argument is prime.+primeCountProperty3 :: Positive Integer -> Bool+primeCountProperty3 (Positive n)+  =  isPrime n && primeCount (n - 1) + 1 == primeCount n+  || not (isPrime n) && primeCount (n - 1) == primeCount n++-- | Check tabulated values.+primeCountSpecialCases :: [Assertion]+primeCountSpecialCases = map a table+  where+  a (n, m) = assertEqual "primeCount" m (primeCount n)+++-- | Check that values of 'nthPrime' are positive.+nthPrimeProperty1 :: Positive Integer -> Bool+nthPrimeProperty1 (Positive n) = n > nthPrimeMaxArg+  || nthPrime n > 0++-- | Check that 'nthPrime' is monotonically increasing function.+nthPrimeProperty2 :: Positive Integer -> Positive Integer -> Bool+nthPrimeProperty2 (Positive n1) (Positive n2)+  =  n1 > nthPrimeMaxArg+  || n2 > nthPrimeMaxArg+  || n1 <= n2 && p1 <= p2+  || n1 >  n2 && p1 >= p2+  where+    p1 = nthPrime n1+    p2 = nthPrime n2++-- | Check that values of 'nthPrime' are prime.+nthPrimeProperty3 :: Positive Integer -> Bool+nthPrimeProperty3 (Positive n) = isPrime $ nthPrime n++-- | Check tabulated values.+nthPrimeSpecialCases :: [Assertion]+nthPrimeSpecialCases = map a table+  where+  a (n, m) = assertBool "nthPrime" $ n > nthPrime m+++-- | Check that values of 'approxPrimeCount' are non-negative.+approxPrimeCountProperty1 :: Integral a => AnySign a -> Bool+approxPrimeCountProperty1 (AnySign a) = approxPrimeCount a >= 0++-- | Check that 'approxPrimeCount' is consistent with 'approxPrimeCountOverestimateLimit'.+approxPrimeCountProperty2 :: Integral a => Positive a -> Bool+approxPrimeCountProperty2 (Positive a) = a >= approxPrimeCountOverestimateLimit+  || toInteger (approxPrimeCount a) >= primeCount (toInteger a)+++-- | Check that values of 'nthPrimeApprox' are positive.+nthPrimeApproxProperty1 :: AnySign Integer -> Bool+nthPrimeApproxProperty1 (AnySign a) = nthPrimeApprox a > 0++-- | Check that 'nthPrimeApprox' is consistent with 'nthPrimeApproxUnderestimateLimit'.+nthPrimeApproxProperty2 :: Positive Integer -> Bool+nthPrimeApproxProperty2 (Positive a) = a >= nthPrimeApproxUnderestimateLimit+  || toInteger (nthPrimeApprox a) <= nthPrime (toInteger a)+++testSuite :: TestTree+testSuite = testGroup "Counting"+  [ testGroup "primeCount"+    ( testSmallAndQuick "non-negative"        primeCountProperty1+    : testSmallAndQuick "monotonic"           primeCountProperty2+    : testSmallAndQuick "increases on primes" primeCountProperty3+    : zipWith (\i a -> testCase ("special case " ++ show i) a) [1..] primeCountSpecialCases+    )+  , testGroup "nthPrime"+    ( testSmallAndQuick "positive"  nthPrimeProperty1+    : testSmallAndQuick "monotonic" nthPrimeProperty2+    : testSmallAndQuick "is prime"  nthPrimeProperty3+    : zipWith (\i a -> testCase ("special case " ++ show i) a) [1..] nthPrimeSpecialCases+    )+  , testGroup "approxPrimeCount"+    [ testIntegralProperty "non-negative"             approxPrimeCountProperty1+    , testIntegralProperty "overestimates primeCount" approxPrimeCountProperty2+    ]+  , testGroup "nthPrimeApprox"+    [ testSmallAndQuick "positive"                nthPrimeApproxProperty1+    , testSmallAndQuick "underestimates nthPrime" nthPrimeApproxProperty2+    ]+  ]

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

@@ -0,0 +1,69 @@+-- |+-- Module:      Math.NumberTheory.Primes.SieveTests+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.Primes.Sieve+--++{-# LANGUAGE CPP #-}++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.Primes.SieveTests+  ( testSuite+  ) where++import Prelude hiding (words)++import Test.Tasty+import Test.Tasty.HUnit++import Math.NumberTheory.Primes.Sieve+import qualified Math.NumberTheory.Primes.Heap as H+import Math.NumberTheory.TestUtils++-- | Check that both 'primes' produce the same.+primesProperty1 :: Assertion+primesProperty1 = do+  assertEqual "Sieve == Heap" (trim primes) (trim H.primes)+  where+    trim = take 100000++-- | Check that both 'sieveFrom' produce the same.+sieveFromProperty1 :: AnySign Integer -> Bool+sieveFromProperty1 (AnySign lowBound)+  = trim (sieveFrom lowBound) == trim (H.sieveFrom lowBound)+  where+    trim = take 1000++-- | Check that 'primeList' from 'primeSieve' matches truncated 'primes'.+primeSieveProperty1 :: AnySign Integer -> Bool+primeSieveProperty1 (AnySign highBound)+  = primeList (primeSieve highBound) == takeWhile (<= (highBound `max` 7)) primes++-- | Check that 'primeList' from 'psieveList' matches 'primes'.+psieveListProperty1 :: Assertion+psieveListProperty1 = do+  assertEqual "primes == primeList . psieveList" (trim primes) (trim $ concatMap primeList psieveList)+  where+    trim = take 100000++-- | Check that 'primeList' from 'psieveFrom' matches 'sieveFrom'.+psieveFromProperty1 :: AnySign Integer -> Bool+psieveFromProperty1 (AnySign lowBound)+  = trim (sieveFrom lowBound) == trim (filter (>= lowBound) (concatMap primeList $ psieveFrom lowBound))+  where+    trim = take 1000+++testSuite :: TestTree+testSuite = testGroup "Sieve"+  [ testCase          "primes"     primesProperty1+  , testSmallAndQuick "sieveFrom"  sieveFromProperty1+  , testSmallAndQuick "primeSieve" primeSieveProperty1+  , testCase          "psieveList" psieveListProperty1+  , testSmallAndQuick "psieveFrom" psieveFromProperty1+  ]

@@ -0,0 +1,40 @@+-- |+-- Module:      Math.NumberTheory.PrimesTests+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.Primes+--++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.PrimesTests+  ( testSuite+  ) where++import Test.Tasty+import Test.Tasty.HUnit++import Math.NumberTheory.Primes+import Math.NumberTheory.TestUtils++primesSumWonk :: Int -> Int+primesSumWonk upto = sum . takeWhile (< upto) . map fromInteger . primeList $ primeSieve (toInteger upto)++primesSum :: Int -> Int+primesSum upto = sum . takeWhile (< upto) . map fromInteger $ primes++primesSumProperty :: NonNegative Int -> Bool+primesSumProperty (NonNegative n) = primesSumWonk n == primesSum n+++sieveFactorSpecialCase1 :: Assertion+sieveFactorSpecialCase1 = assertEqual "sieveFactor" [(29, 1), (73, 1)] $ sieveFactor (factorSieve 2048) (29*73)++testSuite :: TestTree+testSuite = testGroup "Primes"+  [ testSmallAndQuick "primesSum"   primesSumProperty+  , testCase          "sieveFactor" sieveFactorSpecialCase1+  ]

@@ -0,0 +1,226 @@+-- |+-- Module:      Math.NumberTheory.TestUtils+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Utils to test Math.NumberTheory+--++{-# LANGUAGE ConstraintKinds            #-}+{-# LANGUAGE CPP                        #-}+{-# LANGUAGE DataKinds                  #-}+{-# LANGUAGE DeriveFoldable             #-}+{-# LANGUAGE DeriveFunctor              #-}+{-# LANGUAGE DeriveTraversable          #-}+{-# LANGUAGE FlexibleContexts           #-}+{-# LANGUAGE FlexibleInstances          #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE KindSignatures             #-}+{-# LANGUAGE MultiParamTypeClasses      #-}+{-# LANGUAGE RankNTypes                 #-}+{-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE TypeFamilies               #-}+{-# LANGUAGE TypeOperators              #-}+{-# LANGUAGE UndecidableInstances       #-}++#if __GLASGOW_HASKELL__ >= 800+{-# LANGUAGE UndecidableSuperClasses    #-}++{-# OPTIONS_GHC -fconstraint-solver-iterations=0 #-}+#endif++{-# OPTIONS_GHC -fno-warn-orphans #-}+{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.TestUtils+  ( module Math.NumberTheory.TestUtils+  , module Test.SmallCheck.Series+  , Large(..)+  ) where++import Test.Tasty+import Test.Tasty.SmallCheck as SC+import Test.Tasty.QuickCheck as QC hiding (Positive, NonNegative, generate, getNonNegative)++import Test.SmallCheck.Series (Positive(..), NonNegative(..), Serial(..), Series, generate)++import Control.Applicative+import Data.Bits+#if MIN_VERSION_base(4,8,0)+#else+import Data.Foldable (Foldable)+import Data.Traversable (Traversable)+import Data.Word+#endif+import GHC.Exts++import Math.NumberTheory.Primes++newtype AnySign a = AnySign { getAnySign :: a }+  deriving (Eq, Ord, Read, Show, Num, Enum, Bounded, Integral, Real, Functor, Foldable, Traversable, Arbitrary)++instance (Monad m, Serial m a) => Serial m (AnySign a) where+  series = AnySign <$> series++instance (Num a, Ord a, Arbitrary a) => Arbitrary (Positive a) where+  arbitrary = Positive <$> (arbitrary `suchThat` (> 0))+  shrink (Positive x) = Positive <$> filter (> 0) (shrink x)++instance (Num a, Ord a, Arbitrary a) => Arbitrary (NonNegative a) where+  arbitrary = NonNegative <$> (arbitrary `suchThat` (>= 0))+  shrink (NonNegative x) = NonNegative <$> filter (>= 0) (shrink x)++instance (Num a, Bounded a) => Bounded (Positive a) where+  minBound = Positive 1+  maxBound = Positive (maxBound :: a)++instance (Num a, Bounded a) => Bounded (NonNegative a) where+  minBound = NonNegative 0+  maxBound = NonNegative (maxBound :: a)++newtype Huge a = Huge { getHuge :: a }+  deriving (Eq, Ord, Enum, Bounded, Show, Num, Real, Integral)++instance (Num a, Arbitrary a) => Arbitrary (Huge a) where+  arbitrary = do+    Positive l <- arbitrary+    ds <- vector l+    return $ Huge $ foldl1 (\acc n -> acc * 2^63 + n) ds++newtype Power a = Power { getPower :: a }+  deriving (Eq, Ord, Enum, Bounded, Show, Num, Real, Integral)++instance (Monad m, Num a, Ord a, Serial m a) => Serial m (Power a) where+  series = Power <$> series `suchThatSerial` (> 0)++instance (Num a, Ord a, Integral a, Arbitrary a) => Arbitrary (Power a) where+  arbitrary = Power <$> (getSmall <$> arbitrary) `suchThat` (> 0)+  shrink (Power x) = Power <$> filter (> 0) (shrink x)++newtype Prime = Prime { getPrime :: Integer }+  deriving (Eq, Ord, Show)++instance Arbitrary Prime where+  arbitrary = Prime <$> arbitrary `suchThat` (\p -> p > 0 && isPrime p)++instance Monad m => Serial m Prime where+  series = Prime <$> series `suchThatSerial` (\p -> p > 0 && isPrime p)++instance Monad m => Serial m Word where+  series =+    generate (\d -> if d >= 0 then pure 0 else empty) <|> nats+    where+      nats = generate $ \d -> if d > 0 then [1 .. fromInteger (toInteger d)] else empty++suchThatSerial :: Series m a -> (a -> Bool) -> Series m a+suchThatSerial s p = s >>= \x -> if p x then pure x else empty+++-- https://www.cs.ox.ac.uk/projects/utgp/school/andres.pdf, p. 21+-- :k Compose = (k1 -> Constraint) -> (k2 -> k1) -> (k2 -> Constraint)+class    (f (g x)) => (f `Compose` g) x+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]+  , Matrix '[Bounded, Integral] wrapper '[Int, Word]+  , Num (wrapper Integer)+  )+++testIntegralProperty+  :: forall wrapper bool. (TestableIntegral wrapper, SC.Testable IO bool, QC.Testable bool)+  => String -> (forall a. (Integral a, Bits a) => wrapper a -> bool) -> TestTree+testIntegralProperty name f = testGroup name+  [ SC.testProperty "smallcheck Int"     (f :: wrapper Int     -> bool)+  , SC.testProperty "smallcheck Word"    (f :: wrapper Word    -> bool)+  , SC.testProperty "smallcheck Integer" (f :: wrapper Integer -> bool)+  , QC.testProperty "quickcheck Int"     (f :: wrapper Int     -> bool)+  , QC.testProperty "quickcheck Word"    (f :: wrapper Word    -> bool)+  , QC.testProperty "quickcheck Integer" (f :: wrapper Integer -> bool)+  , QC.testProperty "quickcheck Large Int"     ((f :: wrapper Int     -> bool) . getLarge)+  , QC.testProperty "quickcheck Large Word"    ((f :: wrapper Word    -> bool) . getLarge)+  , QC.testProperty "quickcheck Huge  Integer" ((f :: wrapper Integer -> bool) . getHuge)+  ]++testSameIntegralProperty+  :: forall wrapper1 wrapper2 bool. (TestableIntegral wrapper1, TestableIntegral wrapper2, SC.Testable IO bool, QC.Testable bool)+  => String -> (forall a. (Integral a, Bits a) => wrapper1 a -> wrapper2 a -> bool) -> TestTree+testSameIntegralProperty name f = testGroup name+  [ SC.testProperty "smallcheck Int"     (f :: wrapper1 Int     -> wrapper2 Int     -> bool)+  , SC.testProperty "smallcheck Word"    (f :: wrapper1 Word    -> wrapper2 Word    -> bool)+  , SC.testProperty "smallcheck Integer" (f :: wrapper1 Integer -> wrapper2 Integer -> bool)+  , QC.testProperty "quickcheck Int"     (f :: wrapper1 Int     -> wrapper2 Int     -> bool)+  , QC.testProperty "quickcheck Word"    (f :: wrapper1 Word    -> wrapper2 Word    -> bool)+  , QC.testProperty "quickcheck Integer" (f :: wrapper1 Integer -> wrapper2 Integer -> bool)+  , QC.testProperty "quickcheck Large Int"     (\(Large a) (Large b) -> (f :: wrapper1 Int     -> wrapper2 Int     -> bool) a b)+  , QC.testProperty "quickcheck Large Word"    (\(Large a) (Large b) -> (f :: wrapper1 Word    -> wrapper2 Word    -> bool) a b)+  , QC.testProperty "quickcheck Huge  Integer" (\(Huge  a) (Huge  b) -> (f :: wrapper1 Integer -> wrapper2 Integer -> bool) a b)+  ]++testIntegral2Property+  :: forall wrapper1 wrapper2 bool. (TestableIntegral wrapper1, TestableIntegral wrapper2, SC.Testable IO bool, QC.Testable bool)+  => String -> (forall a1 a2. (Integral a1, Integral a2, Bits a1, Bits a2) => wrapper1 a1 -> wrapper2 a2 -> bool) -> TestTree+testIntegral2Property name f = testGroup name+  [ SC.testProperty "smallcheck Int Int"         (f :: wrapper1 Int     -> wrapper2 Int     -> bool)+  , SC.testProperty "smallcheck Int Word"        (f :: wrapper1 Int     -> wrapper2 Word    -> bool)+  , SC.testProperty "smallcheck Int Integer"     (f :: wrapper1 Int     -> wrapper2 Integer -> bool)+  , SC.testProperty "smallcheck Word Int"        (f :: wrapper1 Word    -> wrapper2 Int     -> bool)+  , SC.testProperty "smallcheck Word Word"       (f :: wrapper1 Word    -> wrapper2 Word    -> bool)+  , SC.testProperty "smallcheck Word Integer"    (f :: wrapper1 Word    -> wrapper2 Integer -> bool)+  , SC.testProperty "smallcheck Integer Int"     (f :: wrapper1 Integer -> wrapper2 Int     -> bool)+  , SC.testProperty "smallcheck Integer Word"    (f :: wrapper1 Integer -> wrapper2 Word    -> bool)+  , SC.testProperty "smallcheck Integer Integer" (f :: wrapper1 Integer -> wrapper2 Integer -> bool)++  , QC.testProperty "quickcheck Int Int"         (f :: wrapper1 Int     -> wrapper2 Int     -> bool)+  , QC.testProperty "quickcheck Int Word"        (f :: wrapper1 Int     -> wrapper2 Word    -> bool)+  , QC.testProperty "quickcheck Int Integer"     (f :: wrapper1 Int     -> wrapper2 Integer -> bool)+  , QC.testProperty "quickcheck Word Int"        (f :: wrapper1 Word    -> wrapper2 Int     -> bool)+  , QC.testProperty "quickcheck Word Word"       (f :: wrapper1 Word    -> wrapper2 Word    -> bool)+  , QC.testProperty "quickcheck Word Integer"    (f :: wrapper1 Word    -> wrapper2 Integer -> bool)+  , QC.testProperty "quickcheck Integer Int"     (f :: wrapper1 Integer -> wrapper2 Int     -> bool)+  , QC.testProperty "quickcheck Integer Word"    (f :: wrapper1 Integer -> wrapper2 Word    -> bool)+  , QC.testProperty "quickcheck Integer Integer" (f :: wrapper1 Integer -> wrapper2 Integer -> bool)++  , QC.testProperty "quickcheck Large Int Int"         ((f :: wrapper1 Int     -> wrapper2 Int     -> bool) . getLarge)+  , QC.testProperty "quickcheck Large Int Word"        ((f :: wrapper1 Int     -> wrapper2 Word    -> bool) . getLarge)+  , QC.testProperty "quickcheck Large Int Integer"     ((f :: wrapper1 Int     -> wrapper2 Integer -> bool) . getLarge)+  , QC.testProperty "quickcheck Large Word Int"        ((f :: wrapper1 Word    -> wrapper2 Int     -> bool) . getLarge)+  , QC.testProperty "quickcheck Large Word Word"       ((f :: wrapper1 Word    -> wrapper2 Word    -> bool) . getLarge)+  , QC.testProperty "quickcheck Large Word Integer"    ((f :: wrapper1 Word    -> wrapper2 Integer -> bool) . getLarge)+  , QC.testProperty "quickcheck Huge  Integer Int"     ((f :: wrapper1 Integer -> wrapper2 Int     -> bool) . getHuge)+  , QC.testProperty "quickcheck Huge  Integer Word"    ((f :: wrapper1 Integer -> wrapper2 Word    -> bool) . getHuge)+  , QC.testProperty "quickcheck Huge  Integer Integer" ((f :: wrapper1 Integer -> wrapper2 Integer -> bool) . getHuge)+  ]++testSmallAndQuick+  :: SC.Testable IO a+  => QC.Testable a+  => String -> a -> TestTree+testSmallAndQuick name f = testGroup name+  [ SC.testProperty "smallcheck" f+  , QC.testProperty "quickcheck" f+  ]

@@ -1,1 +0,0 @@-{-# OPTIONS_GHC -F -pgmF hspec-discover #-}

@@ -0,0 +1,66 @@+import Test.Tasty++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.ModuliTests as Moduli++import qualified Math.NumberTheory.MoebiusInversionTests as MoebiusInversion+import qualified Math.NumberTheory.MoebiusInversion.IntTests as MoebiusInversionInt++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.HeapTests as Heap+import qualified Math.NumberTheory.Primes.SieveTests as Sieve++import qualified Math.NumberTheory.GaussianIntegersTests as Gaussian++main :: IO ()+main = defaultMain tests++tests :: TestTree+tests = testGroup "All"+  [ testGroup "Powers"+    [ 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 "Moduli"+    [ Moduli.testSuite+    ]+  , testGroup "MoebiusInversion"+    [ MoebiusInversion.testSuite+    , MoebiusInversionInt.testSuite+    ]+  , testGroup "Primes"+    [ Primes.testSuite+    , Counting.testSuite+    , Heap.testSuite+    , Sieve.testSuite+    ]+  , testGroup "Gaussian"+    [ Gaussian.testSuite+    ]+  ]