arithmoi 0.4.2.0 → 0.4.3.0

raw patch · 40 files changed

+1656/−457 lines, 40 filesdep +natsdep +semigroupsdep +transformersdep ~QuickCheckdep ~arraydep ~basePVP ok

version bump matches the API change (PVP)

Dependencies added: nats, semigroups, transformers, transformers-compat

Dependency ranges changed: QuickCheck, array, base, containers, random

API changes (from Hackage documentation)

Files

• Changes +15/−1
• Math/NumberTheory/ArithmeticFunctions.hs +20/−0
• Math/NumberTheory/ArithmeticFunctions/Class.hs +102/−0
• Math/NumberTheory/ArithmeticFunctions/Standard.hs +299/−0
• Math/NumberTheory/GCD/LowLevel.hs +2/−5
• Math/NumberTheory/GaussianIntegers.hs +2/−1
• Math/NumberTheory/Logarithms.hs +2/−5
• Math/NumberTheory/Logarithms/Internal.hs +0/−155
• Math/NumberTheory/Powers/Cubes.hs +1/−4
• Math/NumberTheory/Powers/Fourth.hs +1/−4
• Math/NumberTheory/Powers/General.hs +1/−4
• Math/NumberTheory/Powers/Integer.hs +1/−4
• Math/NumberTheory/Powers/Squares.hs +46/−76
• Math/NumberTheory/Powers/Squares/Internal.hs +144/−0
• Math/NumberTheory/Primes/Counting/Impl.hs +0/−5
• Math/NumberTheory/Primes/Factorisation.hs +5/−0
• Math/NumberTheory/Primes/Factorisation/Montgomery.hs +2/−4
• Math/NumberTheory/Primes/Factorisation/Utils.hs +2/−0
• Math/NumberTheory/Primes/Heap.hs +1/−28
• Math/NumberTheory/Primes/Sieve/Eratosthenes.hs +0/−5
• Math/NumberTheory/Primes/Sieve/Misc.hs +31/−5
• Math/NumberTheory/Primes/Testing/Probabilistic.hs +0/−3
• Math/NumberTheory/UniqueFactorisation.hs +131/−0
• Math/NumberTheory/Unsafe.hs +0/−3
• Math/NumberTheory/Utils.hs +0/−20
• arithmoi.cabal +43/−25
• benchmark/Bench.hs +11/−1
• benchmark/Math/NumberTheory/ArithmeticFunctionsBench.hs +34/−0
• benchmark/Math/NumberTheory/PowersBench.hs +30/−0
• test-suite/Math/NumberTheory/ArithmeticFunctionsTests.hs +279/−0
• test-suite/Math/NumberTheory/ModuliTests.hs +23/−19
• test-suite/Math/NumberTheory/MoebiusInversion/IntTests.hs +1/−1
• test-suite/Math/NumberTheory/MoebiusInversionTests.hs +1/−1
• test-suite/Math/NumberTheory/Powers/SquaresTests.hs +14/−4
• test-suite/Math/NumberTheory/PrimesTests.hs +8/−2
• test-suite/Math/NumberTheory/TestUtils.hs +41/−72
• test-suite/Math/NumberTheory/TestUtils/Compose.hs +47/−0
• test-suite/Math/NumberTheory/TestUtils/Wrappers.hs +255/−0
• test-suite/Math/NumberTheory/UniqueFactorisationTests.hs +52/−0
• test-suite/Test.hs +9/−0

@@ -1,5 +1,19 @@+0.4.3.0:+    This release supports GHC 7.6, 7.8, 7.10 and 8.0.++    Add Math.NumberTheory.ArithmeticFunctions with brand-new machinery+    for arithmetic functions: divisors, tau, sigma, totient, jordan,+    moebius, liouville, smallOmega, bigOmega, carmichael, expMangoldt (#30).+    Old implementations (exposed via Math.NumberTheory.Primes.Factorisation+    and Math.NumberTheory.Powers.Integer) are deprecated and will be removed+    in the next major release.++    Add Karatsuba sqrt algorithm, improving performance on large integers (#6).++    Fix incorrect indexing of FactorSieve (#35).+ 0.4.2.0:-    This release supports GHC 7.6, 7.8 and 8.0.+    This release supports GHC 7.6, 7.8, 7.10 and 8.0.      Add new cabal flag check-bounds, which replaces all unsafe array functions with safe ones. 

@@ -0,0 +1,20 @@+-- |+-- Module:      Math.NumberTheory.ArithmeticFunctions+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- This module provides an interface for defining and manipulating+-- arithmetic functions. It also defines several most widespreaded+-- arithmetic functions.+--++module Math.NumberTheory.ArithmeticFunctions+  ( module Math.NumberTheory.ArithmeticFunctions.Class+  , module Math.NumberTheory.ArithmeticFunctions.Standard+  ) where++import Math.NumberTheory.ArithmeticFunctions.Class+import Math.NumberTheory.ArithmeticFunctions.Standard

@@ -0,0 +1,102 @@+-- |+-- Module:      Math.NumberTheory.ArithmeticFunctions.Class+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Generic type for arithmetic functions over arbitrary unique+-- factorisation domains.+--++{-# LANGUAGE CPP                 #-}+{-# LANGUAGE GADTs               #-}++{-# OPTIONS_HADDOCK hide #-}++module Math.NumberTheory.ArithmeticFunctions.Class+  ( ArithmeticFunction(..)+  , runFunction+  ) where++import Control.Applicative+import Data.Semigroup++#if MIN_VERSION_base(4,8,0)+#else+import Data.Word+#endif++import Math.NumberTheory.UniqueFactorisation++-- | A typical arithmetic function operates on the canonical factorisation of+-- a number into prime's powers and consists of two rules. The first one+-- determines the values of the function on the powers of primes. The second+-- one determines how to combine these values into final result.+--+-- In the following definition the first argument is the function on prime's+-- powers, the monoid instance determines a rule of combination (typically+-- 'Product' or 'Sum'), and the second argument is convenient for unwrapping+-- (typically, 'getProduct' or 'getSum').+data ArithmeticFunction n a where+  ArithmeticFunction+    :: Monoid m+    => (Prime n -> Word -> m)+    -> (m -> a)+    -> ArithmeticFunction n a++-- | Convert to function. The value on 0 is undefined.+runFunction :: UniqueFactorisation n => ArithmeticFunction n a -> n -> a+runFunction (ArithmeticFunction f g)+  = g+  . mconcat+  . map (uncurry f)+  . factorise++instance Functor (ArithmeticFunction n) where+  fmap f (ArithmeticFunction g h) = ArithmeticFunction g (f . h)++instance Applicative (ArithmeticFunction n) where+  pure x+    = ArithmeticFunction (\_ _ -> ()) (const x)+  (ArithmeticFunction f1 g1) <*> (ArithmeticFunction f2 g2)+    = ArithmeticFunction (\p k -> (f1 p k, f2 p k)) (\(a1, a2) -> g1 a1 (g2 a2))++instance Semigroup a => Semigroup (ArithmeticFunction n a) where+  (<>) = liftA2 (<>)++instance Monoid a => Monoid (ArithmeticFunction n a) where+  mempty  = pure mempty+  mappend = liftA2 mappend++-- | Factorisation is expensive, so it is better to avoid doing it twice.+-- Write 'runFunction (f + g) n' instead of 'runFunction f n + runFunction g n'.+instance Num a => Num (ArithmeticFunction n a) where+  fromInteger = pure . fromInteger+  negate = fmap negate+  signum = fmap signum+  abs    = fmap abs+  (+) = liftA2 (+)+  (-) = liftA2 (-)+  (*) = liftA2 (*)++instance Fractional a => Fractional (ArithmeticFunction n a) where+  fromRational = pure . fromRational+  recip = fmap recip+  (/) = liftA2 (/)++instance Floating a => Floating (ArithmeticFunction n a) where+  pi    = pure pi+  exp   = fmap exp+  log   = fmap log+  sin   = fmap sin+  cos   = fmap cos+  asin  = fmap asin+  acos  = fmap acos+  atan  = fmap atan+  sinh  = fmap sinh+  cosh  = fmap cosh+  asinh = fmap asinh+  acosh = fmap acosh+  atanh = fmap atanh

@@ -0,0 +1,299 @@+-- |+-- Module:      Math.NumberTheory.ArithmeticFunctions.Standard+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Textbook arithmetic functions.+--++{-# LANGUAGE CPP                 #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeFamilies        #-}+{-# LANGUAGE ViewPatterns        #-}++{-# OPTIONS_HADDOCK hide #-}++module Math.NumberTheory.ArithmeticFunctions.Standard+  ( -- * Multiplicative functions+    multiplicative+  , divisors, divisorsA+  , divisorsSmall, divisorsSmallA+  , tau, tauA+  , sigma, sigmaA+  , totient, totientA+  , jordan, jordanA+  , moebius, moebiusA+  , liouville, liouvilleA+    -- * Additive functions+  , additive+  , smallOmega, smallOmegaA+  , bigOmega, bigOmegaA+    -- * Misc+  , carmichael, carmichaelA+  , expMangoldt, expMangoldtA+  ) where++import Data.IntSet (IntSet)+import qualified Data.IntSet as IS+import Data.Set (Set)+import qualified Data.Set as S+import Data.Semigroup++import Math.NumberTheory.ArithmeticFunctions.Class+import Math.NumberTheory.UniqueFactorisation++import Numeric.Natural++#if MIN_VERSION_base(4,8,0)+#else+import Data.Foldable+import Data.Word+#endif++#if MIN_VERSION_base(4,7,0)+import Data.Coerce+#else+import Unsafe.Coerce++coerce :: a -> b+coerce = unsafeCoerce+#endif++wordToInt :: Word -> Int+wordToInt = fromIntegral++-- | Create a multiplicative function from the function on prime's powers. See examples below.+multiplicative :: Num a => (Prime n -> Word -> a) -> ArithmeticFunction n a+multiplicative f = ArithmeticFunction ((Product .) . f) getProduct++divisors :: (UniqueFactorisation n, Num n, Ord n) => n -> Set n+divisors = runFunction divisorsA+{-# SPECIALIZE divisors :: Natural -> Set Natural #-}+{-# SPECIALIZE divisors :: Integer -> Set Integer #-}++-- | The set of all (positive) divisors of an argument.+divisorsA :: forall n. (UniqueFactorisation n, Num n, Ord n) => ArithmeticFunction n (Set n)+divisorsA = ArithmeticFunction (\((unPrime :: Prime n -> n) -> p) k -> SetProduct $ divisorsHelper p k) (S.insert 1 . getSetProduct)++divisorsHelper :: Num n => n -> Word -> Set n+divisorsHelper _ 0 = S.empty+divisorsHelper p 1 = S.singleton p+divisorsHelper p a = S.fromDistinctAscList $ p : p * p : map (p ^) [3 .. wordToInt a]+{-# INLINE divisorsHelper #-}++divisorsSmall :: (UniqueFactorisation n, Prime n ~ Prime Int) => n -> IntSet+divisorsSmall = runFunction divisorsSmallA++-- | Same as 'divisors', but with better performance on cost of type restriction.+divisorsSmallA :: forall n. (Prime n ~ Prime Int) => ArithmeticFunction n IntSet+divisorsSmallA = ArithmeticFunction (\p k -> IntSetProduct $ divisorsHelperSmall (unPrime p) k) (IS.insert 1 . getIntSetProduct)++divisorsHelperSmall :: Int -> Word -> IntSet+divisorsHelperSmall _ 0 = IS.empty+divisorsHelperSmall p 1 = IS.singleton p+divisorsHelperSmall p a = IS.fromDistinctAscList $ p : p * p : map (p ^) [3 .. wordToInt a]+{-# INLINE divisorsHelperSmall #-}++tau :: (UniqueFactorisation n, Num a) => n -> a+tau = runFunction tauA++-- | The number of (positive) divisors of an argument.+--+-- > tauA = multiplicative (\_ k -> k + 1)+tauA :: Num a => ArithmeticFunction n a+tauA = multiplicative $ const (fromIntegral . succ)++sigma :: (UniqueFactorisation n, Integral n) => Word -> n -> n+sigma = runFunction . sigmaA++-- | The sum of the @k@-th powers of (positive) divisors of an argument.+--+-- > sigmaA = multiplicative (\p k -> sum $ map (p ^) [0..k])+-- > sigmaA 0 = tauA+sigmaA :: forall n. (UniqueFactorisation n, Integral n) => Word -> ArithmeticFunction n n+sigmaA 0 = tauA+sigmaA 1 = multiplicative $ \((unPrime :: Prime n -> n) -> p) -> sigmaHelper p+sigmaA a = multiplicative $ \((unPrime :: Prime n -> n) -> p) -> sigmaHelper (p ^ wordToInt a)++sigmaHelper :: Integral n => n -> Word -> n+sigmaHelper pa 1 = pa + 1+sigmaHelper pa 2 = pa * pa + pa + 1+sigmaHelper pa k = (pa ^ wordToInt (k + 1) - 1) `quot` (pa - 1)+{-# INLINE sigmaHelper #-}++totient :: (UniqueFactorisation n, Integral n) => n -> n+totient = runFunction totientA++-- | Calculates the totient of a positive number @n@, i.e.+--   the number of @k@ with @1 <= k <= n@ and @'gcd' n k == 1@,+--   in other words, the order of the group of units in @&#8484;/(n)@.+totientA :: forall n. (UniqueFactorisation n, Integral n) => ArithmeticFunction n n+totientA = multiplicative $ \((unPrime :: Prime n -> n) -> p) -> jordanHelper p++jordan :: (UniqueFactorisation n, Integral n) => Word -> n -> n+jordan = runFunction . jordanA++-- | Calculates the k-th Jordan function of an argument.+--+-- > jordanA 1 = totientA+jordanA :: forall n. (UniqueFactorisation n, Integral n) => Word -> ArithmeticFunction n n+jordanA 0 = multiplicative $ \_ _ -> 0+jordanA 1 = totientA+jordanA a = multiplicative $ \((unPrime :: Prime n -> n) -> p) -> jordanHelper (p ^ wordToInt a)++jordanHelper :: Integral n => n -> Word -> n+jordanHelper pa 1 = pa - 1+jordanHelper pa 2 = (pa - 1) * pa+jordanHelper pa k = (pa - 1) * pa ^ wordToInt (k - 1)+{-# INLINE jordanHelper #-}++moebius :: (UniqueFactorisation n, Num a) => n -> a+moebius = runFunction moebiusA++-- | Calculates the Moebius function of an argument.+moebiusA :: Num a => ArithmeticFunction n a+moebiusA = ArithmeticFunction (const f) runMoebius+  where+    f 1 = MoebiusN+    f 0 = MoebiusP+    f _ = MoebiusZ++liouville :: (UniqueFactorisation n, Num a) => n -> a+liouville = runFunction liouvilleA++-- | Calculates the Liouville function of an argument.+liouvilleA :: Num a => ArithmeticFunction n a+liouvilleA = ArithmeticFunction (const $ Xor . odd) runXor++carmichael :: (UniqueFactorisation n, Integral n) => n -> n+carmichael = runFunction carmichaelA+{- The specializations reflects available specializations of lcm. -}+{-# SPECIALIZE carmichael :: Int -> Int #-}+{-# SPECIALIZE carmichael :: Integer -> Integer #-}++-- | Calculates the Carmichael function for a positive integer, that is,+--   the (smallest) exponent of the group of units in @&#8484;/(n)@.+carmichaelA :: forall n. (UniqueFactorisation n, Integral n) => ArithmeticFunction n n+carmichaelA = ArithmeticFunction (\((unPrime :: Prime n -> n) -> p) k -> LCM $ f p k) getLCM+  where+    f 2 1 = 1+    f 2 2 = 2+    f 2 k = 2 ^ wordToInt (k - 2)+    f p 1 = p - 1+    f p 2 = (p - 1) * p+    f p k = (p - 1) * p ^ wordToInt (k - 1)++-- | Create an additive function from the function on prime's powers. See examples below.+additive :: Num a => (Prime n -> Word -> a) -> ArithmeticFunction n a+additive f = ArithmeticFunction ((Sum .) . f) getSum++smallOmega :: (UniqueFactorisation n, Num a) => n -> a+smallOmega = runFunction smallOmegaA++-- | Number of distinct prime factors.+--+-- > smallOmegaA = additive (\_ _ -> 1)+smallOmegaA :: Num a => ArithmeticFunction n a+smallOmegaA = additive (\_ _ -> 1)++bigOmega :: UniqueFactorisation n => n -> Word+bigOmega = runFunction bigOmegaA++-- | Number of prime factors, counted with multiplicity.+--+-- > bigOmegaA = additive (\_ k -> k)+bigOmegaA :: ArithmeticFunction n Word+bigOmegaA = additive $ const id++expMangoldt :: (UniqueFactorisation n, Num n) => n -> n+expMangoldt = runFunction expMangoldtA++-- | The exponent of von Mangoldt function. Use @log expMangoldtA@ to recover von Mangoldt function itself.+expMangoldtA :: forall n. (UniqueFactorisation n, Num n) => ArithmeticFunction n n+expMangoldtA = ArithmeticFunction (\((unPrime :: Prime n -> n) -> p) _ -> MangoldtOne p) runMangoldt++data Moebius+  = MoebiusZ+  | MoebiusP+  | MoebiusN++runMoebius :: Num a => Moebius -> a+runMoebius m = case m of+  MoebiusZ ->  0+  MoebiusP ->  1+  MoebiusN -> -1++instance Semigroup Moebius where+  MoebiusZ <> _ = MoebiusZ+  _ <> MoebiusZ = MoebiusZ+  MoebiusP <> a = a+  a <> MoebiusP = a+  _ <> _ = MoebiusP++instance Monoid Moebius where+  mempty = MoebiusP+  mappend = (<>)++data Mangoldt a+  = MangoldtZero+  | MangoldtOne a+  | MangoldtMany++runMangoldt :: Num a => Mangoldt a -> a+runMangoldt m = case m of+  MangoldtZero  -> 1+  MangoldtOne a -> a+  MangoldtMany  -> 1++instance Semigroup (Mangoldt a) where+  MangoldtZero <> a = a+  a <> MangoldtZero = a+  _ <> _ = MangoldtMany++instance Monoid (Mangoldt a) where+  mempty  = MangoldtZero+  mappend = (<>)++newtype LCM a = LCM { getLCM :: a }++instance Integral a => Semigroup (LCM a) where+  (<>) = coerce (lcm :: a -> a -> a)++instance Integral a => Monoid (LCM a) where+  mempty  = LCM 1+  mappend = (<>)++newtype Xor = Xor { _getXor :: Bool }++runXor :: Num a => Xor -> a+runXor m = case m of+  Xor False ->  1+  Xor True  -> -1++instance Semigroup Xor where+   (<>) = coerce ((/=) :: Bool -> Bool -> Bool)++instance Monoid Xor where+  mempty  = Xor False+  mappend = (<>)++newtype SetProduct a = SetProduct { getSetProduct :: Set a }++instance (Num a, Ord a) => Semigroup (SetProduct a) where+  SetProduct s1 <> SetProduct s2 = SetProduct $ s1 <> s2 <> foldMap (\n -> S.mapMonotonic (* n) s2) s1++instance (Num a, Ord a) => Monoid (SetProduct a) where+  mempty  = SetProduct mempty+  mappend = (<>)++newtype IntSetProduct = IntSetProduct { getIntSetProduct :: IntSet }++instance Semigroup IntSetProduct where+  IntSetProduct s1 <> IntSetProduct s2 = IntSetProduct $ IS.unions $ s1 : s2 : map (\n -> IS.map (* n) s2) (IS.toAscList s1)++instance Monoid IntSetProduct where+  mempty  = IntSetProduct mempty+  mappend = (<>)

@@ -10,7 +10,8 @@ -- Normally, accessing these via the higher level interface of "Math.NumberTheory.GCD" -- should be sufficient. ---{-# LANGUAGE CPP, MagicHash, UnboxedTuples #-}+{-# LANGUAGE MagicHash     #-}+{-# LANGUAGE UnboxedTuples #-} module Math.NumberTheory.GCD.LowLevel   ( -- * Specialised GCDs     gcdInt@@ -27,10 +28,6 @@   ) where  import GHC.Base--#if __GLASGOW_HASKELL__ < 705-import GHC.Word (Word(..))      -- Moved to GHC.Types-#endif  import Math.NumberTheory.Utils 

@@ -11,6 +11,7 @@ --  {-# LANGUAGE BangPatterns #-}+ module Math.NumberTheory.GaussianIntegers (     GaussianInteger((:+)),     ι,@@ -187,7 +188,7 @@     where     s = a .^ div e 2 --- |Compute the prime factorization of a Gaussian integer. This is unique up to units (+/- 1, +/- i).+-- |Compute the prime factorisation of a Gaussian integer. This is unique up to units (+/- 1, +/- i). factorise :: GaussianInteger -> [(GaussianInteger, Int)] factorise g     | g == 0    = error "0 has no prime factorisation"

@@ -30,14 +30,11 @@  import GHC.Base -#if __GLASGOW_HASKELL__ < 705-import GHC.Word (Word(..))      -- Moved to GHC.Types-#endif- import Data.Bits import Data.Array.Unboxed -import Math.NumberTheory.Logarithms.Internal+import GHC.Integer.Logarithms+ import Math.NumberTheory.Powers.Integer import Math.NumberTheory.Unsafe #if __GLASGOW_HASKELL__ < 707

@@ -1,155 +0,0 @@--- |--- Module:      Math.NumberTheory.Logarithms.Internal--- Copyright:   (c) 2011 Daniel Fischer--- Licence:     MIT--- Maintainer:  Daniel Fischer <daniel.is.fischer@googlemail.com>--- Stability:   Provisional--- Portability: Non-portable (GHC extensions)------ Low level stuff for integer logarithms.-{-# LANGUAGE CPP, MagicHash, UnboxedTuples #-}-{-# OPTIONS_HADDOCK hide #-}-module Math.NumberTheory.Logarithms.Internal-    ( -- * Functions-      integerLog2#-    , wordLog2#-    ) where--#if __GLASGOW_HASKELL__ >= 702---- Stuff is already there-import GHC.Integer.Logarithms--#else---- We have to define it here-#include "MachDeps.h"--import GHC.Base-import GHC.Integer.GMP.Internals--#if (WORD_SIZE_IN_BITS != 32) && (WORD_SIZE_IN_BITS != 64)-#error Only word sizes 32 and 64 are supported.-#endif---#if WORD_SIZE_IN_BITS == 32--#define WSHIFT 5-#define MMASK 31--#else--#define WSHIFT 6-#define MMASK 63--#endif--{---Reference implementation only, the algorithm in M.NT.Logarithms is better.---- | Calculate the integer logarithm for an arbitrary base.---   The base must be greater than 1, the second argument, the number---   whose logarithm is sought; should be positive, otherwise the---   result is meaningless.------ > base ^ integerLogBase# base m <= m < base ^ (integerLogBase# base m + 1)------ for @base > 1@ and @m > 0@.-integerLogBase# :: Integer -> Integer -> Int#-integerLogBase# b m = case step b of-                        (# _, e #) -> e-  where-    step pw =-      if m `ltInteger` pw-        then (# m, 0# #)-        else case step (pw `timesInteger` pw) of-               (# q, e #) ->-                 if q `ltInteger` pw-                   then (# q, 2# *# e #)-                   else (# q `quotInteger` pw, 2# *# e +# 1# #)--}---- | Calculate the integer base 2 logarithm of an 'Integer'.---   The calculation is much more efficient than for the general case.------   The argument must be strictly positive, that condition is /not/ checked.-integerLog2# :: Integer -> Int#-integerLog2# (S# i) = wordLog2# (int2Word# i)-integerLog2# (J# s ba) = check (s -# 1#)-  where-    check i = case indexWordArray# ba i of-                0## -> check (i -# 1#)-                w   -> wordLog2# w +# (uncheckedIShiftL# i WSHIFT#)---- | This function calculates the integer base 2 logarithm of a 'Word#'.---   @'wordLog2#' 0## = -1#@.-{-# INLINE wordLog2# #-}-wordLog2# :: Word# -> Int#-wordLog2# w =-  case leadingZeros of-   BA lz ->-    let zeros u = indexInt8Array# lz (word2Int# u) in-#if WORD_SIZE_IN_BITS == 64-    case uncheckedShiftRL# w 56# of-     a ->-      if a `neWord#` 0##-       then 64# -# zeros a-       else-        case uncheckedShiftRL# w 48# of-         b ->-          if b `neWord#` 0##-           then 56# -# zeros b-           else-            case uncheckedShiftRL# w 40# of-             c ->-              if c `neWord#` 0##-               then 48# -# zeros c-               else-                case uncheckedShiftRL# w 32# of-                 d ->-                  if d `neWord#` 0##-                   then 40# -# zeros d-                   else-#endif-                    case uncheckedShiftRL# w 24# of-                     e ->-                      if e `neWord#` 0##-                       then 32# -# zeros e-                       else-                        case uncheckedShiftRL# w 16# of-                         f ->-                          if f `neWord#` 0##-                           then 24# -# zeros f-                           else-                            case uncheckedShiftRL# w 8# of-                             g ->-                              if g `neWord#` 0##-                               then 16# -# zeros g-                               else 8# -# zeros w---- Lookup table-data BA = BA ByteArray#--leadingZeros :: BA-leadingZeros =-    let mkArr s =-          case newByteArray# 256# s of-            (# s1, mba #) ->-              case writeInt8Array# mba 0# 9# s1 of-                s2 ->-                  let fillA lim val idx st =-                        if idx ==# 256#-                          then st-                          else if idx <# lim-                                then case writeInt8Array# mba idx val st of-                                        nx -> fillA lim val (idx +# 1#) nx-                                else fillA (2# *# lim) (val -# 1#) idx st-                  in case fillA 2# 8# 1# s2 of-                      s3 -> case unsafeFreezeByteArray# mba s3 of-                              (# _, ba #) -> ba-    in case mkArr realWorld# of-        b -> BA b--#endif

@@ -24,15 +24,12 @@ import Data.Array.ST  import Data.Bits-#if __GLASGOW_HASKELL__ < 705-import Data.Word-#endif  import GHC.Base import GHC.Integer import GHC.Integer.GMP.Internals+import GHC.Integer.Logarithms (integerLog2#) -import Math.NumberTheory.Logarithms.Internal (integerLog2#) import Math.NumberTheory.Unsafe #if __GLASGOW_HASKELL__ < 707 import Math.NumberTheory.Utils (isTrue#)

@@ -23,16 +23,13 @@ import GHC.Base import GHC.Integer import GHC.Integer.GMP.Internals+import GHC.Integer.Logarithms (integerLog2#)  import Data.Array.Unboxed import Data.Array.ST  import Data.Bits-#if __GLASGOW_HASKELL__ < 705-import Data.Word-#endif -import Math.NumberTheory.Logarithms.Internal (integerLog2#) import Math.NumberTheory.Unsafe #if __GLASGOW_HASKELL__ < 707 import Math.NumberTheory.Utils (isTrue#)

@@ -25,16 +25,13 @@ import GHC.Base import GHC.Integer import GHC.Integer.GMP.Internals+import GHC.Integer.Logarithms (integerLog2#)  import Data.Bits-#if __GLASGOW_HASKELL__ < 705-import Data.Word-#endif import Data.List (foldl') import qualified Data.Set as Set  import Math.NumberTheory.Logarithms (integerLogBase')-import Math.NumberTheory.Logarithms.Internal (integerLog2#) import Math.NumberTheory.Utils  (shiftToOddCount                                 , splitOff #if __GLASGOW_HASKELL__ < 707

@@ -16,11 +16,8 @@     ) where  import GHC.Base-#if __GLASGOW_HASKELL__ < 705-import GHC.Word-#endif+import GHC.Integer.Logarithms (wordLog2#) -import Math.NumberTheory.Logarithms.Internal ( wordLog2# ) #if __GLASGOW_HASKELL__ < 707 import Math.NumberTheory.Utils (isTrue#) #endif

@@ -8,11 +8,13 @@ -- -- Functions dealing with squares. Efficient calculation of integer square roots -- and efficient testing for squareness.-{-# LANGUAGE MagicHash, BangPatterns, CPP, FlexibleContexts #-}+{-# LANGUAGE MagicHash, BangPatterns, PatternGuards, CPP, FlexibleContexts #-} module Math.NumberTheory.Powers.Squares     ( -- * Square root calculation       integerSquareRoot     , integerSquareRoot'+    , integerSquareRootRem+    , integerSquareRootRem'     , exactSquareRoot       -- * Tests for squares     , isSquare@@ -23,24 +25,17 @@  #include "MachDeps.h" -import GHC.Base-import GHC.Integer-import GHC.Integer.GMP.Internals- import Data.Array.Unboxed import Data.Array.ST  import Data.Bits-#if __GLASGOW_HASKELL__ < 705+#if __GLASGOW_HASKELL__ < 709 import Data.Word        -- Moved to GHC.Types #endif -import Math.NumberTheory.Logarithms.Internal (integerLog2#) import Math.NumberTheory.Unsafe-#if __GLASGOW_HASKELL__ < 707-import Math.NumberTheory.Utils (isTrue#)-#endif +import Math.NumberTheory.Powers.Squares.Internal  -- | Calculate the integer square root of a nonnegative number @n@, --   that is, the largest integer @r@ with @r*r <= n@.@@ -58,13 +53,40 @@ --   that is, the largest integer @r@ with @r*r <= n@. --   The precondition @n >= 0@ is not checked. {-# RULES-"integerSquareRoot'/Int"  integerSquareRoot' = isqrtInt'-"integerSquareRoot'/Word" integerSquareRoot' = isqrtWord+"integerSquareRoot'/Int"     integerSquareRoot' = isqrtInt'+"integerSquareRoot'/Word"    integerSquareRoot' = isqrtWord+"integerSquareRoot'/Integer" integerSquareRoot' = isqrtInteger   #-} {-# INLINE [1] integerSquareRoot' #-} integerSquareRoot' :: Integral a => a -> a integerSquareRoot' = isqrtA +-- | Calculate the integer square root of a nonnegative number as well as+--   the difference of that number with the square of that root, that is if+--   @(s,r) = integerSquareRootRem n@ then @s^2 <= n == s^2+r < (s+1)^2@.+{-# SPECIALISE integerSquareRootRem ::+        Int -> (Int, Int),+        Word -> (Word, Word),+        Integer -> (Integer, Integer)+  #-}+integerSquareRootRem :: Integral a => a -> (a, a)+integerSquareRootRem n+  | n < 0       = error "integerSquareRootRem: negative argument"+  | otherwise   = integerSquareRootRem' n++-- | Calculate the integer square root of a nonnegative number as well as+--   the difference of that number with the square of that root, that is if+--   @(s,r) = integerSquareRootRem' n@ then @s^2 <= n == s^2+r < (s+1)^2@.+--   The precondition @n >= 0@ is not checked.+{-# RULES+"integerSquareRootRem'/Integer" integerSquareRootRem' = karatsubaSqrt+  #-}+{-# INLINE [1] integerSquareRootRem' #-}+integerSquareRootRem' :: Integral a => a -> (a, a)+integerSquareRootRem' n = (s, n - s * s)+  where+    s = integerSquareRoot' n+ -- | Returns 'Nothing' if the argument is not a square, --   @'Just' r@ if @r*r == n@ and @r >= 0@. Avoids the expensive calculation --   of the square root if @n@ is recognized as a non-square@@ -77,11 +99,10 @@   #-} exactSquareRoot :: Integral a => a -> Maybe a exactSquareRoot n-  | n < 0                           = Nothing-  | isPossibleSquare n && r*r == n  = Just r-  | otherwise                       = Nothing-    where-      r = integerSquareRoot' n+  | n >= 0+  , isPossibleSquare n+  , (r, 0) <- integerSquareRootRem' n = Just r+  | otherwise                         = Nothing  -- | Test whether the argument is a square. --   After a number is found to be positive, first 'isPossibleSquare'@@ -104,7 +125,10 @@                             Integer -> Bool   #-} isSquare' :: Integral a => a -> Bool-isSquare' n = isPossibleSquare n && let r = integerSquareRoot' n in r*r == n+isSquare' n+    | isPossibleSquare n+    , (_, 0) <- integerSquareRootRem' n = True+    | otherwise                         = False  -- | Test whether a non-negative number may be a square. --   Non-negativity is not checked, passing negative arguments may@@ -157,63 +181,6 @@ ----------------------------------------------------------------------------- --  Auxiliary Stuff --- Find approximation to square root in 'Integer', then--- find the integer square root by the integer variant--- of Heron's method. Takes only a handful of steps--- unless the input is really large.-{-# SPECIALISE isqrtA :: Integer -> Integer #-}-isqrtA :: Integral a => a -> a-isqrtA 0 = 0-isqrtA n = heron n (fromInteger . appSqrt . fromIntegral $ n)---- Heron's method for integers. First make one step to ensure--- the value we're working on is @>= r@, then we have--- @k == r@ iff @k <= step k@.-{-# SPECIALISE heron :: Integer -> Integer -> Integer #-}-heron :: Integral a => a -> a -> a-heron n a = go (step a)-      where-        step k = (k + n `quot` k) `quot` 2-        go k-            | m < k     = go m-            | otherwise = k-              where-                m = step k---- threshold for shifting vs. direct fromInteger--- we shift when we expect more than 256 bits-#if WORD_SIZE_IN_BITS == 64-#define THRESH 5-#else-#define THRESH 9-#endif---- Find a fairly good approximation to the square root.--- At most one off for small Integers, about 48 bits should be correct--- for large Integers.-appSqrt :: Integer -> Integer-appSqrt (S# i#) = S# (double2Int# (sqrtDouble# (int2Double# i#)))-#if __GLASGOW_HASKELL__ < 709-appSqrt n@(J# s# _)-    | isTrue# (s# <# THRESH#) = floor (sqrt $ fromInteger n :: Double)-#else-appSqrt n@(Jp# bn#)-    | isTrue# ((sizeofBigNat# bn#) <# THRESH#) =-          floor (sqrt $ fromInteger n :: Double)-#endif-    | otherwise = case integerLog2# n of-                    l# -> case uncheckedIShiftRA# l# 1# -# 47# of-                            h# -> case shiftRInteger n (2# *# h#) of-                                    m -> case floor (sqrt $ fromInteger m :: Double) of-                                            r -> shiftLInteger r h#-#if __GLASGOW_HASKELL__ >= 709--- There's already a check for negative in integerSquareRoot,--- but integerSquareRoot' is exported directly too.-appSqrt _ = error "integerSquareRoot': negative argument"-#endif---- Auxiliaries- -- Make an array indicating whether a remainder is a square remainder. sqRemArray :: Int -> UArray Int Bool sqRemArray md = runSTUArray $ do@@ -250,7 +217,7 @@ sr325 :: UArray Int Bool sr325 = sqRemArray 325 --- Specialisations for Int and Word+-- Specialisations for Int, Word, and Integer  -- For @n <= 2^64@, the result of --@@ -286,3 +253,6 @@       where         !r = (fromIntegral :: Int -> Word) . (truncate :: Double -> Int) . sqrt $ fromIntegral n +{-# INLINE isqrtInteger #-}+isqrtInteger :: Integer -> Integer+isqrtInteger = fst . karatsubaSqrt

@@ -0,0 +1,144 @@+-- |+-- Module:      Math.NumberTheory.Powers.Squares.Internal+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Internal functions dealing with square roots. End-users should not import this module.++{-# LANGUAGE MagicHash        #-}+{-# LANGUAGE BangPatterns     #-}+{-# LANGUAGE PatternGuards    #-}+{-# LANGUAGE CPP              #-}+{-# LANGUAGE FlexibleContexts #-}++module Math.NumberTheory.Powers.Squares.Internal+  ( karatsubaSqrt+  , isqrtA+  ) where++#include "MachDeps.h"++import Data.Bits++import GHC.Base+import GHC.Integer+import GHC.Integer.GMP.Internals+import GHC.Integer.Logarithms (integerLog2#)++import Math.NumberTheory.Logarithms (integerLog2)+#if __GLASGOW_HASKELL__ < 707+import Math.NumberTheory.Utils (isTrue#)+#endif++-- Find approximation to square root in 'Integer', then+-- find the integer square root by the integer variant+-- of Heron's method. Takes only a handful of steps+-- unless the input is really large.+{-# SPECIALISE isqrtA :: Integer -> Integer #-}+isqrtA :: Integral a => a -> a+isqrtA 0 = 0+isqrtA n = heron n (fromInteger . appSqrt . fromIntegral $ n)++-- Heron's method for integers. First make one step to ensure+-- the value we're working on is @>= r@, then we have+-- @k == r@ iff @k <= step k@.+{-# SPECIALISE heron :: Integer -> Integer -> Integer #-}+heron :: Integral a => a -> a -> a+heron n a = go (step a)+      where+        step k = (k + n `quot` k) `quot` 2+        go k+            | m < k     = go m+            | otherwise = k+              where+                m = step k++-- threshold for shifting vs. direct fromInteger+-- we shift when we expect more than 256 bits+#if WORD_SIZE_IN_BITS == 64+#define THRESH 5+#else+#define THRESH 9+#endif++-- Find a fairly good approximation to the square root.+-- At most one off for small Integers, about 48 bits should be correct+-- for large Integers.+appSqrt :: Integer -> Integer+appSqrt (S# i#) = S# (double2Int# (sqrtDouble# (int2Double# i#)))+#if __GLASGOW_HASKELL__ < 709+appSqrt n@(J# s# _)+    | isTrue# (s# <# THRESH#) = floor (sqrt $ fromInteger n :: Double)+#else+appSqrt n@(Jp# bn#)+    | isTrue# ((sizeofBigNat# bn#) <# THRESH#) =+          floor (sqrt $ fromInteger n :: Double)+#endif+    | otherwise = case integerLog2# n of+                    l# -> case uncheckedIShiftRA# l# 1# -# 47# of+                            h# -> case shiftRInteger n (2# *# h#) of+                                    m -> case floor (sqrt $ fromInteger m :: Double) of+                                            r -> shiftLInteger r h#+#if __GLASGOW_HASKELL__ >= 709+-- There's already a check for negative in integerSquareRoot,+-- but integerSquareRoot' is exported directly too.+appSqrt _ = error "integerSquareRoot': negative argument"+#endif+++-- Integer square root with remainder, using the Karatsuba Square Root+-- algorithm from+-- Paul Zimmermann. Karatsuba Square Root. [Research Report] RR-3805, 1999,+-- pp.8. <inria-00072854>++karatsubaSqrt :: Integer -> (Integer, Integer)+karatsubaSqrt 0 = (0, 0)+karatsubaSqrt n+    | lgN < 2300 =+        let s = isqrtA n in (s, n - s * s)+    | otherwise =+        if lgN .&. 2 /= 0 then+            karatsubaStep k (karatsubaSplit k n)+        else+            -- before we split n into 4 part we must ensure that the first part+            -- is at least 2^k/4, since this doesn't happen here we scale n by+            -- multiplying it by 4+            let n' = n `unsafeShiftL` 2+                (s, r) = karatsubaStep k (karatsubaSplit k n')+                r' | s .&. 1 == 0 = r+                   | otherwise = r + double s - 1+            in  (s `unsafeShiftR` 1, r' `unsafeShiftR` 2)+  where+    k = lgN `unsafeShiftR` 2 + 1+    lgN = integerLog2 n++karatsubaStep :: Int -> (Integer, Integer, Integer, Integer) -> (Integer, Integer)+karatsubaStep k (a3, a2, a1, a0)+    | r >= 0 = (s, r)+    | otherwise = (s - 1, r + double s - 1)+  where+    r = cat u a0 - q * q+    s = s' `unsafeShiftL` k + q+    (q, u) = cat r' a1 `quotRem` double s'+    (s', r') = karatsubaSqrt (cat a3 a2)+    cat x y = x `unsafeShiftL` k .|. y+    {-# INLINE cat #-}++karatsubaSplit :: Int -> Integer -> (Integer, Integer, Integer, Integer)+karatsubaSplit k n0 = (a3, a2, a1, a0)+  where+    a3 = n3+    n3 = n2 `unsafeShiftR` k+    a2 = n2 .&. m+    n2 = n1 `unsafeShiftR` k+    a1 = n1 .&. m+    n1 = n0 `unsafeShiftR` k+    a0 = n0 .&. m+    m = 1 `unsafeShiftL` k - 1++double :: Bits a => a -> a+double x = x `unsafeShiftL` 1+{-# INLINE double #-}

@@ -9,9 +9,7 @@ -- Number of primes not exceeding @n@, @π(n)@, and @n@-th prime. -- {-# LANGUAGE CPP, BangPatterns, FlexibleContexts #-}-#if __GLASGOW_HASKELL__ >= 700 {-# OPTIONS_GHC -fspec-constr-count=24 #-}-#endif {-# OPTIONS_HADDOCK hide #-} module Math.NumberTheory.Primes.Counting.Impl     ( primeCount@@ -31,9 +29,6 @@ import Math.NumberTheory.Unsafe  import Data.Array.ST-#if !MIN_VERSION_array(0,5,0)-    hiding (unsafeThaw)-#endif import Control.Monad.ST import Data.Bits import Data.Int

@@ -13,6 +13,9 @@ -- on the canonical factorisation, these require that the number be positive -- and in the case of the Carmichael function that the list of prime factors -- with their multiplicities is ascending.++{-# OPTIONS_GHC -fno-warn-warnings-deprecations #-}+ module Math.NumberTheory.Primes.Factorisation     ( -- * Factorisation functions       -- $algorithm@@ -73,6 +76,8 @@ import Math.NumberTheory.Primes.Factorisation.Montgomery import Math.NumberTheory.Primes.Factorisation.TrialDivision import Math.NumberTheory.Primes.Sieve.Misc++{-# DEPRECATED totient, φ, carmichael, λ, moebius, μ, divisors, tau, τ, divisorCount, divisorSum, sigma, σ, divisorPowerSum "Use 'Math.NumberTheory.ArithmeticFunctions'" #-}  -- $algorithm --

@@ -43,9 +43,6 @@ #include "MachDeps.h"  import GHC.Base-#if __GLASGOW_HASKELL__ < 705-import GHC.Word     -- Moved to GHC.Types-#endif  import System.Random import Control.Monad.State.Strict@@ -55,8 +52,9 @@ import Data.Bits import Data.Maybe +import GHC.Integer.Logarithms+ import Math.NumberTheory.Logarithms-import Math.NumberTheory.Logarithms.Internal import Math.NumberTheory.Powers.General     (highestPower, largePFPower) import Math.NumberTheory.Powers.Squares     (integerSquareRoot') import Math.NumberTheory.Primes.Sieve.Eratosthenes

@@ -27,6 +27,8 @@  import Math.NumberTheory.Powers.Integer +{-# DEPRECATED totientFromCanonical, carmichaelFromCanonical, moebiusFromCanonical, divisorsFromCanonical, tauFromCanonical, divisorSumFromCanonical, sigmaFromCanonical "Use 'Math.NumberTheory.ArithmeticFunctions'" #-}+ -- | Totient of a prime power. ppTotient :: (Integer,Int) -> Integer ppTotient (p,1) = p-1

@@ -16,9 +16,7 @@ -- This module is mainly intended for comparison and verification. {-# LANGUAGE BangPatterns, CPP, MonoLocalBinds #-} {-# OPTIONS_GHC -funbox-strict-fields #-}-#if __GLASGOW_HASKELL__ >= 700 {-# OPTIONS_GHC -fno-float-in -fno-spec-constr -fno-full-laziness #-}-#endif module Math.NumberTheory.Primes.Heap (primes, sieveFrom) where  import Data.Array.Unboxed@@ -58,10 +56,6 @@ {-# SPECIALISE push :: Word -> Word -> Int -> Hipp Word -> Hipp Word #-} {-# SPECIALISE push :: Integer -> Integer -> Int -> Hipp Integer -> Hipp Integer #-} push :: Integral a => a -> a -> Int -> Hipp a -> Hipp a--- GHC 7 does not like the old code, so it gets a new implementation.--- That is faster than what it does with the old code, but still slower--- than what GHC 6 did with it. :(-#if __GLASGOW_HASKELL__ >= 700 push !c !p !w = go   where     less = (< c)@@ -69,21 +63,12 @@         | less hc   = H hc hp hw (go r) l         | otherwise = H c p w (push hc hp hw r) l     go _ = H c p w E E-#else-push c p w (H hc hp hw l r)-    | c < hc    = H c p w (push hc hp hw r) l-    | otherwise = H hc hp hw (push c p w r) l-push c p w E = H c p w E E-#endif  -- bubble down increased top to regain heap invariant {-# SPECIALISE bubble :: Hipp Int -> Hipp Int #-} {-# SPECIALISE bubble :: Hipp Word -> Hipp Word #-} {-# SPECIALISE bubble :: Hipp Integer -> Hipp Integer #-} bubble :: Integral a => Hipp a -> Hipp a--- Again, GHC 6 fared better, so new code for GHC 7, still--- not as good as GHC 6 was.-#if __GLASGOW_HASKELL__ >= 700 bubble h@(H c p w l r) =     case r of         E -> case l of@@ -100,19 +85,9 @@                 | rc < c -> H rc rp rw l (mkHipp c p w rl rr)                 | otherwise -> h             _ -> error "Heap invariant violated, left smaller than right!"-#else-bubble h@(H c p w l@(H lc lp lw ll lr) r@(H rc rp rw rl rr))-    | c <= lc && c <= rc = h-    | lc < rc   = H lc lp lw (bubble (H c p w ll lr)) r-    | otherwise = H rc rp rw l (bubble (H c p w rl rr))-bubble h@(H c p w (H lc lp lw _ _) _)-    | c <= lc   = h-    | otherwise = H lc lp lw (H c p w E E) E-#endif bubble h = h --- join two heaps and composite-data, GHC 7 doesn't do well on the old bubble.-#if __GLASGOW_HASKELL__ >= 700+-- join two heaps and composite-data {-# SPECIALISE     mkHipp :: Int -> Int -> Int -> Hipp Int -> Hipp Int -> Hipp Int,               Integer -> Integer -> Int -> Hipp Integer -> Hipp Integer -> Hipp Integer,@@ -138,8 +113,6 @@                 | less rc -> H rc rp rw l (go rl rr)                 | otherwise -> H c p w l r             _ -> error "Heap invariant violated, left smaller than right!"--- {-# INLINE mkHipp #-}-#endif  -- increase the top of the heap and re-heap {-# SPECIALISE inc :: Hipp Int -> Hipp Int #-}

@@ -9,9 +9,7 @@ -- Sieve -- {-# LANGUAGE CPP, BangPatterns, FlexibleContexts #-}-#if __GLASGOW_HASKELL__ >= 700 {-# OPTIONS_GHC -fspec-constr-count=8 #-}-#endif {-# OPTIONS_HADDOCK hide #-} module Math.NumberTheory.Primes.Sieve.Eratosthenes     ( primes@@ -36,9 +34,6 @@  import Control.Monad.ST import Data.Array.ST-#if !MIN_VERSION_array(0,5,0)-                     hiding (unsafeFreeze, unsafeThaw, castSTUArray)-#endif import Control.Monad (when) import Data.Bits #if __GLASGOW_HASKELL__ < 709 || WORD_SIZE_IN_BITS == 32

@@ -6,10 +6,12 @@ -- Stability:   Provisional -- Portability: Non-portable (GHC extensions) ---{-# LANGUAGE CPP, BangPatterns, ScopedTypeVariables, MonoLocalBinds, FlexibleContexts #-}-#if __GLASGOW_HASKELL__ >= 700+{-# LANGUAGE BangPatterns        #-}+{-# LANGUAGE FlexibleContexts    #-}+{-# LANGUAGE MonoLocalBinds      #-}+{-# LANGUAGE ScopedTypeVariables #-} {-# OPTIONS_GHC -fspec-constr-count=8 #-}-#endif+{-# OPTIONS_GHC -fno-warn-warnings-deprecations #-} {-# OPTIONS_HADDOCK hide #-} module Math.NumberTheory.Primes.Sieve.Misc     ( -- * Types@@ -46,6 +48,27 @@ import Math.NumberTheory.Unsafe import Math.NumberTheory.Utils +{-+IMPORTANT NOTICE: Not all sieves use the same layout!++FactorSieve:++   To remain as efficient as possible, FactorSieve omits only even numbers.+   To relate an odd number x to its index i:++       i = (x `div` 2) - 1+       x = i * 2 + 3++TotientSieve, CarmichaelSieve:++   These sieves use a (2,3,5) wheel optimization, sacrificing performance to save+   more memory. The only indices stored are those coprime to 2, 3, and 5.+   To relate such an integer x to its index i:++       i = toIdx x+       x = toPrim i+-}+ -- | A compact store of smallest prime factors. data FactorSieve = FS {-# UNPACK #-} !Word {-# UNPACK #-} !(UArray Int Word16) @@ -148,8 +171,8 @@                                                j | j <= bound -> intLoop (fromIntegral (j `shiftR` 1) - 1)                                                  | otherwise -> tdLoop j (integerSquareRoot' j) (ix+1)           where-            p = toPrim ix-            pix = unsafeAt sve $ fromIntegral p+            p = fromIntegral $ 2 * ix + 3+            pix = unsafeAt sve ix     curve n = stdGenFactorisation (Just (bound*(bound+2))) (mkStdGen $ fromIntegral n `xor` 0xdecaf00d) Nothing n  -- | @'totientSieve' n@ creates a store of the totients of the numbers not exceeding @n@.@@ -277,6 +300,9 @@     curve tt n = tt `lcm` carmichaelFromCanonical (stdGenFactorisation (Just (bound*(bound+2))) (mkStdGen $ fromIntegral n `xor` 0xdecaf00d) Nothing n)  +-- NOTE: This is a legacy implementation of FactorSieve which uses the+--       same (2,3,5) wheel optimization as the other sieves.+--       It is still used to generate the other sieves. spfSieve :: Word -> ST s (STUArray s Int Word) spfSieve bound = do   let (octs,lidx) = idxPr bound

@@ -28,9 +28,6 @@  import GHC.Base -#if __GLASGOW_HASKELL__ < 705-import GHC.Word     -- Moved to GHC.Types-#endif import GHC.Integer.GMP.Internals  -- | @'isPrime' n@ tests whether @n@ is a prime (negative or positive).

@@ -0,0 +1,131 @@+-- |+-- Module:      Math.NumberTheory.UniqueFactorisation+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- An abstract type class for unique factorisation domains.+--++{-# LANGUAGE CPP          #-}+{-# LANGUAGE TypeFamilies #-}++module Math.NumberTheory.UniqueFactorisation+  ( Prime+  , UniqueFactorisation(..)+  ) where++import Control.Arrow++#if MIN_VERSION_base(4,8,0)+#else+import Data.Word+#endif++import Math.NumberTheory.Primes.Factorisation as F (factorise')+import Math.NumberTheory.GaussianIntegers as G++import Numeric.Natural++#if MIN_VERSION_base(4,7,0)+import Data.Coerce+#else+import Unsafe.Coerce++coerce :: a -> b+coerce = unsafeCoerce+#endif++newtype SmallPrime = SmallPrime { _unSmallPrime :: Word }+  deriving (Eq, Ord, Show)++newtype BigPrime = BigPrime { _unBigPrime :: Natural }+  deriving (Eq, Ord, Show)++-- | Type of primes of a given unique factorisation domain.+-- When the domain has exactly one unit, @Prime t = t@,+-- but when units are multiple more restricted types+-- (or at least newtypes) should be specified.+--+-- @abs (unPrime n) == unPrime n@ must hold for all @n@ of type @Prime t@+type family Prime (f :: *) :: *++type instance Prime Int     = SmallPrime+type instance Prime Word    = SmallPrime+type instance Prime Integer = BigPrime+type instance Prime Natural = BigPrime++type instance Prime G.GaussianInteger = GaussianPrime++-- | The following invariant must hold for @n /= 0@:+--+-- > abs n == abs (product (map (\(p, k) -> unPrime p ^ k) (factorise n)))+--+-- The result of 'factorise' should not contain zero powers and should not change after multiplication of the argument by domain's unit.+class UniqueFactorisation a where+  unPrime   :: Prime a -> a+  factorise :: a -> [(Prime a, Word)]++instance UniqueFactorisation Int where+  unPrime   = coerce wordToInt+  factorise m' = if m <= 1+                then []+                else map (coerce integerToWord *** intToWord) . F.factorise' . intToInteger $ m+                  where+                    m = abs m'++instance UniqueFactorisation Word where+  unPrime     = coerce+  factorise m = if m <= 1+                  then []+                  else map (coerce integerToWord *** intToWord) . F.factorise' . wordToInteger $ m+++instance UniqueFactorisation Integer where+  unPrime      = coerce naturalToInteger+  factorise m' = if m <= 1+                then []+                else map (coerce integerToNatural *** intToWord) . F.factorise' $ m+                  where+                    m = abs m'++instance UniqueFactorisation Natural where+  unPrime   = coerce+  factorise m = if m <= 1+                  then []+                  else map (coerce integerToNatural *** intToWord) . F.factorise' . naturalToInteger $ m++newtype GaussianPrime = GaussianPrime { _unGaussianPrime :: G.GaussianInteger }+  deriving (Eq, Show)++instance UniqueFactorisation G.GaussianInteger where+  unPrime = coerce++  factorise 0 = []+  factorise g = map (coerce *** intToWord) $ filter (\(h, _) -> abs h /= 1) $ G.factorise g++-----------+-- Utils++wordToInt :: Word -> Int+wordToInt = fromIntegral++wordToInteger :: Word -> Integer+wordToInteger = fromIntegral++intToWord :: Int -> Word+intToWord = fromIntegral++intToInteger :: Int -> Integer+intToInteger = fromIntegral++naturalToInteger :: Natural -> Integer+naturalToInteger = fromIntegral++integerToNatural :: Integer -> Natural+integerToNatural = fromIntegral++integerToWord :: Integer -> Word+integerToWord = fromIntegral

@@ -34,9 +34,6 @@   , (!)   ) import Data.Array.MArray-#if !MIN_VERSION_array(0,5,0)-  hiding (unsafeFreeze, unsafeThaw)-#endif  unsafeAt :: (IArray a e, Ix i) => a i e -> i -> e unsafeAt = (!)

@@ -28,9 +28,6 @@ #include "MachDeps.h"  import GHC.Base-#if __GLASGOW_HASKELL__ < 705-import GHC.Word     -- Word and its constructor moved to GHC.Types-#endif  import GHC.Integer import GHC.Integer.GMP.Internals@@ -178,24 +175,11 @@  -- | Number of 1-bits in a @'Word'@. bitCountWord :: Word -> Int-#if __GLASGOW_HASKELL__ >= 703 bitCountWord = popCount--- should yield a machine instruction-#else-bitCountWord w = case w - (shiftR w 1 .&. m5) of-                   !w1 -> case (w1 .&. m3) + (shiftR w1 2 .&. m3) of-                            !w2 -> case (w2 + shiftR w2 4) .&. mf of-                                     !w3 -> fromIntegral (shiftR (w3 * m1) sd)-#endif  -- | Number of 1-bits in an @'Int'@. bitCountInt :: Int -> Int-#if __GLASGOW_HASKELL__ >= 703 bitCountInt = popCount--- should yield a machine instruction-#else-bitCountInt (I# i#) = bitCountWord (W# (int2Word# i#))-#endif  -- | Number of trailing zeros in a @'Word#'@, wrong for @0@. {-# INLINE trailZeros# #-}@@ -214,11 +198,7 @@ --                            Int -> Int -> (Int, Int), --                            Word -> Word -> (Int, Word) --   #-}-#if __GLASGOW_HASKELL__ >= 700 {-# INLINABLE splitOff #-}-#else-{-# INLINE splitOff #-}-#endif splitOff :: Integral a => a -> a -> (Int, a) splitOff p n = go 0 n   where

@@ -1,5 +1,5 @@ name                : arithmoi-version             : 0.4.2.0+version             : 0.4.3.0 cabal-version       : >= 1.10 author              : Daniel Fischer copyright           : (c) 2011 Daniel Fischer@@ -21,10 +21,6 @@                       powers (integer roots and tests, modular exponentiation),                       integer logarithms. -                      Note: Requires GHC >= 6.12 with the integer-gmp package-                      for efficiency. Portability is on the to-do list (with-                      low priority, however).- category            : Math, Algorithms, Number Theory  tested-with         : GHC==7.6.3, GHC==7.8.4, GHC==7.10.3, GHC==8.0.1@@ -38,15 +34,22 @@  library     default-language: Haskell2010-    build-depends       : base >= 4 && < 5-                          , array >= 0.3 && < 0.6-                          , ghc-prim < 0.6-                          , integer-gmp < 1.1-                          , containers >= 0.3 && < 0.6-                          , random >= 1.0 && < 1.2-                          , mtl >= 2.0 && < 2.3+    build-depends       : base >= 4.6 && < 5+                        , array >= 0.5 && < 0.6+                        , ghc-prim < 0.6+                        , integer-gmp < 1.1+                        , containers >= 0.5 && < 0.6+                        , random >= 1.0 && < 1.2+                        , mtl >= 2.0 && < 2.3+    if impl(ghc < 7.10)+      build-depends     : nats >= 1 && <1.2+    if impl(ghc < 8.0)+      build-depends     : semigroups >= 0.8 -    exposed-modules     : Math.NumberTheory.Logarithms+    exposed-modules     : Math.NumberTheory.ArithmeticFunctions+                          Math.NumberTheory.ArithmeticFunctions.Class+                          Math.NumberTheory.ArithmeticFunctions.Standard+                          Math.NumberTheory.Logarithms                           Math.NumberTheory.Moduli                           Math.NumberTheory.MoebiusInversion                           Math.NumberTheory.MoebiusInversion.Int@@ -56,6 +59,7 @@                           Math.NumberTheory.GCD.LowLevel                           Math.NumberTheory.Powers                           Math.NumberTheory.Powers.Squares+                          Math.NumberTheory.Powers.Squares.Internal                           Math.NumberTheory.Powers.Cubes                           Math.NumberTheory.Powers.Fourth                           Math.NumberTheory.Powers.General@@ -68,9 +72,9 @@                           Math.NumberTheory.Primes.Testing                           Math.NumberTheory.Primes.Testing.Certificates                           Math.NumberTheory.Primes.Heap+                          Math.NumberTheory.UniqueFactorisation     other-modules       : Math.NumberTheory.Utils                           Math.NumberTheory.Unsafe-                          Math.NumberTheory.Logarithms.Internal                           Math.NumberTheory.Primes.Counting.Impl                           Math.NumberTheory.Primes.Counting.Approximate                           Math.NumberTheory.Primes.Factorisation.Montgomery@@ -96,8 +100,12 @@  benchmark criterion   build-depends:    base-                    ,arithmoi-                    ,criterion+                    , arithmoi+                    , criterion+                    , containers+                    , random+  other-modules:    Math.NumberTheory.ArithmeticFunctionsBench+                  , Math.NumberTheory.PowersBench   hs-source-dirs:   benchmark   main-is:          Bench.hs   type:             exitcode-stdio-1.0@@ -109,15 +117,22 @@   ghc-options:          -Wall   main-is:              Test.hs   default-language: Haskell2010-  build-depends:        base >= 4 && < 5-                        , arithmoi >= 0.4 && < 0.5-                        , tasty >= 0.10 && < 0.12-                        , tasty-smallcheck >= 0.8 && < 0.9-                        , tasty-quickcheck >= 0.8 && < 0.9-                        , tasty-hunit >= 0.9 && < 0.10-                        , QuickCheck >= 2.8 && < 2.9-                        , smallcheck >= 1.1 && < 1.2-  other-modules :   Math.NumberTheory.GaussianIntegersTests+  build-depends:        base >= 4.6 && < 5+                      , containers >= 0.5 && < 0.6+                      , arithmoi >= 0.4 && < 0.5+                      , tasty >= 0.10 && < 0.12+                      , tasty-smallcheck >= 0.8 && < 0.9+                      , tasty-quickcheck >= 0.8 && < 0.9+                      , tasty-hunit >= 0.9 && < 0.10+                      , QuickCheck >= 2.7.6 && < 2.10+                      , smallcheck >= 1.1 && < 1.2+                      , transformers >= 0.3+                      , transformers-compat >= 0.4+  if impl(ghc < 7.10)+    build-depends     : nats >= 1 && <1.2++  other-modules :   Math.NumberTheory.ArithmeticFunctionsTests+                  , Math.NumberTheory.GaussianIntegersTests                   , Math.NumberTheory.GCDTests                   , Math.NumberTheory.GCD.LowLevelTests                   , Math.NumberTheory.LogarithmsTests@@ -135,3 +150,6 @@                   , Math.NumberTheory.Primes.HeapTests                   , Math.NumberTheory.Primes.SieveTests                   , Math.NumberTheory.TestUtils+                  , Math.NumberTheory.TestUtils.Wrappers+                  , Math.NumberTheory.TestUtils.Compose+                  , Math.NumberTheory.UniqueFactorisationTests

@@ -1,3 +1,13 @@+{-# OPTIONS_GHC -fno-warn-warnings-deprecations #-}+ module Main where -main = return ()+import Criterion.Main++import Math.NumberTheory.ArithmeticFunctionsBench as ArithmeticFunctions+import Math.NumberTheory.PowersBench as Powers++main = defaultMain+  [ ArithmeticFunctions.benchSuite+  , Powers.benchSuite+  ]

@@ -0,0 +1,34 @@+{-# OPTIONS_GHC -fno-warn-warnings-deprecations #-}++module Math.NumberTheory.ArithmeticFunctionsBench+  ( benchSuite+  ) where++import Criterion.Main+import Data.Set (Set)++import Math.NumberTheory.ArithmeticFunctions as A+import Math.NumberTheory.Primes.Factorisation as F++compareFunctions :: String -> (Integer -> Integer) -> (Integer -> Integer) -> Benchmark+compareFunctions name old new = bgroup name+  [ bench "old" $ nf (map old) [1..100000]+  , bench "new" $ nf (map new) [1..100000]+  ]++compareSetFunctions :: String -> (Integer -> Set Integer) -> (Integer -> Set Integer) -> Benchmark+compareSetFunctions name old new = bgroup name+  [ bench "old" $ nf (map old) [1..100000]+  , bench "new" $ nf (map new) [1..100000]+  ]++benchSuite = bgroup "ArithmeticFunctions"+  [ compareSetFunctions "divisors" F.divisors A.divisors+  , bench "divisors/int" $ nf (map A.divisorsSmall) [1 :: Int .. 100000]+  , compareFunctions "totient" F.totient A.totient+  , compareFunctions "carmichael" F.carmichael A.carmichael+  , compareFunctions "moebius" F.moebius A.moebius+  , compareFunctions "tau" F.tau A.tau+  , compareFunctions "sigma 1" (F.sigma 1) (A.sigma 1)+  , compareFunctions "sigma 2" (F.sigma 2) (A.sigma 2)+  ]

@@ -0,0 +1,30 @@+{-# OPTIONS_GHC -fno-warn-warnings-deprecations #-}++module Math.NumberTheory.PowersBench+  ( benchSuite+  ) where++import Criterion.Main+import System.Random++import Math.NumberTheory.Logarithms (integerLog2)+import Math.NumberTheory.Powers.Squares.Internal++genInteger :: Int -> Int -> Integer+genInteger salt bits+    = head+    . dropWhile ((< bits) . integerLog2)+    . scanl (\a r -> a * 2^31 + abs r) 1+    . randoms+    . mkStdGen+    $ salt + bits++compareRoots :: Int -> Benchmark+compareRoots bits = bgroup ("sqrt" ++ show bits)+  [ bench "new" $ nf (fst . karatsubaSqrt) n+  , bench "old" $ nf isqrtA n+  ]+  where+    n = genInteger 0 bits++benchSuite = bgroup "Powers" $ map compareRoots [2300, 2400 .. 2600]

@@ -0,0 +1,279 @@+-- |+-- Module:      Math.NumberTheory.ArithmeticFunctionsTests+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.ArithmeticFunctions+--++{-# LANGUAGE CPP       #-}++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.ArithmeticFunctionsTests+  ( testSuite+  ) where++import Test.Tasty+import Test.Tasty.HUnit++#if MIN_VERSION_base(4,8,0)+#else+import Prelude hiding (sum, all, elem)+import Data.Foldable+#endif++import qualified Data.Set as S+import qualified Data.IntSet as IS++import Math.NumberTheory.ArithmeticFunctions+import Math.NumberTheory.Primes.Factorisation+import Math.NumberTheory.TestUtils++import Numeric.Natural++oeisAssertion :: (Eq a, Show a) => String -> ArithmeticFunction Natural a -> [a] -> Assertion+oeisAssertion name f baseline = assertEqual name baseline (map (runFunction f) [1 .. fromIntegral (length baseline)])++-- | tau(n) equals to a number of divisors.+divisorsProperty1 :: Natural -> Bool+divisorsProperty1 n = S.size (runFunction divisorsA n) == runFunction tauA n++-- | sigma(n) equals to a number of divisors.+divisorsProperty2 :: Natural -> Bool+divisorsProperty2 n = sum (runFunction divisorsA n) == runFunction (sigmaA 1) n++-- | All divisors of n truly divides n.+divisorsProperty3 :: Natural -> Bool+divisorsProperty3 n = all (\d -> n `mod` d == 0) (runFunction divisorsA n)++-- | All divisors of n truly divides n.+divisorsProperty4 :: Int -> Bool+divisorsProperty4 n = S.toAscList (runFunction divisorsA n) == IS.toAscList (runFunction divisorsSmallA n)++-- | tau matches baseline from OEIS.+tauOeis :: Assertion+tauOeis = oeisAssertion "A000005" tauA+  [ 1, 2, 2, 3, 2, 4, 2, 4, 3, 4, 2, 6, 2, 4, 4, 5, 2, 6, 2, 6, 4, 4, 2, 8+  , 3, 4, 4, 6, 2, 8, 2, 6, 4, 4, 4, 9, 2, 4, 4, 8, 2, 8, 2, 6, 6, 4, 2, 10+  , 3, 6, 4, 6, 2, 8, 4, 8, 4, 4, 2, 12, 2, 4, 6, 7, 4, 8, 2, 6, 4, 8, 2+  , 12, 2, 4, 6, 6, 4, 8, 2, 10, 5, 4, 2, 12, 4, 4, 4, 8, 2, 12, 4, 6, 4, 4+  , 4, 12, 2, 6, 6, 9, 2, 8, 2, 8+  ]++-- | sigma_0 coincides with tau by definition+sigmaProperty1 :: Natural -> Bool+sigmaProperty1 n = runFunction tauA n == runFunction (sigmaA 0) n++-- | value of totient is bigger than argument+sigmaProperty2 :: Natural -> Bool+sigmaProperty2 n = n <= 1 || runFunction (sigmaA 1) n > n++-- | sigma_1 matches baseline from OEIS.+sigma1Oeis :: Assertion+sigma1Oeis = oeisAssertion "A000203" (sigmaA 1)+  [ 1, 3, 4, 7, 6, 12, 8, 15, 13, 18, 12, 28, 14, 24, 24, 31, 18, 39, 20+  , 42, 32, 36, 24, 60, 31, 42, 40, 56, 30, 72, 32, 63, 48, 54, 48, 91, 38+  , 60, 56, 90, 42, 96, 44, 84, 78, 72, 48, 124, 57, 93, 72, 98, 54, 120+  , 72, 120, 80, 90, 60, 168, 62, 96, 104, 127, 84, 144, 68, 126, 96, 144+  ]++-- | sigma_2 matches baseline from OEIS.+sigma2Oeis :: Assertion+sigma2Oeis = oeisAssertion "A001157" (sigmaA 2)+  [ 1, 5, 10, 21, 26, 50, 50, 85, 91, 130, 122, 210, 170, 250, 260, 341, 290+  , 455, 362, 546, 500, 610, 530, 850, 651, 850, 820, 1050, 842, 1300, 962+  , 1365, 1220, 1450, 1300, 1911, 1370, 1810, 1700, 2210, 1682, 2500, 1850+  , 2562, 2366, 2650, 2210, 3410, 2451, 3255+  ]++-- | value of totient if even, except totient(1) and totient(2)+totientProperty1 :: Natural -> Bool+totientProperty1 n = n <= 2 || even (runFunction totientA n)++-- | value of totient is smaller than argument+totientProperty2 :: Natural -> Bool+totientProperty2 n = n <= 1 || runFunction totientA n < n++totientSieve100 :: TotientSieve+totientSieve100 = totientSieve 100++-- | totient matches sieveTotient+totientProperty3 :: Natural -> Bool+totientProperty3 n = n < 1+  || fromIntegral (runFunction totientA n)+    == sieveTotient totientSieve100 (fromIntegral n)++-- | totient matches baseline from OEIS.+totientOeis :: Assertion+totientOeis = oeisAssertion "A000010" totientA+  [ 1, 1, 2, 2, 4, 2, 6, 4, 6, 4, 10, 4, 12, 6, 8, 8, 16, 6, 18, 8, 12, 10+  , 22, 8, 20, 12, 18, 12, 28, 8, 30, 16, 20, 16, 24, 12, 36, 18, 24, 16, 40+  , 12, 42, 20, 24, 22, 46, 16, 42, 20, 32, 24, 52, 18, 40, 24, 36, 28, 58+  , 16, 60, 30, 36, 32, 48, 20, 66, 32, 44+  ]++-- | jordan_0 is zero for argument > 1+jordanProperty1 :: Natural -> Bool+jordanProperty1 n = n <= 1 || runFunction (jordanA 0) n == 0++-- | jordan_1 coincides with totient by definition+jordanProperty2 :: Natural -> Bool+jordanProperty2 n = runFunction totientA n == runFunction (jordanA 1) n++-- | jordan_2 matches baseline from OEIS.+jordan2Oeis :: Assertion+jordan2Oeis = oeisAssertion "A007434" (jordanA 2)+  [ 1, 3, 8, 12, 24, 24, 48, 48, 72, 72, 120, 96, 168, 144, 192, 192, 288+  , 216, 360, 288, 384, 360, 528, 384, 600, 504, 648, 576, 840, 576, 960+  , 768, 960, 864, 1152, 864, 1368, 1080, 1344, 1152, 1680, 1152, 1848, 1440+  , 1728, 1584, 2208, 1536+  ]++-- | moebius values are [-1, 0, 1]+moebiusProperty1 :: Natural -> Bool+moebiusProperty1 n = runFunction moebiusA n `elem` [-1, 0, 1]++-- | moebius does not require full factorisation+moebiusLazy :: Assertion+moebiusLazy = assertEqual "moebius" 0 (runFunction moebiusA (2^2 * (2^100000-1) :: Natural))++-- | moebius matches baseline from OEIS.+moebiusOeis :: Assertion+moebiusOeis = oeisAssertion "A008683" moebiusA+  [ 1, -1, -1, 0, -1, 1, -1, 0, 0, 1, -1, 0, -1, 1, 1, 0, -1, 0, -1, 0, 1, 1, -1+  , 0, 0, 1, 0, 0, -1, -1, -1, 0, 1, 1, 1, 0, -1, 1, 1, 0, -1, -1, -1, 0, 0, 1+  , -1, 0, 0, 0, 1, 0, -1, 0, 1, 0, 1, 1, -1, 0, -1, 1, 0, 0, 1, -1, -1, 0, 1+  , -1, -1, 0, -1, 1, 0, 0, 1+  ]++-- | liouville values are [-1, 1]+liouvilleProperty1 :: Natural -> Bool+liouvilleProperty1 n = runFunction liouvilleA n `elem` [-1, 1]++-- | moebius is zero or equal to liouville+liouvilleProperty2 :: Natural -> Bool+liouvilleProperty2 n = m == 0 || l == m+  where+    l = runFunction liouvilleA n+    m = runFunction moebiusA   n++-- | liouville matches baseline from OEIS.+liouvilleOeis :: Assertion+liouvilleOeis = oeisAssertion "A008836" liouvilleA+  [ 1, -1, -1, 1, -1, 1, -1, -1, 1, 1, -1, -1, -1, 1, 1, 1, -1, -1, -1, -1, 1, 1+  , -1, 1, 1, 1, -1, -1, -1, -1, -1, -1, 1, 1, 1, 1, -1, 1, 1, 1, -1, -1, -1, -1+  , -1, 1, -1, -1, 1, -1, 1, -1, -1, 1, 1, 1, 1, 1, -1, 1, -1, 1, -1, 1, 1, -1+  , -1, -1, 1, -1, -1, -1, -1, 1, -1, -1, 1, -1, -1, -1, 1, 1, -1, 1, 1, 1, 1, 1+  , -1, 1, 1, -1, 1, 1, 1, 1, -1, -1, -1, 1, -1+  ]++-- | carmichaeil divides totient+carmichaelProperty1 :: Natural -> Bool+carmichaelProperty1 n = runFunction totientA n `mod` runFunction carmichaelA n == 0++carmichaelSieve100 :: CarmichaelSieve+carmichaelSieve100 = carmichaelSieve 100++-- | carmichael matches sieveCarmichael+carmichaelProperty2 :: Natural -> Bool+carmichaelProperty2 n = n < 1+  || fromIntegral (runFunction carmichaelA n)+    == sieveCarmichael carmichaelSieve100 (fromIntegral n)++-- | carmichael matches baseline from OEIS.+carmichaelOeis :: Assertion+carmichaelOeis = oeisAssertion "A002322" carmichaelA+  [ 1, 1, 2, 2, 4, 2, 6, 2, 6, 4, 10, 2, 12, 6, 4, 4, 16, 6, 18, 4, 6, 10, 22, 2+  , 20, 12, 18, 6, 28, 4, 30, 8, 10, 16, 12, 6, 36, 18, 12, 4, 40, 6, 42, 10, 12+  , 22, 46, 4, 42, 20, 16, 12, 52, 18, 20, 6, 18, 28, 58, 4, 60, 30, 6, 16, 12+  , 10, 66, 16, 22, 12, 70, 6, 72, 36, 20, 18, 30, 12, 78, 4, 54+  ]++-- | smallOmega is smaller than bigOmega+omegaProperty1 :: Natural -> Bool+omegaProperty1 n = runFunction smallOmegaA n <= runFunction bigOmegaA n++-- | smallOmega matches baseline from OEIS.+smallOmegaOeis :: Assertion+smallOmegaOeis = oeisAssertion "A001221" smallOmegaA+  [ 0, 1, 1, 1, 1, 2, 1, 1, 1, 2, 1, 2, 1, 2, 2, 1, 1, 2, 1, 2, 2, 2, 1, 2, 1, 2+  , 1, 2, 1, 3, 1, 1, 2, 2, 2, 2, 1, 2, 2, 2, 1, 3, 1, 2, 2, 2, 1, 2, 1, 2, 2, 2+  , 1, 2, 2, 2, 2, 2, 1, 3, 1, 2, 2, 1, 2, 3, 1, 2, 2, 3, 1, 2, 1, 2, 2, 2, 2, 3+  , 1, 2, 1, 2, 1, 3, 2, 2, 2, 2, 1, 3, 2, 2, 2, 2, 2, 2, 1, 2, 2, 2, 1, 3, 1, 2+  , 3, 2, 1, 2, 1, 3, 2+  ]++-- | bigOmega matches baseline from OEIS.+bigOmegaOeis :: Assertion+bigOmegaOeis = oeisAssertion "A001222" bigOmegaA+  [ 0, 1, 1, 2, 1, 2, 1, 3, 2, 2, 1, 3, 1, 2, 2, 4, 1, 3, 1, 3, 2, 2, 1, 4, 2, 2+  , 3, 3, 1, 3, 1, 5, 2, 2, 2, 4, 1, 2, 2, 4, 1, 3, 1, 3, 3, 2, 1, 5, 2, 3, 2, 3+  , 1, 4, 2, 4, 2, 2, 1, 4, 1, 2, 3, 6, 2, 3, 1, 3, 2, 3, 1, 5, 1, 2, 3, 3, 2, 3+  , 1, 5, 4, 2, 1, 4, 2, 2, 2, 4, 1, 4, 2, 3, 2, 2, 2, 6, 1, 3, 3, 4, 1, 3, 1, 4+  , 3, 2, 1, 5, 1, 3, 2+  ]++-- | expMangoldt matches baseline from OEIS.+mangoldtOeis :: Assertion+mangoldtOeis = oeisAssertion "A014963" expMangoldtA+  [ 1, 2, 3, 2, 5, 1, 7, 2, 3, 1, 11, 1, 13, 1, 1, 2, 17, 1, 19, 1, 1, 1, 23, 1+  , 5, 1, 3, 1, 29, 1, 31, 2, 1, 1, 1, 1, 37, 1, 1, 1, 41, 1, 43, 1, 1, 1, 47, 1+  , 7, 1, 1, 1, 53, 1, 1, 1, 1, 1, 59, 1, 61, 1, 1, 2, 1, 1, 67, 1, 1, 1, 71, 1+  , 73, 1, 1, 1, 1, 1, 79, 1, 3, 1, 83, 1, 1, 1, 1, 1, 89, 1, 1, 1, 1, 1, 1+  ]++testSuite :: TestTree+testSuite = testGroup "ArithmeticFunctions"+  [ testGroup "Divisors"+    [ testSmallAndQuick "length . divisors = tau"  divisorsProperty1+    , testSmallAndQuick "sum . divisors = sigma_1" divisorsProperty2+    , testSmallAndQuick "matches definition"       divisorsProperty3+    , testSmallAndQuick "divisors = divisorsSmall" divisorsProperty4+    ]+  , testGroup "Tau"+    [ testCase "OEIS" tauOeis+    ]+  , testGroup "Sigma"+    [ testSmallAndQuick "sigma_0 = tau" sigmaProperty1+    , testSmallAndQuick "sigma_1 n > n" sigmaProperty2+    , testCase          "OEIS sigma_1"  sigma1Oeis+    , testCase          "OEIS sigma_2"  sigma2Oeis+    ]+  , testGroup "Totient"+    [ testSmallAndQuick "totient is even"      totientProperty1+    , testSmallAndQuick "totient n < n"        totientProperty2+    , testSmallAndQuick "matches sieveTotient" totientProperty3+    , testCase          "OEIS"                 totientOeis+    ]+  , testGroup "Jordan"+    [ testSmallAndQuick "jordan_0 = [== 1]"  jordanProperty1+    , testSmallAndQuick "jordan_1 = totient" jordanProperty2+    , testCase          "OEIS jordan_2"      jordan2Oeis+    ]+  , testGroup "Moebius"+    [ testSmallAndQuick "moebius values" moebiusProperty1+    , testCase          "OEIS"           moebiusOeis+    , testCase          "Lazy"           moebiusLazy+    ]+  , testGroup "Liouville"+    [ testSmallAndQuick "liouville values"          liouvilleProperty1+    , testSmallAndQuick "liouville matches moebius" liouvilleProperty2+    , testCase          "OEIS"                      liouvilleOeis+    ]+  , testGroup "Carmichael"+    [ testSmallAndQuick "carmichael divides totient" carmichaelProperty1+    , testSmallAndQuick "matches sieveCarmichael"    carmichaelProperty2+    , testCase          "OEIS"                       carmichaelOeis+    ]+  , testGroup "Omegas"+    [ testSmallAndQuick "smallOmega <= bigOmega" omegaProperty1+    , testCase          "OEIS smallOmega"        smallOmegaOeis+    , testCase          "OEIS bigOmega"          bigOmegaOeis+    ]+  , testGroup "Mangoldt"+    [ testCase "OEIS" mangoldtOeis+    ]+  ]

@@ -8,8 +8,8 @@ -- Tests for Math.NumberTheory.Moduli -- -{-# LANGUAGE CPP          #-}-{-# LANGUAGE ViewPatterns #-}+{-# LANGUAGE CPP             #-}+{-# LANGUAGE ViewPatterns    #-}  {-# OPTIONS_GHC -fno-warn-type-defaults #-} @@ -23,51 +23,51 @@ import Data.Bits import Data.List (tails, nub) import Data.Maybe+import Data.Functor.Compose  import Math.NumberTheory.Moduli import Math.NumberTheory.TestUtils -toOdd :: Num a => a -> a-toOdd n = n * 2 + 1- unwrapPP :: (Prime, Power Int) -> (Integer, Int) unwrapPP (Prime p, Power e) = (p, e)  -- | Check that 'jacobi' matches 'jacobi''.-jacobiProperty1 :: (Integral a, Bits a) => AnySign a -> NonNegative a -> Bool-jacobiProperty1 (AnySign a) (NonNegative (toOdd -> n)) = n == 1 && j == 1 || n > 1 && j == j'+jacobiProperty1 :: (Integral a, Bits a) => AnySign a -> (Compose Positive Odd) a -> Bool+jacobiProperty1 (AnySign a) (Compose (Positive (Odd n))) = n == 1 && j == 1 || n > 1 && j == j'   where     j = jacobi a n     j' = jacobi' a n  -- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 2-jacobiProperty2 :: (Integral a, Bits a) => AnySign a -> NonNegative a -> Bool-jacobiProperty2 (AnySign a) (NonNegative (toOdd -> n)) = jacobi a n == jacobi (a + n) n+jacobiProperty2 :: (Integral a, Bits a) => AnySign a -> (Compose Positive Odd) a -> Bool+jacobiProperty2 (AnySign a) (Compose (Positive (Odd n)))+  =  a + n < a -- check overflow+  || jacobi a n == jacobi (a + n) n  -- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 3-jacobiProperty3 :: (Integral a, Bits a) => AnySign a -> NonNegative a -> Bool-jacobiProperty3 (AnySign a) (NonNegative (toOdd -> n)) = j == 0 && g /= 1 || abs j == 1 && g == 1+jacobiProperty3 :: (Integral a, Bits a) => AnySign a -> (Compose Positive Odd) a -> Bool+jacobiProperty3 (AnySign a) (Compose (Positive (Odd n))) = j == 0 && g /= 1 || abs j == 1 && g == 1   where     j = jacobi a n     g = gcd a n  -- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 4-jacobiProperty4 :: (Integral a, Bits a) => AnySign a -> AnySign a -> NonNegative a -> Bool-jacobiProperty4 (AnySign a) (AnySign b) (NonNegative (toOdd -> n)) = jacobi (a * b) n == jacobi a n * jacobi b n+jacobiProperty4 :: (Integral a, Bits a) => AnySign a -> AnySign a -> (Compose Positive Odd) a -> Bool+jacobiProperty4 (AnySign a) (AnySign b) (Compose (Positive (Odd n))) = jacobi (a * b) n == jacobi a n * jacobi b n -jacobiProperty4_Integer :: AnySign Integer -> AnySign Integer -> NonNegative Integer -> Bool+jacobiProperty4_Integer :: AnySign Integer -> AnySign Integer -> (Compose Positive Odd) Integer -> Bool jacobiProperty4_Integer = jacobiProperty4  -- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 5-jacobiProperty5 :: (Integral a, Bits a) => AnySign a -> NonNegative a -> NonNegative a -> Bool-jacobiProperty5 (AnySign a) (NonNegative (toOdd -> m)) (NonNegative (toOdd -> n)) = jacobi a (m * n) == jacobi a m * jacobi a n+jacobiProperty5 :: (Integral a, Bits a) => AnySign a -> (Compose Positive Odd) a -> (Compose Positive Odd) a -> Bool+jacobiProperty5 (AnySign a) (Compose (Positive (Odd m))) (Compose (Positive (Odd n))) = jacobi a (m * n) == jacobi a m * jacobi a n -jacobiProperty5_Integer :: AnySign Integer -> NonNegative Integer -> NonNegative Integer -> Bool+jacobiProperty5_Integer :: AnySign Integer -> (Compose Positive Odd) Integer -> (Compose Positive Odd) Integer -> Bool jacobiProperty5_Integer = jacobiProperty5  -- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 6-jacobiProperty6 :: (Integral a, Bits a) => NonNegative a -> NonNegative a -> Bool-jacobiProperty6 (NonNegative (toOdd -> m)) (NonNegative (toOdd -> n)) = gcd m n /= 1 || jacobi m n * jacobi n m == (if m `mod` 4 == 1 || n `mod` 4 == 1 then 1 else -1)+jacobiProperty6 :: (Integral a, Bits a) => (Compose Positive Odd) a -> (Compose Positive Odd) a -> Bool+jacobiProperty6 (Compose (Positive (Odd m))) (Compose (Positive (Odd n))) = gcd m n /= 1 || jacobi m n * jacobi n m == (if m `mod` 4 == 1 || n `mod` 4 == 1 then 1 else -1)  -- | Check that 'invertMod' inverts numbers modulo. invertModProperty :: AnySign Integer -> Positive Integer -> Bool@@ -98,6 +98,10 @@ powerModProperty3 :: (Integral a, Bits a) => AnySign a -> AnySign a -> AnySign Integer -> Positive Integer -> Bool powerModProperty3 (AnySign e1) (AnySign e2) (AnySign b) (Positive m)   =  (e1 < 0 || e2 < 0) && isNothing (invertMod b m)+  || e2 >= 0 && e1 + e2 < e1 -- check overflow+  || e1 >= 0 && e1 + e2 < e2 -- check overflow+  || e2 <= 0 && e1 + e2 > e1 -- check overflow+  || e1 <= 0 && e1 + e2 > e2 -- check overflow   || pm1 * pm2 `mod` m == pm12   where     pm1  = powerMod b e1 m

@@ -19,7 +19,7 @@ import Test.Tasty.QuickCheck as QC hiding (Positive)  import Math.NumberTheory.MoebiusInversion.Int-import Math.NumberTheory.Primes.Factorisation+import Math.NumberTheory.ArithmeticFunctions import Math.NumberTheory.TestUtils  totientSumProperty :: Positive Int -> Bool

@@ -19,7 +19,7 @@ import Test.Tasty.QuickCheck as QC hiding (Positive)  import Math.NumberTheory.MoebiusInversion-import Math.NumberTheory.Primes.Factorisation+import Math.NumberTheory.ArithmeticFunctions import Math.NumberTheory.TestUtils  totientSumProperty :: Positive Int -> Bool

@@ -49,6 +49,10 @@ integerSquareRootProperty_Word :: NonNegative Word -> Bool integerSquareRootProperty_Word = integerSquareRootProperty +-- | Specialized to trigger 'isqrtInteger'.+integerSquareRootProperty_Integer :: NonNegative Integer -> Bool+integerSquareRootProperty_Integer = integerSquareRootProperty+ -- | Check that 'integerSquareRoot' returns the largest integer @m@ with @m*m <= n@, where @n@ has form @k@^2-1. integerSquareRootProperty2 :: Integral a => NonNegative a -> Bool integerSquareRootProperty2 (NonNegative k) = n < 0@@ -65,6 +69,10 @@ integerSquareRootProperty2_Word :: NonNegative Word -> Bool integerSquareRootProperty2_Word = integerSquareRootProperty2 +-- | Specialized to trigger 'isqrtInteger'.+integerSquareRootProperty2_Integer :: NonNegative Integer -> Bool+integerSquareRootProperty2_Integer = integerSquareRootProperty2+ #if WORD_SIZE_IN_BITS == 64  -- | Check that 'integerSquareRoot' of 2^62-1 is 2^31-1, not 2^31.@@ -130,14 +138,16 @@  testSuite :: TestTree testSuite = testGroup "Squares"-  [ testGroup "integerSquareRoor"+  [ testGroup "integerSquareRoot"     [ testIntegralProperty "generic"          integerSquareRootProperty     , testSmallAndQuick    "generic Int"      integerSquareRootProperty_Int     , testSmallAndQuick    "generic Word"     integerSquareRootProperty_Word+    , testSmallAndQuick    "generic Integer"  integerSquareRootProperty_Integer -    , testIntegralProperty "almost square"      integerSquareRootProperty2-    , testSmallAndQuick    "almost square Int"  integerSquareRootProperty2_Int-    , testSmallAndQuick    "almost square Word" integerSquareRootProperty2_Word+    , testIntegralProperty "almost square"         integerSquareRootProperty2+    , testSmallAndQuick    "almost square Int"     integerSquareRootProperty2_Int+    , testSmallAndQuick    "almost square Word"    integerSquareRootProperty2_Word+    , testSmallAndQuick    "almost square Integer" integerSquareRootProperty2_Integer  #if WORD_SIZE_IN_BITS == 64     , testCase             "maxBound / 2 :: Int"  integerSquareRootSpecialCase1_Int

@@ -31,10 +31,16 @@   sieveFactorSpecialCase1 :: Assertion-sieveFactorSpecialCase1 = assertEqual "sieveFactor" [(29, 1), (73, 1)] $ sieveFactor (factorSieve 2048) (29*73)+sieveFactorSpecialCase1 = do+    assertEqual "sieveFactor 2048" [(29, 1), (73, 1)] $ sieveFactor (factorSieve 2048) (29*73)+    assertEqual "sieveFactor 15"   [(3, 1),  (5, 2)]  $ sieveFactor (factorSieve 15) (75) +sieveFactorProperty :: Integer -> Bool+sieveFactorProperty n = (n==0) || factorise n == sieveFactor (factorSieve 25) n+ testSuite :: TestTree testSuite = testGroup "Primes"   [ testSmallAndQuick "primesSum"   primesSumProperty-  , testCase          "sieveFactor" sieveFactorSpecialCase1+  , testCase          "sieveFactor special"  sieveFactorSpecialCase1+  , testSmallAndQuick "sieveFactor property" sieveFactorProperty   ]

@@ -12,12 +12,8 @@ {-# LANGUAGE ConstraintKinds            #-} {-# LANGUAGE CPP                        #-} {-# LANGUAGE DataKinds                  #-}-{-# LANGUAGE DeriveFoldable             #-}-{-# LANGUAGE DeriveFunctor              #-}-{-# LANGUAGE DeriveTraversable          #-} {-# LANGUAGE FlexibleContexts           #-} {-# LANGUAGE FlexibleInstances          #-}-{-# LANGUAGE GeneralizedNewtypeDeriving #-} {-# LANGUAGE KindSignatures             #-} {-# LANGUAGE MultiParamTypeClasses      #-} {-# LANGUAGE RankNTypes                 #-}@@ -36,11 +32,16 @@ {-# OPTIONS_GHC -fno-warn-type-defaults #-}  module Math.NumberTheory.TestUtils-  ( module Math.NumberTheory.TestUtils+  ( module Math.NumberTheory.TestUtils.Wrappers   , module Test.SmallCheck.Series   , Large(..)+  , testIntegralProperty+  , testSameIntegralProperty+  , testIntegral2Property+  , testSmallAndQuick   ) where +import Test.SmallCheck.Series (cons2) import Test.Tasty import Test.Tasty.SmallCheck as SC import Test.Tasty.QuickCheck as QC hiding (Positive, NonNegative, generate, getNonNegative)@@ -49,75 +50,41 @@  import Control.Applicative import Data.Bits-#if MIN_VERSION_base(4,8,0)-#else-import Data.Foldable (Foldable)-import Data.Traversable (Traversable)+#if !(MIN_VERSION_base(4,8,0)) import Data.Word #endif import GHC.Exts--import Math.NumberTheory.Primes--newtype AnySign a = AnySign { getAnySign :: a }-  deriving (Eq, Ord, Read, Show, Num, Enum, Bounded, Integral, Real, Functor, Foldable, Traversable, Arbitrary)--instance (Monad m, Serial m a) => Serial m (AnySign a) where-  series = AnySign <$> series--instance (Num a, Ord a, Arbitrary a) => Arbitrary (Positive a) where-  arbitrary = Positive <$> (arbitrary `suchThat` (> 0))-  shrink (Positive x) = Positive <$> filter (> 0) (shrink x)--instance (Num a, Ord a, Arbitrary a) => Arbitrary (NonNegative a) where-  arbitrary = NonNegative <$> (arbitrary `suchThat` (>= 0))-  shrink (NonNegative x) = NonNegative <$> filter (>= 0) (shrink x)--instance (Num a, Bounded a) => Bounded (Positive a) where-  minBound = Positive 1-  maxBound = Positive (maxBound :: a)--instance (Num a, Bounded a) => Bounded (NonNegative a) where-  minBound = NonNegative 0-  maxBound = NonNegative (maxBound :: a)--newtype Huge a = Huge { getHuge :: a }-  deriving (Eq, Ord, Enum, Bounded, Show, Num, Real, Integral)--instance (Num a, Arbitrary a) => Arbitrary (Huge a) where-  arbitrary = do-    Positive l <- arbitrary-    ds <- vector l-    return $ Huge $ foldl1 (\acc n -> acc * 2^63 + n) ds--newtype Power a = Power { getPower :: a }-  deriving (Eq, Ord, Enum, Bounded, Show, Num, Real, Integral)--instance (Monad m, Num a, Ord a, Serial m a) => Serial m (Power a) where-  series = Power <$> series `suchThatSerial` (> 0)+import Numeric.Natural -instance (Num a, Ord a, Integral a, Arbitrary a) => Arbitrary (Power a) where-  arbitrary = Power <$> (getSmall <$> arbitrary) `suchThat` (> 0)-  shrink (Power x) = Power <$> filter (> 0) (shrink x)+import Math.NumberTheory.GaussianIntegers (GaussianInteger(..)) -newtype Prime = Prime { getPrime :: Integer }-  deriving (Eq, Ord, Show)+import Math.NumberTheory.TestUtils.Compose ()+import Math.NumberTheory.TestUtils.Wrappers -instance Arbitrary Prime where-  arbitrary = Prime <$> arbitrary `suchThat` (\p -> p > 0 && isPrime p)+instance Monad m => Serial m Word where+  series =+    generate (\d -> if d >= 0 then pure 0 else empty) <|> nats+    where+      nats = generate $ \d -> if d > 0 then [1 .. fromInteger (toInteger d)] else empty -instance Monad m => Serial m Prime where-  series = Prime <$> series `suchThatSerial` (\p -> p > 0 && isPrime p)+#if !(MIN_VERSION_base(4,8,0)) && !(MIN_VERSION_QuickCheck(2,9,0))+instance Arbitrary Natural where+  arbitrary = fromInteger <$> (arbitrary `suchThat` (>= 0))+  shrink = map fromInteger . filter (>= 0) . shrink . toInteger+#endif -instance Monad m => Serial m Word where+instance Monad m => Serial m Natural where   series =     generate (\d -> if d >= 0 then pure 0 else empty) <|> nats     where       nats = generate $ \d -> if d > 0 then [1 .. fromInteger (toInteger d)] else empty -suchThatSerial :: Series m a -> (a -> Bool) -> Series m a-suchThatSerial s p = s >>= \x -> if p x then pure x else empty+instance Arbitrary GaussianInteger where+  arbitrary = (:+) <$> arbitrary <*> arbitrary+  shrink (x :+ y) = (:+) <$> shrink x <*> shrink y +instance Monad m => Serial m GaussianInteger where+  series = cons2 (:+)  -- https://www.cs.ox.ac.uk/projects/utgp/school/andres.pdf, p. 21 -- :k Compose = (k1 -> Constraint) -> (k2 -> k1) -> (k2 -> Constraint)@@ -146,8 +113,10 @@  type TestableIntegral wrapper =   ( Matrix '[Arbitrary, Show, Serial IO] wrapper '[Int, Word, Integer]+  , Matrix '[Arbitrary, Show] wrapper '[Large Int, Large Word, Huge Integer]   , Matrix '[Bounded, Integral] wrapper '[Int, Word]   , Num (wrapper Integer)+  , Functor wrapper   )  @@ -176,9 +145,9 @@   , QC.testProperty "quickcheck Int"     (f :: wrapper1 Int     -> wrapper2 Int     -> bool)   , QC.testProperty "quickcheck Word"    (f :: wrapper1 Word    -> wrapper2 Word    -> bool)   , QC.testProperty "quickcheck Integer" (f :: wrapper1 Integer -> wrapper2 Integer -> bool)-  , QC.testProperty "quickcheck Large Int"     (\(Large a) (Large b) -> (f :: wrapper1 Int     -> wrapper2 Int     -> bool) a b)-  , QC.testProperty "quickcheck Large Word"    (\(Large a) (Large b) -> (f :: wrapper1 Word    -> wrapper2 Word    -> bool) a b)-  , QC.testProperty "quickcheck Huge  Integer" (\(Huge  a) (Huge  b) -> (f :: wrapper1 Integer -> wrapper2 Integer -> bool) a b)+  , QC.testProperty "quickcheck Large Int"     (\a b -> (f :: wrapper1 Int     -> wrapper2 Int     -> bool) (getLarge <$> a) (getLarge <$> b))+  , QC.testProperty "quickcheck Large Word"    (\a b -> (f :: wrapper1 Word    -> wrapper2 Word    -> bool) (getLarge <$> a) (getLarge <$> b))+  , QC.testProperty "quickcheck Huge  Integer" (\a b -> (f :: wrapper1 Integer -> wrapper2 Integer -> bool) (getHuge  <$> a) (getHuge  <$> b))   ]  testIntegral2Property@@ -205,15 +174,15 @@   , QC.testProperty "quickcheck Integer Word"    (f :: wrapper1 Integer -> wrapper2 Word    -> bool)   , QC.testProperty "quickcheck Integer Integer" (f :: wrapper1 Integer -> wrapper2 Integer -> bool) -  , QC.testProperty "quickcheck Large Int Int"         ((f :: wrapper1 Int     -> wrapper2 Int     -> bool) . getLarge)-  , QC.testProperty "quickcheck Large Int Word"        ((f :: wrapper1 Int     -> wrapper2 Word    -> bool) . getLarge)-  , QC.testProperty "quickcheck Large Int Integer"     ((f :: wrapper1 Int     -> wrapper2 Integer -> bool) . getLarge)-  , QC.testProperty "quickcheck Large Word Int"        ((f :: wrapper1 Word    -> wrapper2 Int     -> bool) . getLarge)-  , QC.testProperty "quickcheck Large Word Word"       ((f :: wrapper1 Word    -> wrapper2 Word    -> bool) . getLarge)-  , QC.testProperty "quickcheck Large Word Integer"    ((f :: wrapper1 Word    -> wrapper2 Integer -> bool) . getLarge)-  , QC.testProperty "quickcheck Huge  Integer Int"     ((f :: wrapper1 Integer -> wrapper2 Int     -> bool) . getHuge)-  , QC.testProperty "quickcheck Huge  Integer Word"    ((f :: wrapper1 Integer -> wrapper2 Word    -> bool) . getHuge)-  , QC.testProperty "quickcheck Huge  Integer Integer" ((f :: wrapper1 Integer -> wrapper2 Integer -> bool) . getHuge)+  , QC.testProperty "quickcheck Large Int Int"         ((f :: wrapper1 Int     -> wrapper2 Int     -> bool) . fmap getLarge)+  , QC.testProperty "quickcheck Large Int Word"        ((f :: wrapper1 Int     -> wrapper2 Word    -> bool) . fmap getLarge)+  , QC.testProperty "quickcheck Large Int Integer"     ((f :: wrapper1 Int     -> wrapper2 Integer -> bool) . fmap getLarge)+  , QC.testProperty "quickcheck Large Word Int"        ((f :: wrapper1 Word    -> wrapper2 Int     -> bool) . fmap getLarge)+  , QC.testProperty "quickcheck Large Word Word"       ((f :: wrapper1 Word    -> wrapper2 Word    -> bool) . fmap getLarge)+  , QC.testProperty "quickcheck Large Word Integer"    ((f :: wrapper1 Word    -> wrapper2 Integer -> bool) . fmap getLarge)+  , QC.testProperty "quickcheck Huge  Integer Int"     ((f :: wrapper1 Integer -> wrapper2 Int     -> bool) . fmap getHuge)+  , QC.testProperty "quickcheck Huge  Integer Word"    ((f :: wrapper1 Integer -> wrapper2 Word    -> bool) . fmap getHuge)+  , QC.testProperty "quickcheck Huge  Integer Integer" ((f :: wrapper1 Integer -> wrapper2 Integer -> bool) . fmap getHuge)   ]  testSmallAndQuick

@@ -0,0 +1,47 @@+-- |+-- Module:      Math.NumberTheory.TestUtils.Compose+-- 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 CPP                        #-}+{-# LANGUAGE DeriveGeneric              #-}+{-# LANGUAGE FlexibleContexts           #-}+{-# LANGUAGE FlexibleInstances          #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE MultiParamTypeClasses      #-}+{-# LANGUAGE StandaloneDeriving         #-}+{-# LANGUAGE UndecidableInstances       #-}++{-# OPTIONS_GHC -fno-warn-orphans       #-}++module Math.NumberTheory.TestUtils.Compose where++import Data.Functor.Compose+#if MIN_VERSION_transformers(0,5,0)+#else+import GHC.Generics+#endif++import Test.Tasty.QuickCheck (Arbitrary)+import Test.SmallCheck.Series (Serial)++deriving instance Num (f (g a))     => Num (Compose f g a)+deriving instance Enum (f (g a))    => Enum (Compose f g a)+deriving instance Bounded (f (g a)) => Bounded (Compose f g a)++deriving instance (Ord (Compose f g a), Real (f (g a)))     => Real (Compose f g a)+deriving instance (Ord (Compose f g a), Integral (f (g a))) => Integral (Compose f g a)++deriving instance Arbitrary (f (g a)) => Arbitrary (Compose f g a)++#if MIN_VERSION_transformers(0,5,0)+#else+deriving instance Generic (Compose f g a)+#endif+instance (Monad m, Serial m (f (g a))) => Serial m (Compose f g a)

@@ -0,0 +1,255 @@+-- |+-- Module:      Math.NumberTheory.TestUtils.Wrappers+-- 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 CPP                        #-}+{-# LANGUAGE DeriveFoldable             #-}+{-# LANGUAGE DeriveFunctor              #-}+{-# LANGUAGE DeriveTraversable          #-}+{-# LANGUAGE FlexibleInstances          #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE MultiParamTypeClasses      #-}+{-# LANGUAGE ScopedTypeVariables        #-}+{-# LANGUAGE StandaloneDeriving         #-}++{-# OPTIONS_GHC -fno-warn-orphans       #-}+{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.TestUtils.Wrappers where++import Control.Applicative+import Data.Functor.Classes+#if MIN_VERSION_base(4,8,0)+#else+import Data.Foldable (Foldable)+import Data.Traversable (Traversable)+#endif++import Test.Tasty.QuickCheck as QC hiding (Positive, NonNegative, generate, getNonNegative, getPositive)+import Test.SmallCheck.Series (Positive(..), NonNegative(..), Serial(..), Series)++import Math.NumberTheory.Primes (isPrime)++-------------------------------------------------------------------------------+-- AnySign++newtype AnySign a = AnySign { getAnySign :: a }+  deriving (Eq, Ord, Read, Show, Num, Enum, Bounded, Integral, Real, Functor, Foldable, Traversable, Arbitrary)++instance (Monad m, Serial m a) => Serial m (AnySign a) where+  series = AnySign <$> series++instance Eq1 AnySign where+#if MIN_VERSION_transformers(0,5,0)+  liftEq eq (AnySign a) (AnySign b) = a `eq` b+#else+  (AnySign a) `eq1` (AnySign b) = a == b+#endif++instance Ord1 AnySign where+#if MIN_VERSION_transformers(0,5,0)+  liftCompare cmp (AnySign a) (AnySign b) = a `cmp` b+#else+  (AnySign a) `compare1` (AnySign b) = a `compare` b+#endif++instance Show1 AnySign where+#if MIN_VERSION_transformers(0,5,0)+  liftShowsPrec shw _ p (AnySign a) = shw p a+#else+  showsPrec1 p (AnySign a) = showsPrec p a+#endif++-------------------------------------------------------------------------------+-- Positive from smallcheck++deriving instance Functor Positive++instance (Num a, Ord a, Arbitrary a) => Arbitrary (Positive a) where+  arbitrary = Positive <$> (arbitrary `suchThat` (> 0))+  shrink (Positive x) = Positive <$> filter (> 0) (shrink x)++instance (Num a, Bounded a) => Bounded (Positive a) where+  minBound = Positive 1+  maxBound = Positive (maxBound :: a)++instance Eq1 Positive where+#if MIN_VERSION_transformers(0,5,0)+  liftEq eq (Positive a) (Positive b) = a `eq` b+#else+  (Positive a) `eq1` (Positive b) = a == b+#endif++instance Ord1 Positive where+#if MIN_VERSION_transformers(0,5,0)+  liftCompare cmp (Positive a) (Positive b) = a `cmp` b+#else+  (Positive a) `compare1` (Positive b) = a `compare` b+#endif++instance Show1 Positive where+#if MIN_VERSION_transformers(0,5,0)+  liftShowsPrec shw _ p (Positive a) = shw p a+#else+  showsPrec1 p (Positive a) = showsPrec p a+#endif++-------------------------------------------------------------------------------+-- NonNegative from smallcheck++deriving instance Functor NonNegative++instance (Num a, Ord a, Arbitrary a) => Arbitrary (NonNegative a) where+  arbitrary = NonNegative <$> (arbitrary `suchThat` (>= 0))+  shrink (NonNegative x) = NonNegative <$> filter (>= 0) (shrink x)++instance (Num a, Bounded a) => Bounded (NonNegative a) where+  minBound = NonNegative 0+  maxBound = NonNegative (maxBound :: a)++instance Eq1 NonNegative where+#if MIN_VERSION_transformers(0,5,0)+  liftEq eq (NonNegative a) (NonNegative b) = a `eq` b+#else+  (NonNegative a) `eq1` (NonNegative b) = a == b+#endif++instance Ord1 NonNegative where+#if MIN_VERSION_transformers(0,5,0)+  liftCompare cmp (NonNegative a) (NonNegative b) = a `cmp` b+#else+  (NonNegative a) `compare1` (NonNegative b) = a `compare` b+#endif++instance Show1 NonNegative where+#if MIN_VERSION_transformers(0,5,0)+  liftShowsPrec shw _ p (NonNegative a) = shw p a+#else+  showsPrec1 p (NonNegative a) = showsPrec p a+#endif++-------------------------------------------------------------------------------+-- Huge++newtype Huge a = Huge { getHuge :: a }+  deriving (Eq, Ord, Read, Show, Num, Enum, Bounded, Integral, Real, Functor, Foldable, Traversable)++instance (Num a, Arbitrary a) => Arbitrary (Huge a) where+  arbitrary = do+    Positive l <- arbitrary+    ds <- vector l+    return $ Huge $ foldl1 (\acc n -> acc * 2^63 + n) ds++instance Eq1 Huge where+#if MIN_VERSION_transformers(0,5,0)+  liftEq eq (Huge a) (Huge b) = a `eq` b+#else+  (Huge a) `eq1` (Huge b) = a == b+#endif++instance Ord1 Huge where+#if MIN_VERSION_transformers(0,5,0)+  liftCompare cmp (Huge a) (Huge b) = a `cmp` b+#else+  (Huge a) `compare1` (Huge b) = a `compare` b+#endif++instance Show1 Huge where+#if MIN_VERSION_transformers(0,5,0)+  liftShowsPrec shw _ p (Huge a) = shw p a+#else+  showsPrec1 p (Huge a) = showsPrec p a+#endif++-------------------------------------------------------------------------------+-- Power++newtype Power a = Power { getPower :: a }+  deriving (Eq, Ord, Read, Show, Num, Enum, Bounded, Integral, Real, Functor, Foldable, Traversable)++instance (Monad m, Num a, Ord a, Serial m a) => Serial m (Power a) where+  series = Power <$> series `suchThatSerial` (> 0)++instance (Num a, Ord a, Integral a, Arbitrary a) => Arbitrary (Power a) where+  arbitrary = Power <$> (getSmall <$> arbitrary) `suchThat` (> 0)+  shrink (Power x) = Power <$> filter (> 0) (shrink x)++instance Eq1 Power where+#if MIN_VERSION_transformers(0,5,0)+  liftEq eq (Power a) (Power b) = a `eq` b+#else+  (Power a) `eq1` (Power b) = a == b+#endif++instance Ord1 Power where+#if MIN_VERSION_transformers(0,5,0)+  liftCompare cmp (Power a) (Power b) = a `cmp` b+#else+  (Power a) `compare1` (Power b) = a `compare` b+#endif++instance Show1 Power where+#if MIN_VERSION_transformers(0,5,0)+  liftShowsPrec shw _ p (Power a) = shw p a+#else+  showsPrec1 p (Power a) = showsPrec p a+#endif++-------------------------------------------------------------------------------+-- Odd++newtype Odd a = Odd { getOdd :: a }+  deriving (Eq, Ord, Read, Show, Num, Enum, Bounded, Integral, Real, Functor, Foldable, Traversable)++instance (Monad m, Serial m a, Integral a) => Serial m (Odd a) where+  series = Odd <$> series `suchThatSerial` odd++instance (Integral a, Arbitrary a) => Arbitrary (Odd a) where+  arbitrary = Odd <$> (arbitrary `suchThat` odd)+  shrink (Odd x) = Odd <$> filter odd (shrink x)++instance Eq1 Odd where+#if MIN_VERSION_transformers(0,5,0)+  liftEq eq (Odd a) (Odd b) = a `eq` b+#else+  (Odd a) `eq1` (Odd b) = a == b+#endif++instance Ord1 Odd where+#if MIN_VERSION_transformers(0,5,0)+  liftCompare cmp (Odd a) (Odd b) = a `cmp` b+#else+  (Odd a) `compare1` (Odd b) = a `compare` b+#endif++instance Show1 Odd where+#if MIN_VERSION_transformers(0,5,0)+  liftShowsPrec shw _ p (Odd a) = shw p a+#else+  showsPrec1 p (Odd a) = showsPrec p a+#endif++-------------------------------------------------------------------------------+-- Prime++newtype Prime = Prime { getPrime :: Integer }+  deriving (Eq, Ord, Show)++instance Arbitrary Prime where+  arbitrary = Prime <$> arbitrary `suchThat` (\p -> p > 0 && isPrime p)++instance Monad m => Serial m Prime where+  series = Prime <$> series `suchThatSerial` (\p -> p > 0 && isPrime p)++-------------------------------------------------------------------------------+-- Utils++suchThatSerial :: Series m a -> (a -> Bool) -> Series m a+suchThatSerial s p = s >>= \x -> if p x then pure x else empty

@@ -0,0 +1,52 @@+-- |+-- Module:      Math.NumberTheory.UniqueFactorisationTests+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.ArithmeticFunctions+--++{-# LANGUAGE CPP                 #-}+{-# LANGUAGE ScopedTypeVariables #-}++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.UniqueFactorisationTests+  ( testSuite+  ) where++import Test.Tasty++#if MIN_VERSION_base(4,8,0)+#else+import Data.Word+#endif++import Math.NumberTheory.GaussianIntegers hiding (factorise)+import Math.NumberTheory.UniqueFactorisation+import Math.NumberTheory.TestUtils hiding (Prime)++import Numeric.Natural++testRules :: forall a. (UniqueFactorisation a, Num a, Eq a) => a -> Bool+testRules n+  = n == 0+  || all (\(p, _) -> unP p == abs (unP p)) fs+  && abs n == abs (product (map (\(p, k) -> unP p ^ k) fs))+  where+    fs = factorise n++    unP :: Prime a -> a+    unP = unPrime++testSuite :: TestTree+testSuite = testGroup "UniqueFactorisation"+  [ testSmallAndQuick "Int"     (testRules :: Int     -> Bool)+  , testSmallAndQuick "Word"    (testRules :: Word    -> Bool)+  , testSmallAndQuick "Integer" (testRules :: Integer -> Bool)+  , testSmallAndQuick "Natural" (testRules :: Natural -> Bool)++  , testSmallAndQuick "GaussianInteger" (testRules :: GaussianInteger -> Bool)+  ]

@@ -25,6 +25,9 @@  import qualified Math.NumberTheory.GaussianIntegersTests as Gaussian +import qualified Math.NumberTheory.ArithmeticFunctionsTests as ArithmeticFunctions+import qualified Math.NumberTheory.UniqueFactorisationTests as UniqueFactorisation+ main :: IO () main = defaultMain tests @@ -62,5 +65,11 @@     ]   , testGroup "Gaussian"     [ Gaussian.testSuite+    ]+  , testGroup "ArithmeticFunctions"+    [ ArithmeticFunctions.testSuite+    ]+  , testGroup "UniqueFactorisation"+    [ UniqueFactorisation.testSuite     ]   ]