arithmoi 0.5.0.1 → 0.6.0.0

raw patch · 31 files changed

+1845/−957 lines, 31 filesdep ~arithmoidep ~basedep ~tastyPVP ok

version bump matches the API change (PVP)

Dependency ranges changed: arithmoi, base, tasty, tasty-hunit

API changes (from Hackage documentation)

Files

• Changes +35/−0
• GHC/TypeNats/Compat.hs +36/−0
• Math/NumberTheory/ArithmeticFunctions/Standard.hs +23/−0
• Math/NumberTheory/Curves/Montgomery.hs +172/−0
• Math/NumberTheory/Moduli.hs +10/−512
• Math/NumberTheory/Moduli/Chinese.hs +59/−0
• Math/NumberTheory/Moduli/Class.hs +306/−0
• Math/NumberTheory/Moduli/Jacobi.hs +150/−0
• Math/NumberTheory/Moduli/Sqrt.hs +229/−0
• Math/NumberTheory/MoebiusInversion.hs +6/−3
• Math/NumberTheory/MoebiusInversion/Int.hs +6/−3
• Math/NumberTheory/Powers.hs +0/−3
• Math/NumberTheory/Primes/Factorisation/Certified.hs +6/−3
• Math/NumberTheory/Primes/Factorisation/Montgomery.hs +112/−153
• Math/NumberTheory/Primes/Testing/Certificates/Internal.hs +12/−9
• Math/NumberTheory/Primes/Testing/Probabilistic.hs +20/−11
• Math/NumberTheory/Recurrencies/Linear.hs +34/−12
• arithmoi.cabal +17/−5
• benchmark/Math/NumberTheory/PrimesBench.hs +28/−5
• test-suite/Math/NumberTheory/ArithmeticFunctionsTests.hs +7/−1
• test-suite/Math/NumberTheory/CurvesTests.hs +99/−0
• test-suite/Math/NumberTheory/Moduli/ChineseTests.hs +48/−0
• test-suite/Math/NumberTheory/Moduli/ClassTests.hs +168/−0
• test-suite/Math/NumberTheory/Moduli/JacobiTests.hs +66/−0
• test-suite/Math/NumberTheory/Moduli/SqrtTests.hs +121/−0
• test-suite/Math/NumberTheory/ModuliTests.hs +0/−225
• test-suite/Math/NumberTheory/MoebiusInversion/IntTests.hs +8/−0
• test-suite/Math/NumberTheory/MoebiusInversionTests.hs +8/−0
• test-suite/Math/NumberTheory/Primes/FactorisationTests.hs +42/−7
• test-suite/Math/NumberTheory/TestUtils/Wrappers.hs +5/−3
• test-suite/Test.hs +12/−2

@@ -1,3 +1,38 @@+0.6.0.0:+    This release supports GHC 7.8, 7.10, 8.0 and 8.2.++    Breaking changes:++        'Math.NumberTheory.Moduli' was split into+        'Math.NumberTheory.Moduli.{Chinese,Class,Jacobi,Sqrt}'.++        Functions 'jacobi' and 'jacobi'' return 'JacobiSymbol'+        instead of 'Int'.++        Functions 'invertMod', 'powerMod' and 'powerModInteger' were removed,+        as well as their unchecked counterparts. Use new interface to+        modular computations, provided by 'Math.NumberTheory.Moduli.Class'.++    New functions:++        Brand new 'Math.NumberTheory.Moduli.Class' (#56), providing+        flexible and type safe modular arithmetic. Due to use of GMP built-ins+        it is also significantly faster.++        New function 'divisorsList', which is lazier than 'divisors' and+        does not require 'Ord' constraint (#64). Thus, it can be used+        for 'GaussianInteger'.++    Improvements:++        Speed up factorisation over elliptic curve up to 15x (#65).++        Polymorphic 'fibonacci' and 'lucas' functions, which previously+        were restricted to 'Integer' only (#63). This is especially useful+        for modular computations, e. g., 'map fibonacci [1..10] :: [Mod 7]'.++        Make 'totientSum' more robust and idiomatic (#58).+ 0.5.0.1:     Switch to QuickCheck 2.10. 

@@ -0,0 +1,36 @@+{-# LANGUAGE CPP         #-}++{-# OPTIONS_HADDOCK hide #-}++#if MIN_VERSION_base(4,10,0)++module GHC.TypeNats.Compat+  ( module GHC.TypeNats+  ) where++import GHC.TypeNats++#else++module GHC.TypeNats.Compat+  ( Nat+  , KnownNat+  , SomeNat(..)+  , natVal+  , someNatVal+  , sameNat+  ) where++import GHC.TypeLits (Nat, KnownNat, SomeNat(..), sameNat)+import qualified GHC.TypeLits as TL+import Numeric.Natural++natVal :: KnownNat n => proxy n -> Natural+natVal = fromInteger . TL.natVal++someNatVal :: Natural -> SomeNat+someNatVal n = case TL.someNatVal (toInteger n) of+  Nothing -> error "someNatVal: impossible negative argument"+  Just sn -> sn++#endif

@@ -20,6 +20,7 @@   ( -- * Multiplicative functions     multiplicative   , divisors, divisorsA+  , divisorsList, divisorsListA   , divisorsSmall, divisorsSmallA   , tau, tauA   , sigma, sigmaA@@ -76,6 +77,19 @@ divisorsHelper p a = S.fromDistinctAscList $ p : p * p : map (p ^) [3 .. wordToInt a] {-# INLINE divisorsHelper #-} +divisorsList :: (UniqueFactorisation n, Num n) => n -> [n]+divisorsList = runFunction divisorsListA++-- | The unsorted list of all (positive) divisors of an argument, produced in lazy fashion.+divisorsListA :: forall n. (UniqueFactorisation n, Num n) => ArithmeticFunction n [n]+divisorsListA = ArithmeticFunction (\((unPrime :: Prime n -> n) -> p) k -> ListProduct $ divisorsListHelper p k) ((1 :) . getListProduct)++divisorsListHelper :: Num n => n -> Word -> [n]+divisorsListHelper _ 0 = []+divisorsListHelper p 1 = [p]+divisorsListHelper p a = p : p * p : map (p ^) [3 .. wordToInt a]+{-# INLINE divisorsListHelper #-}+ divisorsSmall :: (UniqueFactorisation n, Prime n ~ Prime Int) => n -> IntSet divisorsSmall = runFunction divisorsSmallA @@ -279,6 +293,15 @@  instance (Num a, Ord a) => Monoid (SetProduct a) where   mempty  = SetProduct mempty+  mappend = (<>)++newtype ListProduct a = ListProduct { getListProduct :: [a] }++instance Num a => Semigroup (ListProduct a) where+  ListProduct s1 <> ListProduct s2 = ListProduct $ s1 <> s2 <> foldMap (\n -> map (* n) s2) s1++instance Num a => Monoid (ListProduct a) where+  mempty  = ListProduct mempty   mappend = (<>)  newtype IntSetProduct = IntSetProduct { getIntSetProduct :: IntSet }

@@ -0,0 +1,172 @@+-- |+-- Module:      Math.NumberTheory.Curves.Montgomery+-- Copyright:   (c) 2017 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Arithmetic on Montgomery elliptic curve.+--++{-# LANGUAGE BangPatterns        #-}+{-# LANGUAGE DataKinds           #-}+{-# LANGUAGE GADTs               #-}+{-# LANGUAGE KindSignatures      #-}+{-# LANGUAGE MagicHash           #-}+{-# LANGUAGE ScopedTypeVariables #-}++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.Curves.Montgomery+  ( Point+  , pointX+  , pointZ+  , pointN+  , pointA24+  , SomePoint(..)+  , newPoint+  , add+  , double+  , multiply+  ) where++import Data.Proxy+import GHC.Exts+import GHC.Integer.GMP.Internals+import GHC.Integer.Logarithms+import GHC.TypeNats.Compat++-- | We use the Montgomery form of elliptic curve:+-- b Y² = X³ + a X² + X (mod n).+-- See Eq. (10.3.1.1) at p. 260 of <http://www.ams.org/journals/mcom/1987-48-177/S0025-5718-1987-0866113-7/S0025-5718-1987-0866113-7.pdf Speeding the Pollard and Elliptic Curve Methods of Factorization> by P. L. Montgomery.+--+-- Switching to projective space by substitutions Y = y \/ z, X = x \/ z,+-- we get b y² z = x³ + a x² z + x z² (mod n).+-- The point on projective elliptic curve is characterized by three coordinates,+-- but it appears that only x- and z-components matter for computations.+-- By the same reason there is no need to store coefficient b.+--+-- That said, the chosen curve is represented by a24 = (a + 2) \/ 4+-- and modulo n at type level, making points on different curves+-- incompatible.+data Point (a24 :: Nat) (n :: Nat) = Point+  { pointX :: !Integer -- ^ Extract x-coordinate.+  , pointZ :: !Integer -- ^ Extract z-coordinate.+  }++pointA24 :: forall a24 n. KnownNat a24 => Point a24 n -> Integer+pointA24 _ = toInteger $ natVal (Proxy :: Proxy a24)++pointN :: forall a24 n. KnownNat n => Point a24 n -> Integer+pointN _ = toInteger $ natVal (Proxy :: Proxy n)++-- | In projective space 'Point's are equal, if they are both at infinity+-- or if respective ratios 'pointX' \/ 'pointZ' are equal.+instance KnownNat n => Eq (Point a24 n) where+  Point _ 0 == Point _ 0 = True+  Point _ 0 == _         = False+  _         == Point _ 0 = False+  p@(Point x1 z1) == Point x2 z2 = let n = pointN p in x1 * z2 `mod` n == x2 * z1 `mod` n++-- | For debugging.+instance (KnownNat a24, KnownNat n) => Show (Point a24 n) where+  show p = "(" ++ show (pointX p) ++ ", " ++ show (pointZ p) ++ ") (a24 "+    ++ show (pointA24 p) ++ ", mod "+    ++ show (pointN p) ++ ")"++-- | Point on unknown curve.+data SomePoint where+  SomePoint :: (KnownNat a24, KnownNat n) => Point a24 n -> SomePoint++instance Show SomePoint where+  show (SomePoint p) = show p++-- | 'newPoint' @s@ @n@ creates a point on an elliptic curve modulo @n@, uniquely determined by seed @s@.+-- Some choices of @s@ and @n@ produce ill-parametrized curves, which is reflected by return value 'Nothing'.+--+-- We choose a curve by Suyama's parametrization. See Eq. (3)-(4) at p. 4+-- of <http://www.hyperelliptic.org/tanja/SHARCS/talks06/Gaj.pdf Implementing the Elliptic Curve Method of Factoring in Reconfigurable Hardware>+-- by K. Gaj, S. Kwon et al.+newPoint :: Integer -> Integer -> Maybe SomePoint+newPoint s n = do+    a24denRecip <- case recipModInteger a24den n of+      0 -> Nothing+      t -> Just t+    a24 <- case a24num * a24denRecip `rem` n of+      -- (a+2)/4 = 0 corresponds to singular curve with A = -2+      0 -> Nothing+      -- (a+2)/4 = 1 corresponds to singular curve with A = 2+      1 -> Nothing+      t -> Just t+    SomeNat (_ :: Proxy a24Ty) <- if a24 < 0+                                  then Nothing+                                  else Just $ someNatVal $ fromInteger a24+    SomeNat (_ :: Proxy nTy)   <- if n < 0+                                  then Nothing+                                  else Just $ someNatVal $ fromInteger n+    return $ SomePoint (Point x z :: Point a24Ty nTy)+  where+    u = s * s `rem` n - 5+    v = 4 * s+    d = v - u+    x = u * u * u `mod` n+    z = v * v * v `mod` n+    a24num = d * d * d * (3 * u + v) `mod` n+    a24den = 16 * x * v `rem` n++-- | If @p0@ + @p1@ = @p2@, then 'add' @p0@ @p1@ @p2@ equals to @p1@ + @p2@.+-- It is also required that z-coordinates of @p0@, @p1@ and @p2@ are coprime with modulo+-- of elliptic curve; and x-coordinate of @p0@ is non-zero.+-- If preconditions do not hold, return value is undefined.+--+-- Remarkably such addition does not require 'KnownNat' @a24@ constraint.+--+-- Computations follow Algorithm 3 at p. 4+-- of <http://www.hyperelliptic.org/tanja/SHARCS/talks06/Gaj.pdf Implementing the Elliptic Curve Method of Factoring in Reconfigurable Hardware>+-- by K. Gaj, S. Kwon et al.+add :: KnownNat n => Point a24 n -> Point a24 n -> Point a24 n -> Point a24 n+add p0@(Point x0 z0) (Point x1 z1) (Point x2 z2) = Point x3 z3+  where+    n = pointN p0+    a = (x1 - z1) * (x2 + z2) `rem` n+    b = (x1 + z1) * (x2 - z2) `rem` n+    apb = a + b+    amb = a - b+    c = apb * apb `rem` n+    d = amb * amb `rem` n+    x3 = c * z0 `mod` n+    z3 = d * x0 `mod` n++-- | Multiply by 2.+--+-- Computations follow Algorithm 3 at p. 4+-- of <http://www.hyperelliptic.org/tanja/SHARCS/talks06/Gaj.pdf Implementing the Elliptic Curve Method of Factoring in Reconfigurable Hardware>+-- by K. Gaj, S. Kwon et al.+double :: (KnownNat a24, KnownNat n) => Point a24 n -> Point a24 n+double p@(Point x z) = Point x' z'+  where+    n = pointN p+    a24 = pointA24 p+    r = x + z+    s = x - z+    u = r * r `rem` n+    v = s * s `rem` n+    t = u - v+    x' = u * v `mod` n+    z' = (v + a24 * t `rem` n) * t `mod` n++-- | Multiply by given number, using binary algorithm.+multiply :: (KnownNat a24, KnownNat n) => Word -> Point a24 n -> Point a24 n+multiply 0 _ = Point 0 0+multiply 1 p = p+multiply (W# w##) p =+    case wordLog2# w## of+      l# -> go (l# -# 1#) p (double p)+  where+    go 0# !p0 !p1 = case w## `and#` 1## of+                      0## -> double p0+                      _   -> add p p0 p1+    go i# p0 p1 = case uncheckedShiftRL# w## i# `and#` 1## of+                    0## -> go (i# -# 1#) (double p0) (add p p0 p1)+                    _   -> go (i# -# 1#) (add p p0 p1) (double p1)

@@ -8,517 +8,15 @@ -- -- Miscellaneous functions related to modular arithmetic. ---{-# LANGUAGE CPP, BangPatterns #-}-module Math.NumberTheory.Moduli-    ( -- * Functions with input check-      jacobi-    , invertMod-    , powerMod-    , powerModInteger-    , chineseRemainder-      -- ** Partially checked input-    , sqrtModP-      -- * Unchecked functions-    , jacobi'-    , powerMod'-    , powerModInteger'-    , sqrtModPList-    , sqrtModP'-    , tonelliShanks-    , sqrtModPP-    , sqrtModPPList-    , sqrtModF-    , sqrtModFList-    , chineseRemainder2-    ) where -#include "MachDeps.h"--#if __GLASGOW_HASKELL__ < 709 || WORD_SIZE_IN_BITS == 32-import Data.Word-#endif-import Data.Bits-import Data.Array.Unboxed-import Data.List (nub)-import Control.Monad (foldM, liftM2)--import Math.NumberTheory.Utils (shiftToOddCount, splitOff)-import Math.NumberTheory.GCD (extendedGCD)-import Math.NumberTheory.Primes.Heap (sieveFrom)-import Math.NumberTheory.Unsafe---- Guesstimated startup time for the Heap algorithm is lower than--- the cost to sieve an entire chunk.---- | Invert a number relative to a positive modulus.---   If @number@ and @modulus@ are coprime, the result is---   @Just inverse@ where------ >    (number * inverse) `mod` modulus == 1--- >    0 <= inverse < modulus------   If @number `mod` modulus == 0@ or @gcd number modulus > 1@, the result is @Nothing@.-invertMod :: Integer -> Integer -> Maybe Integer-invertMod k m-  | m <= 0 = error "Math.NumberTheory.Moduli.invertMod: non-positive modulus"-  | otherwise = wrap $ go False 1 0 m k'-  where-    k' | r < 0     = r+m-       | otherwise = r-         where-           r = k `rem` m-    wrap x = case (x*k') `rem` m of-               1 -> Just x-               _ -> Nothing-    -- Calculate modular inverse of k' modulo m by continued fraction expansion-    -- of m/k', say [a_0,a_1,...,a_s]. Let the convergents be p_j/q_j.-    -- Starting from j = -2, the arguments of go are-    -- (p_j/q_j) > m/k', p_{j+1}, p_j, and n, d with n/d = [a_{j+2},...,a_s].-    -- Since m/k' = p_s/q_s, and p_j*q_{j+1} - p_{j+1}*q_j = (-1)^(j+1), we have-    -- p_{s-1}*k' - q_{s-1}*m = (-1)^s * gcd m k', so if the inverse exists,-    -- it is either p_{s-1} or -p_{s-1}, depending on whether s is even or odd.-    go !b _ po _ 0 = if b then po else (m-po)-    go b !pn po n d = case n `quotRem` d of-                        (q,r) -> go (not b) (q*pn+po) pn d r---- | Jacobi symbol of two numbers.---   The \"denominator\" must be odd and positive, this condition is checked.------   If both numbers have a common prime factor, the result---   is @0@, otherwise it is &#177;1.-{-# SPECIALISE jacobi :: Integer -> Integer -> Int,-                         Int -> Int -> Int,-                         Word -> Word -> Int-  #-}-jacobi :: (Integral a, Bits a) => a -> a -> Int-jacobi a b-  | b < 0       = error "Math.NumberTheory.Moduli.jacobi: negative denominator"-  | evenI b     = error "Math.NumberTheory.Moduli.jacobi: even denominator"-  | b == 1      = 1-  | a == 0      = 0-  | a == 1      = 1-  | otherwise   = jacobi' a b   -- b odd, > 1, a neither 0 or 1---- Invariant: b > 1 and odd--- | Jacobi symbol of two numbers without validity check of---   the \"denominator\".-{-# SPECIALISE jacobi' :: Integer -> Integer -> Int,-                          Int -> Int -> Int,-                          Word -> Word -> Int-  #-}-jacobi' :: (Integral a, Bits a) => a -> a -> Int-jacobi' a b-  | a == 0      = 0-  | a == 1      = 1-  | a < 0       = let n | rem4 b == 1 = 1-                        | otherwise   = -1-                      -- Blech, minBound may pose problems-                      (z,o) = shiftToOddCount (abs $ toInteger a)-                      s | evenI z || unsafeAt jac2 (rem8 b) == 1 = n-                        | otherwise                              = (-n)-                  in s*jacobi' (fromInteger o) b-  | a >= b      = case a `rem` b of-                    0 -> 0-                    r -> jacPS 1 r b-  | evenI a     = case shiftToOddCount a of-                    (z,o) -> let r = 2 - (rem4 o .&. rem4 b)-                                 s | evenI z || unsafeAt jac2 (rem8 b) == 1 = r-                                   | otherwise                              = (-r)-                             in jacOL s b o-  | otherwise   = case rem4 a .&. rem4 b of-                    3 -> jacOL (-1) b a-                    _ -> jacOL 1 b a---- numerator positive and smaller than denominator-{-# SPECIALISE jacPS :: Int -> Integer -> Integer -> Int,-                        Int -> Int -> Int -> Int,-                        Int -> Word -> Word -> Int-  #-}-jacPS :: (Integral a, Bits a) => Int -> a -> a -> Int-jacPS !j a b-  | evenI a     = case shiftToOddCount a of-                    (z,o) | evenI z || unsafeAt jac2 (rem8 b) == 1 ->-                              jacOL (if rem4 o .&. rem4 b == 3 then (-j) else j) b o-                          | otherwise ->-                              jacOL (if rem4 o .&. rem4 b == 3 then j else (-j)) b o-  | otherwise   = jacOL (if rem4 a .&. rem4 b == 3 then (-j) else j) b a---- numerator odd, positive and larger than denominator-{-# SPECIALISE jacOL :: Int -> Integer -> Integer -> Int,-                        Int -> Int -> Int -> Int,-                        Int -> Word -> Word -> Int-  #-}-jacOL :: (Integral a, Bits a) => Int -> a -> a -> Int-jacOL !j a b-  | b == 1    = j-  | otherwise = case a `rem` b of-                 0 -> 0-                 r -> jacPS j r b---- | Modular power.------ > powerMod base exponent modulus------   calculates @(base ^ exponent) \`mod\` modulus@ by repeated squaring and reduction. Modulus must be positive.---   If @exponent < 0@ and @base@ is invertible modulo @modulus@, @(inverse ^ |exponent|) \`mod\` modulus@---   is calculated. This function does some input checking and sanitation before calling the unsafe worker.-{-# RULES-"powerMod/Integer" powerMod = powerModInteger-  #-}-{-# INLINE [1] powerMod #-}-powerMod :: (Integral a, Bits a) => Integer -> a -> Integer -> Integer-powerMod = powerModImpl--{-# SPECIALISE powerModImpl :: Integer -> Int -> Integer -> Integer,-                               Integer -> Word -> Integer -> Integer-  #-}-powerModImpl :: (Integral a, Bits a) => Integer -> a -> Integer -> Integer-powerModImpl base expo md-  | md <= 0     = error "Math.NumberTheory.Moduli.powerMod: non-positive modulus"-  | md == 1     = 0-  | expo == 0   = 1-  | bse' == 1   = 1-  | expo < 0    = case invertMod bse' md of-                    Just i  -> powerMod'Impl i (negate expo) md-                    Nothing -> error "Math.NumberTheory.Moduli.powerMod: Base isn't invertible with respect to modulus"-  | bse' == 0   = 0-  | otherwise   = powerMod'Impl bse' expo md-    where-      bse' = if base < 0 || md <= base then base `mod` md else base---- | Modular power worker without input checking.---   Assumes all arguments strictly positive and modulus greater than 1.-{-# RULES-"powerMod'/Integer" powerMod' = powerModInteger'-  #-}-{-# INLINE [1] powerMod' #-}-powerMod' :: (Integral a, Bits a) => Integer -> a -> Integer -> Integer-powerMod' = powerMod'Impl---{-# SPECIALISE powerMod'Impl :: Integer -> Int -> Integer -> Integer,-                                Integer -> Word -> Integer -> Integer-  #-}-powerMod'Impl :: (Integral a, Bits a) => Integer -> a -> Integer -> Integer-powerMod'Impl base expo md = go expo 1 base-  where-    go 1 !a !s  = (a*s) `rem` md-    go e a s-      | testBit e 0 = go (e `shiftR` 1) ((a*s) `rem` md) ((s*s) `rem` md)-      | otherwise   = go (e `shiftR` 1) a ((s*s) `rem` md)---- | Specialised version of 'powerMod' for 'Integer' exponents.---   Reduces the number of shifts of the exponent since shifting---   large 'Integer's is expensive. Call this function directly---   if you don't want or can't rely on rewrite rules. Modulus must be positive.-powerModInteger :: Integer -> Integer -> Integer -> Integer-powerModInteger base ex mdl-  | mdl <= 0     = error "Math.NumberTheory.Moduli.powerModInteger: non-positive modulus"-  | mdl == 1    = 0-  | ex == 0     = 1-  | ex < 0      = case invertMod bse' mdl of-                    Just i  -> powerModInteger' i (negate ex) mdl-                    Nothing -> error "Math.NumberTheory.Moduli.powerMod: Base isn't invertible with respect to modulus"-  | bse' == 0   = 0-  | bse' == 1   = 1-  | otherwise   = powerModInteger' bse' ex mdl-    where-      bse' = if base < 0 || mdl <= base then base `mod` mdl else base---- | Specialised worker without input checks. Makes the same assumptions---   as the general version 'powerMod''.-powerModInteger' :: Integer -> Integer -> Integer -> Integer-powerModInteger' base expo md = go w1 1 base e1-  where-    w1 = fromInteger expo-    e1 = expo `shiftR` 64-#if WORD_SIZE_IN_BITS == 32-  -- Shifting large Integers is expensive, hence we reduce the-  -- number of shifts by processing in 64-bit chunks.-  -- On 32-bit systems, every testBit on a Word64 would be a C-call,-  -- thus it is faster to split each Word64 into the constituent 32-bit-  -- Words and process those separately.-  -- The code becomes ugly, unfortunately.-    go :: Word64 -> Integer -> Integer -> Integer -> Integer-    go !w !a !s 0  = end a s w-    go w a s e = inner1 a s 0-      where-        wl :: Word-        !wl = fromIntegral w-        wh :: Word-        !wh = fromIntegral (w `shiftR` 32)-        inner1 !au !sq 32 = inner2 au sq 0-        inner1 au sq i-          | testBit wl i = inner1 ((au*sq) `rem` md) ((sq*sq) `rem` md) (i+1)-          | otherwise    = inner1 au ((sq*sq) `rem` md) (i+1)-        inner2 !au !sq 32 = go (fromInteger e) au sq (e `shiftR` 64)-        inner2 au sq i-          | testBit wh i = inner2 ((au*sq) `rem` md) ((sq*sq) `rem` md) (i+1)-          | otherwise    = inner2 au ((sq*sq) `rem` md) (i+1)-    end !a !s w-      | wh == 0   = fin a s wl-      | otherwise = innerE a s 0-        where-          wl :: Word-          !wl = fromIntegral w-          wh :: Word-          !wh = fromIntegral (w `shiftR` 32)-          innerE !au !sq 32 = fin au sq wh-          innerE au sq i-            | testBit wl i = innerE ((au*sq) `rem` md) ((sq*sq) `rem` md) (i+1)-            | otherwise    = innerE au ((sq*sq) `rem` md) (i+1)-    fin :: Integer -> Integer -> Word -> Integer-    fin !a !s 1 = (a*s) `rem` md-    fin a s w-      | testBit w 0 = fin ((a*s) `rem` md) ((s*s) `rem` md) (w `shiftR` 1)-      | otherwise   = fin a ((s*s) `rem` md) (w `shiftR` 1)--#else-  -- WORD_SIZE_IN_BITS == 64, otherwise things wouldn't compile anyway-  -- Shorter code since we need not split each 64-bit word.-    go :: Word -> Integer -> Integer -> Integer -> Integer-    go !w !a !s 0  = end a s w-    go w a s e = inner a s 0-      where-        inner !au !sq 64 = go (fromInteger e) au sq (e `shiftR` 64)-        inner au sq i-          | testBit w i = inner ((au*sq) `rem` md) ((sq*sq) `rem` md) (i+1)-          | otherwise   = inner au ((sq*sq) `rem` md) (i+1)-    end !a !s 1 = (a*s) `rem` md-    end a s w-      | testBit w 0 = end ((a*s) `rem` md) ((s*s) `rem` md) (w `shiftR` 1)-      | otherwise   = end a ((s*s) `rem` md) (w `shiftR` 1)--#endif---- | @sqrtModP n prime@ calculates a modular square root of @n@ modulo @prime@---   if that exists. The second argument /must/ be a (positive) prime, otherwise---   the computation may not terminate and if it does, may yield a wrong result.---   The precondition is /not/ checked.------   If @prime@ is a prime and @n@ a quadratic residue modulo @prime@, the result---   is @Just r@ where @r^2 ≡ n (mod prime)@, if @n@ is a quadratic nonresidue,---   the result is @Nothing@.-sqrtModP :: Integer -> Integer -> Maybe Integer-sqrtModP n 2 = Just (n `mod` 2)-sqrtModP n prime = case jacobi' n prime of-                     0 -> Just 0-                     1 -> Just (sqrtModP' (n `mod` prime) prime)-                     _ -> Nothing---- | @sqrtModPList n prime@ computes the list of all square roots of @n@---   modulo @prime@. @prime@ /must/ be a (positive) prime.---   The precondition is /not/ checked.-sqrtModPList :: Integer -> Integer -> [Integer]-sqrtModPList n prime-    | prime == 2    = [n `mod` 2]-    | otherwise     = case sqrtModP n prime of-                        Just 0 -> [0]-                        Just r -> [r,prime-r] -- The group of units in Z/(p) is cyclic-                        _      -> []---- | @sqrtModP' square prime@ finds a square root of @square@ modulo---   prime. @prime@ /must/ be a (positive) prime, and @square@ /must/ be a positive---   quadratic residue modulo @prime@, i.e. @'jacobi square prime == 1@.---   The precondition is /not/ checked.-sqrtModP' :: Integer -> Integer -> Integer-sqrtModP' square prime-    | prime == 2    = square-    | rem4 prime == 3 = powerModInteger' square ((prime + 1) `quot` 4) prime-    | otherwise     = tonelliShanks square prime---- | @tonelliShanks square prime@ calculates a square root of @square@---   modulo @prime@, where @prime@ is a prime of the form @4*k + 1@ and---   @square@ is a positive quadratic residue modulo @prime@, using the---   Tonelli-Shanks algorithm.---   No checks on the input are performed.-tonelliShanks :: Integer -> Integer -> Integer-tonelliShanks square prime = loop rc t1 generator log2-  where-    (log2,q) = shiftToOddCount (prime-1)-    nonSquare = findNonSquare prime-    generator = powerModInteger' nonSquare q prime-    rc = powerModInteger' square ((q+1) `quot` 2) prime-    t1 = powerModInteger' square q prime-    msqr x = (x*x) `rem` prime-    msquare 0 x = x-    msquare k x = msquare (k-1) (msqr x)-    findPeriod per 1 = per-    findPeriod per x = findPeriod (per+1) (msqr x)-    loop !r t c m-        | t == 1    = r-        | otherwise = loop nextR nextT nextC nextM-          where-            nextM = findPeriod 0 t-            b     = msquare (m - 1 - nextM) c-            nextR = (r*b) `rem` prime-            nextC = msqr b-            nextT = (t*nextC) `rem` prime---- | @sqrtModPP n (prime,expo)@ calculates a square root of @n@---   modulo @prime^expo@ if one exists. @prime@ /must/ be a---   (positive) prime. @expo@ must be positive, @n@ must be coprime---   to @prime@-sqrtModPP :: Integer -> (Integer,Int) -> Maybe Integer-sqrtModPP n (2,e) = sqM2P n e-sqrtModPP n (prime,expo) = case sqrtModP n prime of-                             Just r -> fixup r-                             _      -> Nothing-  where-    fixup r = let diff' = r*r-n-              in if diff' == 0-                   then Just r-                   else case splitOff prime diff' of-                          (e,q) | expo <= e -> Just r-                                | otherwise -> fmap (\inv -> hoist inv r (q `mod` prime) (prime^e)) (invertMod (2*r) prime)-                      ---    hoist inv root elim pp-        | diff' == 0    = root'-        | expo <= ex    = root'-        | otherwise     = hoist inv root' (nelim `mod` prime) (prime^ex)-          where-            root' = (root + (inv*(prime-elim))*pp) `mod` (prime*pp)-            diff' = root'*root' - n-            (ex, nelim) = splitOff prime diff'---- dirty, dirty-sqM2P :: Integer -> Int -> Maybe Integer-sqM2P n e-    | e < 2     = Just (n `mod` 2)-    | n' == 0   = Just 0-    | e <= k    = Just 0-    | odd k     = Nothing-    | otherwise = fmap ((`mod` mdl) . (`shiftL` k2)) $ solve s e2-      where-        mdl = 1 `shiftL` e-        n' = n `mod` mdl-        (k,s) = shiftToOddCount n'-        k2 = k `quot` 2-        e2 = e-k-        solve _ 1 = Just 1-        solve 1 _ = Just 1-        solve r p-            | rem4 r == 3   = Nothing  -- otherwise r ≡ 1 (mod 4)-            | p == 2        = Just 1   -- otherwise p >= 3-            | rem8 r == 5   = Nothing  -- otherwise r ≡ 1 (mod 8)-            | otherwise     = fixup r (fst $ shiftToOddCount (r-1))-              where-                fixup x pw-                    | pw >= e2  = Just x-                    | otherwise = fixup x' pw'-                      where-                        x' = x + (1 `shiftL` (pw-1))-                        d = x'*x' - r-                        pw' = if d == 0 then e2 else fst (shiftToOddCount d)---- | @sqrtModF n primePowers@ calculates a square root of @n@ modulo---   @product [p^k | (p,k) <- primePowers]@ if one exists and all primes---   are distinct.---   The list must be non-empty, @n@ must be coprime with all primes.-sqrtModF :: Integer -> [(Integer,Int)] -> Maybe Integer-sqrtModF _ []  = Nothing-sqrtModF n pps = do roots <- mapM (sqrtModPP n) pps-                    chineseRemainder $ zip roots (map (uncurry (^)) pps)---- | @sqrtModFList n primePowers@ calculates all square roots of @n@ modulo---   @product [p^k | (p,k) <- primePowers]@ if all primes are distinct.---   The list must be non-empty, @n@ must be coprime with all primes.-sqrtModFList :: Integer -> [(Integer,Int)] -> [Integer]-sqrtModFList _ []  = []-sqrtModFList n pps = map fst $ foldl1 (liftM2 comb) cs-  where-    ms :: [Integer]-    ms = map (uncurry (^)) pps-    rs :: [[Integer]]-    rs = map (sqrtModPPList n) pps-    cs :: [[(Integer,Integer)]]-    cs = zipWith (\l m -> map (\x -> (x,m)) l) rs ms-    comb t1@(_,m1) t2@(_,m2) = (chineseRemainder2 t1 t2,m1*m2)---- | @sqrtModPPList n (prime,expo)@ calculates the list of all---   square roots of @n@ modulo @prime^expo@. The same restriction---   as in 'sqrtModPP' applies to the arguments.-sqrtModPPList :: Integer -> (Integer,Int) -> [Integer]-sqrtModPPList n (2,1) = [n `mod` 2]-sqrtModPPList n (2,expo)-    = case sqM2P n expo of-        Just r -> let m = 1 `shiftL` (expo-1)-                  in nub [r, (r+m) `mod` (2*m), (m-r) `mod` (2*m), 2*m-r]-        _ -> []-sqrtModPPList n pe@(prime,expo)-    = case sqrtModPP n pe of-        Just 0 -> [0]-        Just r -> [prime^expo - r, r] -- The group of units in Z/(p^e) is cyclic-        _      -> []---- | Given a list @[(r_1,m_1), ..., (r_n,m_n)]@ of @(residue,modulus)@---   pairs, @chineseRemainder@ calculates the solution to the simultaneous---   congruences------ >--- > r ≡ r_k (mod m_k)--- >------   if all moduli are positive and pairwise coprime. Otherwise---   the result is @Nothing@ regardless of whether---   a solution exists.-chineseRemainder :: [(Integer,Integer)] -> Maybe Integer-chineseRemainder remainders = foldM addRem 0 remainders-  where-    !modulus = product (map snd remainders)-    addRem acc (_,1) = Just acc-    addRem acc (r,m) = do-        let cf = modulus `quot` m-        inv <- invertMod cf m-        Just $! (acc + inv*cf*r) `mod` modulus---- | @chineseRemainder2 (r_1,m_1) (r_2,m_2)@ calculates the solution of------ >--- > r ≡ r_k (mod m_k)------   if @m_1@ and @m_2@ are coprime.-chineseRemainder2 :: (Integer,Integer) -> (Integer,Integer) -> Integer-chineseRemainder2 (r1, md1) (r2,md2)-    = case extendedGCD md1 md2 of-        (_,u,v) -> ((1 - u*md1)*r1 + (1 - v*md2)*r2) `mod` (md1*md2)---- Utilities---- For large Integers, going via Int is much faster than bit-fiddling--- on the Integer, so we do that.-{-# SPECIALISE evenI :: Integer -> Bool,-                        Int -> Bool,-                        Word -> Bool-  #-}-evenI :: Integral a => a -> Bool-evenI n = fromIntegral n .&. 1 == (0 :: Int)--{-# SPECIALISE rem4 :: Integer -> Int,-                       Int -> Int,-                       Word -> Int-  #-}-rem4 :: Integral a => a -> Int-rem4 n = fromIntegral n .&. 3--{-# SPECIALISE rem8 :: Integer -> Int,-                       Int -> Int,-                       Word -> Int-  #-}-rem8 :: Integral a => a -> Int-rem8 n = fromIntegral n .&. 7--jac2 :: UArray Int Int-jac2 = array (0,7) [(0,0),(1,1),(2,0),(3,-1),(4,0),(5,-1),(6,0),(7,1)]+module Math.NumberTheory.Moduli+  ( module Math.NumberTheory.Moduli.Class+  , module Math.NumberTheory.Moduli.Chinese+  , module Math.NumberTheory.Moduli.Jacobi+  , module Math.NumberTheory.Moduli.Sqrt+  ) where -findNonSquare :: Integer -> Integer-findNonSquare n-    | rem8 n == 5 || rem8 n == 3  = 2-    | otherwise = search primelist-      where-        primelist = [3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67]-                        ++ sieveFrom (68 + n `rem` 4) -- prevent sharing-        search (p:ps)-            | jacobi' p n == -1 = p-            | otherwise         = search ps-        search _ = error "Should never have happened, prime list exhausted."+import Math.NumberTheory.Moduli.Chinese+import Math.NumberTheory.Moduli.Class+import Math.NumberTheory.Moduli.Jacobi+import Math.NumberTheory.Moduli.Sqrt

@@ -0,0 +1,59 @@+-- |+-- Module:      Math.NumberTheory.Moduli.Chinese+-- Copyright:   (c) 2011 Daniel Fischer+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Chinese remainder theorem+--++{-# LANGUAGE BangPatterns #-}++module Math.NumberTheory.Moduli.Chinese+  ( chineseRemainder+  , chineseRemainder2+  ) where++import Control.Monad (foldM)+import GHC.Integer.GMP.Internals++import Math.NumberTheory.GCD (extendedGCD)++-- | Given a list @[(r_1,m_1), ..., (r_n,m_n)]@ of @(residue,modulus)@+--   pairs, @chineseRemainder@ calculates the solution to the simultaneous+--   congruences+--+-- >+-- > r ≡ r_k (mod m_k)+-- >+--+--   if all moduli are positive and pairwise coprime. Otherwise+--   the result is @Nothing@ regardless of whether+--   a solution exists.+chineseRemainder :: [(Integer,Integer)] -> Maybe Integer+chineseRemainder remainders = foldM addRem 0 remainders+  where+    !modulus = product (map snd remainders)+    addRem acc (_,1) = Just acc+    addRem acc (r,m) = do+        let cf = modulus `quot` m+        inv <- recipMod cf m+        Just $! (acc + inv*cf*r) `mod` modulus++-- | @chineseRemainder2 (r_1,m_1) (r_2,m_2)@ calculates the solution of+--+-- >+-- > r ≡ r_k (mod m_k)+--+--   if @m_1@ and @m_2@ are coprime.+chineseRemainder2 :: (Integer,Integer) -> (Integer,Integer) -> Integer+chineseRemainder2 (r1, md1) (r2,md2)+    = case extendedGCD md1 md2 of+        (_,u,v) -> ((1 - u*md1)*r1 + (1 - v*md2)*r2) `mod` (md1*md2)++recipMod :: Integer -> Integer -> Maybe Integer+recipMod x m = case recipModInteger x m of+  0 -> Nothing+  y -> Just y

@@ -0,0 +1,306 @@+-- |+-- Module:      Math.NumberTheory.Moduli.Class+-- Copyright:   (c) 2017 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Safe modular arithmetic with modulo on type level.+--++{-# LANGUAGE CPP                 #-}+{-# LANGUAGE DataKinds           #-}+{-# LANGUAGE GADTs               #-}+{-# LANGUAGE KindSignatures      #-}+{-# LANGUAGE LambdaCase          #-}+{-# LANGUAGE RankNTypes          #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE StandaloneDeriving  #-}++module Math.NumberTheory.Moduli.Class+  ( -- * Known modulo+    Mod+  , getVal+  , getNatVal+  , getMod+  , getNatMod+  , invertMod+  , powMod+  , (^%)+  -- * Unknown modulo+  , SomeMod(..)+  , modulo+  , invertSomeMod+  , powSomeMod+  -- * Re-exported from GHC.TypeNats.Compat+  , KnownNat+  ) where++import Data.Proxy+import Data.Ratio+import Data.Type.Equality+#if __GLASGOW_HASKELL__ < 709+import Data.Word+#endif+import GHC.Integer.GMP.Internals+import GHC.TypeNats.Compat+import Numeric.Natural++-- | Wrapper for residues modulo @m@.+--+-- @Mod 3 :: Mod 10@ stands for the class of integers, congruent to 3 modulo 10 (…−17, −7, 3, 13, 23…).+-- The modulo is stored on type level, so it is impossible, for example, to add up by mistake+-- residues with different moduli.+--+-- > > (3 :: Mod 10) + (4 :: Mod 12)+-- > error: Couldn't match type ‘12’ with ‘10’...+-- > > (3 :: Mod 10) + 8+-- > (1 `modulo` 10)+--+-- Note that modulo cannot be negative.+newtype Mod (m :: Nat) = Mod Natural+  deriving (Eq, Ord)++instance KnownNat m => Show (Mod m) where+  show m = "(" ++ show (getVal m) ++ " `modulo` " ++ show (getMod m) ++ ")"++instance KnownNat m => Num (Mod m) where+  mx@(Mod x) + Mod y =+    Mod $ if xy >= m then xy - m else xy+    where+      xy = x + y+      m = getNatMod mx+  {-# INLINE (+) #-}+  mx@(Mod x) - Mod y =+    Mod $ if x >= y then x - y else m + x - y+    where+      m = getNatMod mx+  {-# INLINE (-) #-}+  negate mx@(Mod x) =+    Mod $ if x == 0 then 0 else getNatMod mx - x+  {-# INLINE negate #-}+  mx@(Mod x) * Mod y =+    Mod $ x * y `rem` getNatMod mx -- `rem` is slightly faster than `mod`+  {-# INLINE (*) #-}+  abs = id+  {-# INLINE abs #-}+  signum = const $ Mod 1+  {-# INLINE signum #-}+  fromInteger x = mx+    where+      mx = Mod $ fromInteger $ x `mod` getMod mx+  {-# INLINE fromInteger #-}++-- | Beware that division by residue, which is not coprime with the modulo,+-- will result in runtime error. Consider using 'invertMod' instead.+instance KnownNat m => Fractional (Mod m) where+  fromRational r = case denominator r of+    1   -> num+    den -> num / fromInteger den+    where+      num = fromInteger (numerator r)+  {-# INLINE fromRational #-}+  recip mx = case invertMod mx of+    Nothing -> error $ "recip{Mod}: residue is not coprime with modulo"+    Just y  -> y+  {-# INLINE recip #-}++-- | Linking type and value levels: extract modulo @m@ as a value.+getMod :: KnownNat m => Mod m -> Integer+getMod = toInteger . natVal+{-# INLINE getMod #-}++-- | Linking type and value levels: extract modulo @m@ as a value.+getNatMod :: KnownNat m => Mod m -> Natural+getNatMod = natVal+{-# INLINE getNatMod #-}++-- | The canonical representative of the residue class, always between 0 and @m-1@ inclusively.+getVal :: KnownNat m => Mod m -> Integer+getVal (Mod x) = toInteger x+{-# INLINE getVal #-}++-- | The canonical representative of the residue class, always between 0 and @m-1@ inclusively.+getNatVal :: KnownNat m => Mod m -> Natural+getNatVal (Mod x) = x+{-# INLINE getNatVal #-}++-- | Computes the modular inverse, if the residue is coprime with the modulo.+--+-- > > invertMod (3 :: Mod 10)+-- > Just (7 `modulo` 10) -- because 3 * 7 = 1 :: Mod 10+-- > > invertMod (4 :: Mod 10)+-- > Nothing+invertMod :: KnownNat m => Mod m -> Maybe (Mod m)+invertMod mx+  = if y <= 0+    then Nothing+    else Just $ Mod $ fromInteger y+  where+    y = recipModInteger (getVal mx) (getMod mx)+{-# INLINABLE invertMod #-}++-- | Drop-in replacement for '^', with much better performance.+--+-- > > powMod (3 :: Mod 10) 4+-- > (1 `modulo` 10)+powMod :: (KnownNat m, Integral a) => Mod m -> a -> Mod m+powMod mx a+  | a < 0     = error $ "^{Mod}: negative exponent"+  | otherwise = Mod $ fromInteger $ powModInteger (getVal mx) (toInteger a) (getMod mx)+{-# INLINABLE [1] powMod #-}++{-# SPECIALISE [1] powMod ::+  KnownNat m => Mod m -> Integer -> Mod m,+  KnownNat m => Mod m -> Natural -> Mod m,+  KnownNat m => Mod m -> Int     -> Mod m,+  KnownNat m => Mod m -> Word    -> Mod m #-}++{-# RULES+"powMod/2/Integer"     forall x. powMod x (2 :: Integer) = let u = x in u*u+"powMod/3/Integer"     forall x. powMod x (3 :: Integer) = let u = x in u*u*u+"powMod/2/Int"         forall x. powMod x (2 :: Int)     = let u = x in u*u+"powMod/3/Int"         forall x. powMod x (3 :: Int)     = let u = x in u*u*u #-}++-- | Infix synonym of 'powMod'.+(^%) :: (KnownNat m, Integral a) => Mod m -> a -> Mod m+(^%) = powMod+{-# INLINE (^%) #-}++infixr 8 ^%++-- Unfortunately, such rule never fires due to technical details+-- of type classes in Core.+-- {-# RULES "^%Mod" forall (x :: KnownNat m => Mod m) p. x ^ p = x ^% p #-}++-- | This type represents residues with unknown modulo and rational numbers.+-- One can freely combine them in arithmetic expressions, but each operation+-- will spend time on modulo's recalculation:+--+-- > > 2 `modulo` 10 + 4 `modulo` 15+-- > (1 `modulo` 5)+-- > > 2 `modulo` 10 * 4 `modulo` 15+-- > (3 `modulo` 5)+-- > > 2 `modulo` 10 + fromRational (3 % 7)+-- > (1 `modulo` 10)+-- > > 2 `modulo` 10 * fromRational (3 % 7)+-- > (8 `modulo` 10)+--+-- If performance is crucial, it is recommended to extract @Mod m@ for further processing+-- by pattern matching. E. g.,+--+-- > case modulo n m of+-- >   SomeMod k -> process k -- Here k has type Mod m+-- >   InfMod{}  -> error "impossible"+data SomeMod where+  SomeMod :: KnownNat m => Mod m -> SomeMod+  InfMod  :: Rational -> SomeMod++instance Eq SomeMod where+  SomeMod mx == SomeMod my = getMod mx == getMod my && getVal mx == getVal my+  InfMod rx  == InfMod ry  = rx == ry+  _          == _          = False++instance Show SomeMod where+  show = \case+    SomeMod m -> show m+    InfMod  r -> show r++-- | Create modular value by representative of residue class and modulo.+-- One can use the result either directly (via functions from 'Num' and 'Fractional'),+-- or deconstruct it by pattern matching. Note that 'modulo' never returns 'InfMod'.+modulo :: Integer -> Natural -> SomeMod+modulo n m = case someNatVal m of+  SomeNat (_ :: Proxy t) -> SomeMod (fromInteger n :: Mod t)+{-# INLINABLE modulo #-}+infixl 7 `modulo`++liftUnOp+  :: (forall k. KnownNat k => Mod k -> Mod k)+  -> (Rational -> Rational)+  -> SomeMod+  -> SomeMod+liftUnOp fm fr = \case+  SomeMod m -> SomeMod (fm m)+  InfMod  r -> InfMod  (fr r)+{-# INLINEABLE liftUnOp #-}++liftBinOpMod+  :: (KnownNat m, KnownNat n)+  => (forall k. KnownNat k => Mod k -> Mod k -> Mod k)+  -> Mod m+  -> Mod n+  -> SomeMod+liftBinOpMod f mx@(Mod x) my@(Mod y) = case someNatVal m of+  SomeNat (_ :: Proxy t) -> SomeMod (Mod (x `mod` m) `f` Mod (y `mod` m) :: Mod t)+  where+    m = natVal mx `gcd` natVal my++liftBinOp+  :: (forall k. KnownNat k => Mod k -> Mod k -> Mod k)+  -> (Rational -> Rational -> Rational)+  -> SomeMod+  -> SomeMod+  -> SomeMod+liftBinOp _ fr (InfMod rx)  (InfMod ry)  = InfMod  (rx `fr` ry)+liftBinOp fm _ (InfMod rx)  (SomeMod my) = SomeMod (fromRational rx `fm` my)+liftBinOp fm _ (SomeMod mx) (InfMod ry)  = SomeMod (mx `fm` fromRational ry)+liftBinOp fm _ (SomeMod (mx :: Mod m)) (SomeMod (my :: Mod n))+  = case (Proxy :: Proxy m) `sameNat` (Proxy :: Proxy n) of+    Nothing   -> liftBinOpMod fm mx my+    Just Refl -> SomeMod (mx `fm` my)++instance Num SomeMod where+  (+)    = liftBinOp (+) (+)+  (-)    = liftBinOp (-) (+)+  negate = liftUnOp negate negate+  {-# INLINE negate #-}+  (*)    = liftBinOp (*) (*)+  abs    = id+  {-# INLINE abs #-}+  signum = const 1+  {-# INLINE signum #-}+  fromInteger = InfMod . fromInteger+  {-# INLINE fromInteger #-}++-- | Beware that division by residue, which is not coprime with the modulo,+-- will result in runtime error. Consider using 'invertSomeMod' instead.+instance Fractional SomeMod where+  fromRational = InfMod+  {-# INLINE fromRational #-}+  recip x = case invertSomeMod x of+    Nothing -> error $ "recip{SomeMod}: residue is not coprime with modulo"+    Just y  -> y++-- | Computes the inverse value, if it exists.+--+-- > > invertSomeMod (3 `modulo` 10)+-- > Just (7 `modulo` 10) -- because 3 * 7 = 1 :: Mod 10+-- > > invertMod (4 `modulo` 10)+-- > Nothing+-- > > invertSomeMod (fromRational (2 % 5))+-- > Just 5 % 2+invertSomeMod :: SomeMod -> Maybe SomeMod+invertSomeMod = \case+  SomeMod m -> fmap SomeMod (invertMod m)+  InfMod  r -> Just (InfMod (recip r))+{-# INLINABLE [1] invertSomeMod #-}++{-# SPECIALISE [1] powSomeMod ::+  SomeMod -> Integer -> SomeMod,+  SomeMod -> Natural -> SomeMod,+  SomeMod -> Int     -> SomeMod,+  SomeMod -> Word    -> SomeMod #-}++-- | Drop-in replacement for '^', with much better performance.+-- When -O is enabled, there is a rewrite rule, which specialises '^' to 'powSomeMod'.+--+-- > > powSomeMod (3 `modulo` 10) 4+-- > (1 `modulo` 10)+powSomeMod :: Integral a => SomeMod -> a -> SomeMod+powSomeMod (SomeMod m) a = SomeMod (m ^% a)+powSomeMod (InfMod  r) a = InfMod  (r ^  a)+{-# INLINABLE [1] powSomeMod #-}++{-# RULES "^%SomeMod" forall x p. x ^ p = powSomeMod x p #-}

@@ -0,0 +1,150 @@+-- |+-- Module:      Math.NumberTheory.Moduli.Jacobi+-- Copyright:   (c) 2011 Daniel Fischer, 2017 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Jacobi symbol.+--++{-# LANGUAGE CPP        #-}+{-# LANGUAGE LambdaCase #-}++module Math.NumberTheory.Moduli.Jacobi+  ( JacobiSymbol(..)+  , jacobi+  , jacobi'+  ) where++import Data.Array.Unboxed+import Data.Bits+import Data.Semigroup+#if __GLASGOW_HASKELL__ < 709+import Data.Word+#endif++import Math.NumberTheory.Unsafe+import Math.NumberTheory.Utils++-- | Type for result of 'jacobi'.+data JacobiSymbol = MinusOne | Zero | One+  deriving (Eq, Ord, Show)++instance Semigroup JacobiSymbol where+  (<>) = \case+    MinusOne -> negJS+    Zero     -> const Zero+    One      -> id++instance Monoid JacobiSymbol where+  mempty = One+  mappend = (<>)++negJS :: JacobiSymbol -> JacobiSymbol+negJS = \case+  MinusOne -> One+  Zero     -> Zero+  One      -> MinusOne++-- | Jacobi symbol of two numbers.+--   The \"denominator\" must be odd and positive, this condition is checked.+--+--   If both numbers have a common prime factor, the result+--   is @0@, otherwise it is &#177;1.+{-# SPECIALISE jacobi :: Integer -> Integer -> JacobiSymbol,+                         Int -> Int -> JacobiSymbol,+                         Word -> Word -> JacobiSymbol+  #-}+jacobi :: (Integral a, Bits a) => a -> a -> JacobiSymbol+jacobi a b+  | b < 0       = error "Math.NumberTheory.Moduli.jacobi: negative denominator"+  | evenI b     = error "Math.NumberTheory.Moduli.jacobi: even denominator"+  | b == 1      = One+  | otherwise   = jacobi' a b   -- b odd, > 1++-- Invariant: b > 1 and odd+-- | Jacobi symbol of two numbers without validity check of+--   the \"denominator\".+{-# SPECIALISE jacobi' :: Integer -> Integer -> JacobiSymbol,+                          Int -> Int -> JacobiSymbol,+                          Word -> Word -> JacobiSymbol+  #-}+jacobi' :: (Integral a, Bits a) => a -> a -> JacobiSymbol+jacobi' a b+  | a == 0      = Zero+  | a == 1      = One+  | a < 0       = let n | rem4 b == 1 = One+                        | otherwise   = MinusOne+                      -- Blech, minBound may pose problems+                      (z,o) = shiftToOddCount (abs $ toInteger a)+                      s | evenI z || unsafeAt jac2 (rem8 b) == 1 = n+                        | otherwise                              = negJS n+                  in s <> jacobi' (fromInteger o) b+  | a >= b      = case a `rem` b of+                    0 -> Zero+                    r -> jacPS One r b+  | evenI a     = case shiftToOddCount a of+                    (z,o) -> let r | rem4 o .&. rem4 b == 1 = One+                                   | otherwise              = MinusOne+                                 s | evenI z || unsafeAt jac2 (rem8 b) == 1 = r+                                   | otherwise                              = negJS r+                             in jacOL s b o+  | otherwise   = case rem4 a .&. rem4 b of+                    3 -> jacOL MinusOne b a+                    _ -> jacOL One      b a++-- numerator positive and smaller than denominator+{-# SPECIALISE jacPS :: JacobiSymbol -> Integer -> Integer -> JacobiSymbol,+                        JacobiSymbol -> Int -> Int -> JacobiSymbol,+                        JacobiSymbol -> Word -> Word -> JacobiSymbol+  #-}+jacPS :: (Integral a, Bits a) => JacobiSymbol -> a -> a -> JacobiSymbol+jacPS j a b+  | evenI a     = case shiftToOddCount a of+                    (z,o) | evenI z || unsafeAt jac2 (rem8 b) == 1 ->+                              jacOL (if rem4 o .&. rem4 b == 3 then (negJS j) else j) b o+                          | otherwise ->+                              jacOL (if rem4 o .&. rem4 b == 3 then j else (negJS j)) b o+  | otherwise   = jacOL (if rem4 a .&. rem4 b == 3 then (negJS j) else j) b a++-- numerator odd, positive and larger than denominator+{-# SPECIALISE jacOL :: JacobiSymbol -> Integer -> Integer -> JacobiSymbol,+                        JacobiSymbol -> Int -> Int -> JacobiSymbol,+                        JacobiSymbol -> Word -> Word -> JacobiSymbol+  #-}+jacOL :: (Integral a, Bits a) => JacobiSymbol -> a -> a -> JacobiSymbol+jacOL j a b+  | b == 1    = j+  | otherwise = case a `rem` b of+                 0 -> Zero+                 r -> jacPS j r b++-- Utilities++-- For large Integers, going via Int is much faster than bit-fiddling+-- on the Integer, so we do that.+{-# SPECIALISE evenI :: Integer -> Bool,+                        Int -> Bool,+                        Word -> Bool+  #-}+evenI :: Integral a => a -> Bool+evenI n = fromIntegral n .&. 1 == (0 :: Int)++{-# SPECIALISE rem4 :: Integer -> Int,+                       Int -> Int,+                       Word -> Int+  #-}+rem4 :: Integral a => a -> Int+rem4 n = fromIntegral n .&. 3++{-# SPECIALISE rem8 :: Integer -> Int,+                       Int -> Int,+                       Word -> Int+  #-}+rem8 :: Integral a => a -> Int+rem8 n = fromIntegral n .&. 7++jac2 :: UArray Int Int+jac2 = array (0,7) [(0,0),(1,1),(2,0),(3,-1),(4,0),(5,-1),(6,0),(7,1)]

@@ -0,0 +1,229 @@+-- |+-- Module:      Math.NumberTheory.Moduli.Sqrt+-- Copyright:   (c) 2011 Daniel Fischer+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+-- Portability: Non-portable (GHC extensions)+--+-- Modular square roots.+--++{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE CPP          #-}++module Math.NumberTheory.Moduli.Sqrt+  ( sqrtModP+  , sqrtModPList+  , sqrtModP'+  , tonelliShanks+  , sqrtModPP+  , sqrtModPPList+  , sqrtModF+  , sqrtModFList+  ) where++import Control.Monad (liftM2)+import Data.Bits+import Data.List (nub)+#if __GLASGOW_HASKELL__ < 709+import Data.Word+#endif+import GHC.Integer.GMP.Internals++import Math.NumberTheory.Moduli.Chinese+import Math.NumberTheory.Moduli.Jacobi+import Math.NumberTheory.Primes.Sieve (sieveFrom)+import Math.NumberTheory.Utils (shiftToOddCount, splitOff)++-- | @sqrtModP n prime@ calculates a modular square root of @n@ modulo @prime@+--   if that exists. The second argument /must/ be a (positive) prime, otherwise+--   the computation may not terminate and if it does, may yield a wrong result.+--   The precondition is /not/ checked.+--+--   If @prime@ is a prime and @n@ a quadratic residue modulo @prime@, the result+--   is @Just r@ where @r^2 ≡ n (mod prime)@, if @n@ is a quadratic nonresidue,+--   the result is @Nothing@.+sqrtModP :: Integer -> Integer -> Maybe Integer+sqrtModP n 2 = Just (n `mod` 2)+sqrtModP n prime = case jacobi' n prime of+                     MinusOne -> Nothing+                     Zero     -> Just 0+                     One      -> Just (sqrtModP' (n `mod` prime) prime)++-- | @sqrtModPList n prime@ computes the list of all square roots of @n@+--   modulo @prime@. @prime@ /must/ be a (positive) prime.+--   The precondition is /not/ checked.+sqrtModPList :: Integer -> Integer -> [Integer]+sqrtModPList n prime+    | prime == 2    = [n `mod` 2]+    | otherwise     = case sqrtModP n prime of+                        Just 0 -> [0]+                        Just r -> [r,prime-r] -- The group of units in Z/(p) is cyclic+                        _      -> []++-- | @sqrtModP' square prime@ finds a square root of @square@ modulo+--   prime. @prime@ /must/ be a (positive) prime, and @square@ /must/ be a positive+--   quadratic residue modulo @prime@, i.e. @'jacobi square prime == 1@.+--   The precondition is /not/ checked.+sqrtModP' :: Integer -> Integer -> Integer+sqrtModP' square prime+    | prime == 2    = square+    | rem4 prime == 3 = powModInteger square ((prime + 1) `quot` 4) prime+    | otherwise     = tonelliShanks square prime++-- | @tonelliShanks square prime@ calculates a square root of @square@+--   modulo @prime@, where @prime@ is a prime of the form @4*k + 1@ and+--   @square@ is a positive quadratic residue modulo @prime@, using the+--   Tonelli-Shanks algorithm.+--   No checks on the input are performed.+tonelliShanks :: Integer -> Integer -> Integer+tonelliShanks square prime = loop rc t1 generator log2+  where+    (log2,q) = shiftToOddCount (prime-1)+    nonSquare = findNonSquare prime+    generator = powModInteger nonSquare q prime+    rc = powModInteger square ((q+1) `quot` 2) prime+    t1 = powModInteger square q prime+    msqr x = (x*x) `rem` prime+    msquare 0 x = x+    msquare k x = msquare (k-1) (msqr x)+    findPeriod per 1 = per+    findPeriod per x = findPeriod (per+1) (msqr x)+    loop !r t c m+        | t == 1    = r+        | otherwise = loop nextR nextT nextC nextM+          where+            nextM = findPeriod 0 t+            b     = msquare (m - 1 - nextM) c+            nextR = (r*b) `rem` prime+            nextC = msqr b+            nextT = (t*nextC) `rem` prime++-- | @sqrtModPP n (prime,expo)@ calculates a square root of @n@+--   modulo @prime^expo@ if one exists. @prime@ /must/ be a+--   (positive) prime. @expo@ must be positive, @n@ must be coprime+--   to @prime@+sqrtModPP :: Integer -> (Integer,Int) -> Maybe Integer+sqrtModPP n (2,e) = sqM2P n e+sqrtModPP n (prime,expo) = case sqrtModP n prime of+                             Just r -> fixup r+                             _      -> Nothing+  where+    fixup r = let diff' = r*r-n+              in if diff' == 0+                   then Just r+                   else case splitOff prime diff' of+                          (e,q) | expo <= e -> Just r+                                | otherwise -> fmap (\inv -> hoist inv r (q `mod` prime) (prime^e)) (recipMod (2*r) prime)++    hoist inv root elim pp+        | diff' == 0    = root'+        | expo <= ex    = root'+        | otherwise     = hoist inv root' (nelim `mod` prime) (prime^ex)+          where+            root' = (root + (inv*(prime-elim))*pp) `mod` (prime*pp)+            diff' = root'*root' - n+            (ex, nelim) = splitOff prime diff'++-- dirty, dirty+sqM2P :: Integer -> Int -> Maybe Integer+sqM2P n e+    | e < 2     = Just (n `mod` 2)+    | n' == 0   = Just 0+    | odd k     = Nothing+    | otherwise = fmap ((`mod` mdl) . (`shiftL` k2)) $ solve s e2+      where+        mdl = 1 `shiftL` e+        n' = n `mod` mdl+        (k,s) = shiftToOddCount n'+        k2 = k `quot` 2+        e2 = e-k+        solve _ 1 = Just 1+        solve 1 _ = Just 1+        solve r _+            | rem4 r == 3   = Nothing  -- otherwise r ≡ 1 (mod 4)+            | rem8 r == 5   = Nothing  -- otherwise r ≡ 1 (mod 8)+            | otherwise     = fixup r (fst $ shiftToOddCount (r-1))+              where+                fixup x pw+                    | pw >= e2  = Just x+                    | otherwise = fixup x' pw'+                      where+                        x' = x + (1 `shiftL` (pw-1))+                        d = x'*x' - r+                        pw' = if d == 0 then e2 else fst (shiftToOddCount d)++-- | @sqrtModF n primePowers@ calculates a square root of @n@ modulo+--   @product [p^k | (p,k) <- primePowers]@ if one exists and all primes+--   are distinct.+--   The list must be non-empty, @n@ must be coprime with all primes.+sqrtModF :: Integer -> [(Integer,Int)] -> Maybe Integer+sqrtModF _ []  = Nothing+sqrtModF n pps = do roots <- mapM (sqrtModPP n) pps+                    chineseRemainder $ zip roots (map (uncurry (^)) pps)++-- | @sqrtModFList n primePowers@ calculates all square roots of @n@ modulo+--   @product [p^k | (p,k) <- primePowers]@ if all primes are distinct.+--   The list must be non-empty, @n@ must be coprime with all primes.+sqrtModFList :: Integer -> [(Integer,Int)] -> [Integer]+sqrtModFList _ []  = []+sqrtModFList n pps = map fst $ foldl1 (liftM2 comb) cs+  where+    ms :: [Integer]+    ms = map (uncurry (^)) pps+    rs :: [[Integer]]+    rs = map (sqrtModPPList n) pps+    cs :: [[(Integer,Integer)]]+    cs = zipWith (\l m -> map (\x -> (x,m)) l) rs ms+    comb t1@(_,m1) t2@(_,m2) = (chineseRemainder2 t1 t2,m1*m2)++-- | @sqrtModPPList n (prime,expo)@ calculates the list of all+--   square roots of @n@ modulo @prime^expo@. The same restriction+--   as in 'sqrtModPP' applies to the arguments.+sqrtModPPList :: Integer -> (Integer,Int) -> [Integer]+sqrtModPPList n (2,1) = [n `mod` 2]+sqrtModPPList n (2,expo)+    = case sqM2P n expo of+        Just r -> let m = 1 `shiftL` (expo-1)+                  in nub [r, (r+m) `mod` (2*m), (m-r) `mod` (2*m), 2*m-r]+        _ -> []+sqrtModPPList n pe@(prime,expo)+    = case sqrtModPP n pe of+        Just 0 -> [0]+        Just r -> [prime^expo - r, r] -- The group of units in Z/(p^e) is cyclic+        _      -> []+++-- Utilities++{-# SPECIALISE rem4 :: Integer -> Int,+                       Int -> Int,+                       Word -> Int+  #-}+rem4 :: Integral a => a -> Int+rem4 n = fromIntegral n .&. 3++{-# SPECIALISE rem8 :: Integer -> Int,+                       Int -> Int,+                       Word -> Int+  #-}+rem8 :: Integral a => a -> Int+rem8 n = fromIntegral n .&. 7++findNonSquare :: Integer -> Integer+findNonSquare n+    | rem8 n == 5 || rem8 n == 3  = 2+    | otherwise = search primelist+      where+        primelist = [3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67]+                        ++ sieveFrom (68 + n `rem` 4) -- prevent sharing+        search (p:ps) = case jacobi' p n of+          MinusOne -> p+          _        -> search ps+        search _ = error "Should never have happened, prime list exhausted."++recipMod :: Integer -> Integer -> Maybe Integer+recipMod x m = case recipModInteger x m of+  0 -> Nothing+  y -> Just y

@@ -23,11 +23,14 @@  -- | @totientSum n@ is, for @n > 0@, the sum of @[totient k | k <- [1 .. n]]@, --   computed via generalised Moebius inversion.---   Arguments less than 1 cause an error to be raised.+--   See <http://mathworld.wolfram.com/TotientSummatoryFunction.html> for the+--   formula used for @totientSum@. totientSum :: Int -> Integer-totientSum = (+1) . generalInversion (triangle . fromIntegral)+totientSum n+  | n < 1 = 0+  | otherwise = generalInversion (triangle . fromIntegral) n   where-    triangle n = (n*(n-1)) `quot` 2+    triangle k = (k*(k+1)) `quot` 2  -- | @generalInversion g n@ evaluates the generalised Moebius inversion of @g@ --   at the argument @n@.

@@ -24,11 +24,14 @@  -- | @totientSum n@ is, for @n > 0@, the sum of @[totient k | k <- [1 .. n]]@, --   computed via generalised Moebius inversion.---   Arguments less than 1 cause an error to be raised.+--   See <http://mathworld.wolfram.com/TotientSummatoryFunction.html> for the+--   formula used for @totientSum@. totientSum :: Int -> Int-totientSum = (+1) . generalInversion triangle+totientSum n+  | n < 1 = 0+  | otherwise = generalInversion (triangle . fromIntegral) n   where-    triangle n = (n*(n-1)) `quot` 2+    triangle k = (k*(k+1)) `quot` 2  -- | @generalInversion g n@ evaluates the generalised Moebius inversion of @g@ --   at the argument @n@.

@@ -32,12 +32,9 @@   , exactRoot   , isPerfectPower   , highestPower-    -- Modular powers-  , powerMod   ) where  import Math.NumberTheory.Powers.Squares import Math.NumberTheory.Powers.Cubes import Math.NumberTheory.Powers.Fourth import Math.NumberTheory.Powers.General-import Math.NumberTheory.Moduli

@@ -25,6 +25,7 @@ import Data.Maybe import Data.Bits +import Math.NumberTheory.Moduli.Class import Math.NumberTheory.Primes.Factorisation.Montgomery import Math.NumberTheory.Primes.Testing.Certificates.Internal import Math.NumberTheory.Primes.Testing.Probabilistic@@ -118,9 +119,11 @@             | count < 0 = return ([],[(m,1)])             | otherwise = do                 s <- rndR m-                case montgomeryFactorisation m b1 b2 s of-                  Nothing -> repFact m b1 b2 (count-1)-                  Just d  -> do+                case s `modulo` fromInteger m of+                  InfMod{} -> error "impossible case"+                  SomeMod sm -> case montgomeryFactorisation b1 b2 sm of+                    Nothing -> repFact m b1 b2 (count-1)+                    Just d  -> do                       let !cof = m `quot` d                       case gcd cof d of                         1 -> do

@@ -20,8 +20,16 @@ -- -- Given enough time, the algorithm should be able to factor numbers of 100-120 digits, but it -- is best suited for numbers of up to 50-60 digits.-{-# LANGUAGE CPP, BangPatterns, MagicHash #-}++{-# LANGUAGE BangPatterns   #-}+{-# LANGUAGE CPP            #-}+{-# LANGUAGE DataKinds      #-}+{-# LANGUAGE KindSignatures #-}+{-# LANGUAGE LambdaCase     #-}++{-# OPTIONS_GHC -fno-warn-type-defaults #-} {-# OPTIONS_HADDOCK hide #-}+ module Math.NumberTheory.Primes.Factorisation.Montgomery   ( -- *  Complete factorisation functions     -- ** Functions with input checking@@ -42,19 +50,22 @@  #include "MachDeps.h" -import GHC.Base- import System.Random import Control.Monad.State.Strict #if __GLASGOW_HASKELL__ < 709 import Control.Applicative+import Data.Word #endif import Data.Bits+import Data.IntMap (IntMap)+import qualified Data.IntMap as IM+import Data.List (foldl') import Data.Maybe -import GHC.Integer.Logarithms+import GHC.TypeNats.Compat -import Math.NumberTheory.Logarithms+import Math.NumberTheory.Curves.Montgomery+import Math.NumberTheory.Moduli.Class import Math.NumberTheory.Powers.General     (highestPower, largePFPower) import Math.NumberTheory.Powers.Squares     (integerSquareRoot') import Math.NumberTheory.Primes.Sieve.Eratosthenes@@ -169,7 +180,7 @@                            then return pfs                            else do                                nfs <- forM cfs $ \(k,j) ->-                                   mult j <$> fact k (if null pfs then digs+4 else digs)+                                   mult j <$> fact k (if null pfs then digs+5 else digs)                                return (mergeAll $ pfs:nfs)         repFact m b1 b2 count = case perfPw m of                                   (_,1) -> workFact m b1 b2 count@@ -179,12 +190,14 @@                                       (as,bs) <- workFact b b1 b2 count                                       return $ (mult e as, mult e bs)         workFact m b1 b2 count-            | count < 0 = return ([],[(m,1)])+            | count == 0 = return ([],[(m,1)])             | otherwise = do                 s <- rndR m-                case montgomeryFactorisation m b1 b2 s of-                  Nothing -> workFact m b1 b2 (count-1)-                  Just d  -> do+                case s `modulo` fromInteger m of+                  InfMod{} -> error "impossible case"+                  SomeMod sm -> case montgomeryFactorisation b1 b2 sm of+                    Nothing -> workFact m b1 b2 (count-1)+                    Just d  -> do                       let !cof = m `quot` d                       case gcd cof d of                         1 -> do@@ -226,132 +239,71 @@ --   It is assumed that @n@ has no small prime factors. -- --   The result is maybe a nontrivial divisor of @n@.-montgomeryFactorisation :: Integer -> Word -> Word -> Integer -> Maybe Integer-montgomeryFactorisation n b1 b2 s = go p5 (list primeStore)+montgomeryFactorisation :: KnownNat n => Word -> Word -> Mod n -> Maybe Integer+montgomeryFactorisation b1 b2 s = case newPoint (getVal s) n of+  Nothing             -> Nothing+  Just (SomePoint p0) -> do+    -- Small step: for each prime p <= b1+    -- multiply point 'p0' by the highest power p^k <= b1.+    let q = foldl (flip multiply) p0 smallPowers+        z = pointZ q++    fromIntegral <$> case gcd n z of+      -- If small step did not succeed, perform a big step.+      1 -> case gcd n (bigStep q b1 b2) of+        1 -> Nothing+        g -> Just g+      g -> Just g   where-    l2 = wordLog2' b1-    b1i = toInteger b1-    (^~) :: Word -> Int -> Word-    w ^~ i = w ^ i-    (e, p0) = montgomeryData n s-    dbl pt = double n e pt-    dbln 0 !pt = pt-    dbln k pt = dbln (k-1) (dbl pt)-    p2 = dbln l2 p0-#if WORD_SIZE_IN_BITS == 64-    mul a b c = (a*b) `quot` c       -- can't overflow, work on Int-#else-    mul a b c = fromInteger ((toInteger a * b) `quot` c) -- might overflow if Int is used-#endif-    adjust bd ml w-      | w <= bd     = adjust bd ml (w*ml)-      | otherwise   = w-    l3 = mul l2 190537 301994-    w3 = 3 ^~ l3-    pw3 = adjust (b1 `quot` 3) 3 w3-    p3 = multiply n e pw3 p2-    l5 = mul l2 1936274 4495889-    w5 = 5 ^~ l5-    pw5 = adjust (b1 `quot` 5) 5 w5-    p5 = multiply n e pw5 p3-    go (P _ 0) _ = Nothing-    go !pt@(P _ z) (pr:prs)-      | pr <= b1    = let !lp = integerLogBase' (fromIntegral pr) b1i-                      in go (multiply n e (pr ^~ lp) pt) prs-      | otherwise   = case gcd n z of-                        1 -> lgo (multiply n e pr pt) prs-                        g -> Just g-    go (P _ z) _    = case gcd n z of-                        1 -> Nothing-                        g -> Just g-    lgo (P _ 0) _ = Nothing-    lgo !pt@(P _ z) (pr:prs)-      | pr <= b2    = lgo (multiply n e pr pt) prs-      | otherwise   = case gcd n z of-                        1 -> Nothing-                        g -> Just g-    lgo (P _ z) _   = case gcd n z of-                        1 -> Nothing-                        g -> Just g+    n = getMod s+    smallPrimes = takeWhile (<= b1) (2 : 3 : 5 : list primeStore)+    smallPowers = map findPower smallPrimes+    findPower p = go p+      where+        go acc+          | acc <= b1 `quot` p = go (acc * p)+          | otherwise          = acc --------------------------------------------------------------------------------------------------------                            Helpers, Curves and elliptic arithmetics                            -------------------------------------------------------------------------------------------------------+-- | The implementation follows the algorithm at p. 6-7+-- of <http://www.hyperelliptic.org/tanja/SHARCS/talks06/Gaj.pdf Implementing the Elliptic Curve Method of Factoring in Reconfigurable Hardware>+-- by K. Gaj, S. Kwon et al.+bigStep :: (KnownNat a24, KnownNat n) => Point a24 n -> Word -> Word -> Integer+bigStep q b1 b2 = rs+  where+    n = pointN q --- A Montgomery curve is given by y^2 = x^3 + (A_n / A_d - 2)*x^2 + x (mod n).--- We store A_n and 4*A_d, since A_n occurs with the factor 4 in all formulae.-data Curve = C !Integer !Integer+    b0 = b1 - b1 `rem` wheel+    qks = zip [0..] $ map (\k -> multiply k q) wheelCoprimes+    qs = enumAndMultiplyFromThenTo q b0 (b0 + wheel) b2 --- Point in the projective plane, will be on the curve--- A coordinate transformation eliminates the y-coordinate, hence--- we store only x and z-data Point = P !Integer !Integer+    rs = foldl' (\ts (_cHi, p) -> foldl' (\us (_cLo, pq) ->+        us * (pointZ p * pointX pq - pointX p * pointZ pq) `rem` n+        ) ts qks) 1 qs --- Get curve and point to start--- Input should satisfy 6 <= s < n-1-montgomeryData :: Integer -> Integer -> (Curve, Point)-montgomeryData n s = (C an ad4, P x z)-  where-    u = (s*s-5) `mod` n-    v = (4*s) `mod` n-    d = (v-u)-    x = (u*u*u) `mod` n-    z = (v*v*v) `mod` n-    an = ((d*d)*(d*(3*u+v))) `mod` n-    ad4 = (16*x*v) `mod` n+wheel :: Word+wheel = 210 --- Addition on the curve, given the modulus n and three points,--- p0, p1 and p2, with p0 = p2 - p1, calculate the point p1 + p2.--- Note that the addition does not depend on the curve.-add :: Integer -> Point -> Point -> Point -> Point-add n (P x0 z0) (P x1 z1) (P x2 z2) = P x3 z3-  where-    a = (x1-z1)*(x2+z2)-    b = (x1+z1)*(x2-z2)-    c = a+b-    d = a-b-    x3 = (c*c*z0) `rem` n-    z3 = (d*d*x0) `rem` n+wheelCoprimes :: [Word]+wheelCoprimes = [ k | k <- [1 .. wheel `div` 2], k `gcd` wheel == 1 ] --- Double a point on the curve.-double :: Integer -> Curve -> Point -> Point-double n (C an ad4) (P x z) = P x' z'+-- | Same as map (id *** flip multiply p) [from, thn .. to],+-- but calculated in more efficient way.+enumAndMultiplyFromThenTo+  :: (KnownNat a24, KnownNat n)+  => Point a24 n+  -> Word+  -> Word+  -> Word+  -> [(Word, Point a24 n)]+enumAndMultiplyFromThenTo p from thn to = zip [from, thn .. to] progression   where-    r = x+z-    s = x-z-    u = r*r-    v = s*s-    t = u-v-    x' = (ad4*u*v) `rem` n-    z' = ((ad4*v+t*an)*t) `rem` n+    step = thn - from --- Multiply a point on the curve by a Word.--- Within Word range, we can use the faster variant going--- from high-order bits to low-order.-multiply :: Integer -> Curve -> Word -> Point -> Point-multiply n cve (W# w##) p =-    case wordLog2# w## of-      l# -> go (l# -# 1#) p (double n cve p)-  where-    go 0# !p0 !p1 = case w## `and#` 1## of-                      0## -> double n cve p0-                      _   -> add n p p0 p1-    go i# p0 p1 = case (uncheckedShiftRL# w## i#) `and#` 1## of-                    0## -> go (i# -# 1#) (double n cve p0) (add n p p0 p1)-                    _   -> go (i# -# 1#) (add n p p0 p1) (double n cve p1)+    pFrom = multiply from p+    pThen = multiply thn  p+    pStep = multiply step p -{-  Not (yet) needed--- Multiply a point on the curve by an Integer.-multIgr :: Integer -> Curve -> Integer -> Point -> Point-multIgr n cve k p = go k-  where-    go 1 = (p, double n cve p)-    go m = case m `quotRem` 2 of-             (q,r) -> let !(!s, l) = go q-                      in case r of-                           0 -> (double n cve s, add n p s l)-                           _ -> (add n p s l, double n cve l)--}+    progression = pFrom : pThen : zipWith (\x0 x1 -> add x0 pStep x1) progression (tail progression)  -- primes, compactly stored as a bit sieve primeStore :: [PrimeSieve]@@ -381,34 +333,41 @@     go m [] = ([(m,1)], Nothing)  -- helpers: merge sorted lists-merge :: [(Integer,Int)] -> [(Integer,Int)] -> [(Integer,Int)]-merge xxs@(x@(p,k):xs) yys@(y@(q,m):ys) = case compare p q of-                                            LT -> x : merge xs yys-                                            EQ -> (p,k+m) : merge xs ys-                                            GT -> y : merge xxs ys+merge :: (Ord a, Num b) => [(a, b)] -> [(a, b)] -> [(a, b)] merge xs [] = xs-merge _ ys = ys+merge [] ys = ys+merge xxs@(x@(p, k) : xs) yys@(y@(q, m) : ys)+  = case p `compare` q of+    LT -> x          : merge xs yys+    EQ -> (p, k + m) : merge xs  ys+    GT -> y          : merge xxs ys -mergeAll :: [[(Integer,Int)]] -> [(Integer,Int)]-mergeAll [] = []-mergeAll [xs] = xs-mergeAll (xs:ys:zss) = merge (merge xs ys) (mergeAll zss)+mergeAll :: (Ord a, Num b) => [[(a, b)]] -> [(a, b)]+mergeAll = \case+  []              -> []+  [xs]            -> xs+  (xs : ys : zss) -> merge (merge xs ys) (mergeAll zss) --- Parameters for the factorisation, the two b-parameters for montgomery and the number of tries--- to use these, depending on the size of the factor we are looking for.--- The numbers are roughly based on the parameters listed on Dario Alpern's ECM site.-testParms :: [(Int,Word,Word,Int)]-testParms = [ (12, 400, 10000, 10), (15, 2000, 50000, 25), (20, 11000, 150000, 90)-            , (25, 50000, 500000, 300), (30, 250000, 1500000, 700)-            , (35, 1000000, 4000000, 1800), (40, 3000000, 12000000, 5100)-            , (45, 11000000, 45000000, 10600), (50, 43000000, 200000000, 19300)-            , (55, 80000000, 400000000,30000), (60, 120000000, 700000000, 50000)-            ]+-- | For a given estimated decimal length of the smallest prime factor+-- ("tier") return parameters B1, B2 and the number of curves to try+-- before next "tier".+-- Roughly based on http://www.mersennewiki.org/index.php/Elliptic_Curve_Method#Choosing_the_best_parameters_for_ECM+testParms :: IntMap (Word, Word, Word)+testParms = IM.fromList+  [ (12, (       400,        40000,     10))+  , (15, (      2000,       200000,     25))+  , (20, (     11000,      1100000,     90))+  , (25, (     50000,      5000000,    300))+  , (30, (    250000,     25000000,    700))+  , (35, (   1000000,    100000000,   1800))+  , (40, (   3000000,    300000000,   5100))+  , (45, (  11000000,   1100000000,  10600))+  , (50, (  43000000,   4300000000,  19300))+  , (55, ( 110000000,  11000000000,  49000))+  , (60, ( 260000000,  26000000000, 124000))+  , (65, ( 850000000,  85000000000, 210000))+  , (70, (2900000000, 290000000000, 340000))+  ] -findParms :: Int -> (Word, Word, Int)-findParms digs = go (100, 1000, 7) testParms-  where-    go p ((d,b1,b2,ct):rest)-      | digs < d    = p-      | otherwise   = go (b1,b2,ct) rest-    go p [] = p+findParms :: Int -> (Word, Word, Word)+findParms digs = maybe (wheel, 1000, 7) snd (IM.lookupLT digs testParms)

@@ -35,8 +35,9 @@ #endif import Data.Bits import Data.Maybe+import GHC.Integer.GMP.Internals -import Math.NumberTheory.Moduli+import Math.NumberTheory.Moduli.Class import Math.NumberTheory.Utils import Math.NumberTheory.Primes.Factorisation.TrialDivision import Math.NumberTheory.Primes.Factorisation.Montgomery@@ -201,8 +202,8 @@     verify (pf,_,base,proof) = pf == cprime proof && crit pf base && checkPrimalityProof proof     crit pf base = gcd p (x-1) == 1 && y == 1       where-        x = powerModInteger' base (pm1 `quot` pf) p-        y = powerModInteger' x pf p+        x = powModInteger base (pm1 `quot` pf) p+        y = powModInteger x pf p  -- | @'trivial'@ records a trivially known prime. --   If the argument is not one of them, an error is raised.@@ -290,8 +291,8 @@                             Prime ppr ->(p,e,bs,ppr)               where                 q = nm1 `quot` p-                x = powerModInteger' bs q n-                y = powerModInteger' x p n+                x = powModInteger bs q n+                y = powModInteger x p n                 g = gcd n (x-1)  -- | Find a decomposition of p-1 for the pocklington certificate.@@ -327,13 +328,15 @@     (lo,hi,count) = findParms digits  -- | Find a factor or say with which curve to continue.-findLoop :: Integer -> Word -> Word -> Int -> Integer -> Either Integer Integer+findLoop :: Integer -> Word -> Word -> Word -> Integer -> Either Integer Integer findLoop _ _  _  0  s = Left s findLoop n lo hi ct s     | n <= s+2  = Left 6-    | otherwise = case montgomeryFactorisation n lo hi s of-                    Nothing -> findLoop n lo hi (ct-1) (s+1)-                    Just fct+    | otherwise = case s `modulo` fromInteger n of+                    InfMod{}   -> error "impossible case"+                    SomeMod sn -> case montgomeryFactorisation lo hi sn of+                      Nothing -> findLoop n lo hi (ct-1) (s+1)+                      Just fct                         | bailliePSW fct -> Right fct                         | otherwise -> Right (findFactor fct 8 (s+1)) 

@@ -23,8 +23,10 @@ import Data.Bits import GHC.Base import GHC.Integer.GMP.Internals+import GHC.TypeNats.Compat -import Math.NumberTheory.Moduli+import Math.NumberTheory.Moduli.Class+import Math.NumberTheory.Moduli.Jacobi import Math.NumberTheory.Utils import Math.NumberTheory.Powers.Squares @@ -83,14 +85,19 @@   | n < 0          = error "isStrongFermatPP: negative argument"   | n <= 1         = False   | n == 2         = True-  | b `mod` n == 0 = True-  | otherwise      = a == 1 || go t a+  | otherwise      = case b `modulo` fromInteger n of+                       SomeMod b' -> isStrongFermatPPMod b'+                       InfMod{}   -> True++isStrongFermatPPMod :: KnownNat n => Mod n -> Bool+isStrongFermatPPMod b = b == 0 || a == 1 || go t a   where-    m = n-1-    (t,u) = shiftToOddCount m-    a = powerModInteger' (b `mod` n) u n+    m = -1+    (t, u) = shiftToOddCount $ getVal m+    a = b ^% u+     go 0 _ = False-    go k x = x == m || go (k-1) ((x*x) `rem` n)+    go k x = x == m || go (k - 1) (x * x)  -- | @'isFermatPP' n b@ tests whether @n@ is a Fermat probable prime --   for the base @b@, that is, whether @b^(n-1) `mod` n == 1@.@@ -109,7 +116,9 @@ --   of prime bases is reasonable to find out whether it's worth the --   effort to undertake the prime factorisation). isFermatPP :: Integer -> Integer -> Bool-isFermatPP n b = powerModInteger' b (n-1) n == 1+isFermatPP n b = case b `modulo` fromInteger n of+  SomeMod b' -> b' ^% (n-1) == 1+  InfMod{}   -> True  -- | Primality test after Baillie, Pomerance, Selfridge and Wagstaff. --   The Baillie-PSW test consists of a strong Fermat probable primality@@ -148,9 +157,9 @@       r = integerSquareRoot n       d = find True 5       find !pos cd = case jacobi' (n `rem` cd) cd of-                       0 -> if cd == n then 1 else 0-                       1 -> find (not pos) (cd+2)-                       _ -> if pos then cd else (-cd)+                       MinusOne -> if pos then cd else (-cd)+                       Zero     -> if cd == n then 1 else 0+                       One      -> find (not pos) (cd+2)       q = (1-d) `quot` 4       (t,o) = shiftToOddCount (n+1)       (uo, vo, qo) = testLucas n q o

@@ -49,15 +49,19 @@ --   close proximity, it is better to use the simple addition --   formula starting from an appropriate pair of successive --   Fibonacci numbers.-fibonacci :: Int -> Integer+fibonacci :: Num a => Int -> a fibonacci = fst . fibonacciPair+{-# SPECIALIZE fibonacci :: Int -> Int     #-}+{-# SPECIALIZE fibonacci :: Int -> Word    #-}+{-# SPECIALIZE fibonacci :: Int -> Integer #-}+{-# SPECIALIZE fibonacci :: Int -> Natural #-}  -- | @'fibonacciPair' k@ returns the pair @(F(k), F(k+1))@ of the @k@-th --   Fibonacci number and its successor, thus it can be used to calculate --   the Fibonacci numbers from some index on without needing to compute --   the previous. The pair is efficiently calculated --   in /O/(@log (abs k)@) steps. The index may be negative.-fibonacciPair :: Int -> (Integer, Integer)+fibonacciPair :: Num a => Int -> (a, a) fibonacciPair n   | n < 0     = let (f,g) = fibonacciPair (-(n+1)) in if testBit n 0 then (g, -f) else (-g, f)   | n == 0    = (0, 1)@@ -68,17 +72,25 @@         | otherwise   = look (k-1)       go k g f         | k < 0       = (f, f+g)-        | testBit n k = go (k-1) (f*(f+shiftL g 1)) ((f+g)*shiftL f 1 + g*g)-        | otherwise   = go (k-1) (f*f+g*g) (f*(f+shiftL g 1))+        | testBit n k = go (k-1) (f*(f+shiftL1 g)) ((f+g)*shiftL1 f + g*g)+        | otherwise   = go (k-1) (f*f+g*g) (f*(f+shiftL1 g))+{-# SPECIALIZE fibonacciPair :: Int -> (Int, Int)         #-}+{-# SPECIALIZE fibonacciPair :: Int -> (Word, Word)       #-}+{-# SPECIALIZE fibonacciPair :: Int -> (Integer, Integer) #-}+{-# SPECIALIZE fibonacciPair :: Int -> (Natural, Natural) #-}  -- | @'lucas' k@ computes the @k@-th Lucas number. Very similar --   to @'fibonacci'@.-lucas :: Int -> Integer+lucas :: Num a => Int -> a lucas = fst . lucasPair+{-# SPECIALIZE lucas :: Int -> Int     #-}+{-# SPECIALIZE lucas :: Int -> Word    #-}+{-# SPECIALIZE lucas :: Int -> Integer #-}+{-# SPECIALIZE lucas :: Int -> Natural #-}  -- | @'lucasPair' k@ computes the pair @(L(k), L(k+1))@ of the @k@-th --   Lucas number and its successor. Very similar to @'fibonacciPair'@.-lucasPair :: Int -> (Integer, Integer)+lucasPair :: Num a => Int -> (a, a) lucasPair n   | n < 0     = let (f,g) = lucasPair (-(n+1)) in if testBit n 0 then (-g, f) else (g, -f)   | n == 0    = (2, 1)@@ -88,13 +100,16 @@         | testBit n k = go (k-1) 0 1         | otherwise   = look (k-1)       go k g f-        | k < 0       = (shiftL g 1 + f,g+3*f)+        | k < 0       = (shiftL1 g + f,g+3*f)         | otherwise   = go (k-1) g' f'           where             (f',g')-              | testBit n k = (shiftL (f*(f+g)) 1 + g*g,f*(shiftL g 1 + f))-              | otherwise   = (f*(shiftL g 1 + f),f*f+g*g)-+              | testBit n k = (shiftL1 (f*(f+g)) + g*g,f*(shiftL1 g + f))+              | otherwise   = (f*(shiftL1 g + f),f*f+g*g)+{-# SPECIALIZE lucasPair :: Int -> (Int, Int)         #-}+{-# SPECIALIZE lucasPair :: Int -> (Word, Word)       #-}+{-# SPECIALIZE lucasPair :: Int -> (Integer, Integer) #-}+{-# SPECIALIZE lucasPair :: Int -> (Natural, Natural) #-}  -- | @'generalLucas' p q k@ calculates the quadruple @(U(k), U(k+1), V(k), V(k+1))@ --   where @U(i)@ is the Lucas sequence of the first kind and @V(i)@ the Lucas@@ -106,7 +121,7 @@ --   the second kind for these parameters. --   Here, the index must be non-negative, since the terms of the sequence for --   negative indices are in general not integers.-generalLucas :: Integer -> Integer -> Int -> (Integer, Integer, Integer, Integer)+generalLucas :: Num a => a -> a -> Int -> (a, a, a, a) generalLucas p q k   | k < 0       = error "generalLucas: negative index"   | k == 0      = (0,1,2,p)@@ -116,6 +131,13 @@         | testBit k i   = go (i-1) 1 p p q         | otherwise     = look (i-1)       go i un un1 vn qn-        | i < 0         = (un, un1, vn, p*un1 - shiftL (q*un) 1)+        | i < 0         = (un, un1, vn, p*un1 - shiftL1 (q*un))         | testBit k i   = go (i-1) (un1*vn-qn) ((p*un1-q*un)*vn - p*qn) ((p*un1 - (2*q)*un)*vn - p*qn) (qn*qn*q)         | otherwise     = go (i-1) (un*vn) (un1*vn-qn) (vn*vn - 2*qn) (qn*qn)+{-# SPECIALIZE generalLucas :: Int     -> Int     -> Int -> (Int, Int, Int, Int)                 #-}+{-# SPECIALIZE generalLucas :: Word    -> Word    -> Int -> (Word, Word, Word, Word)             #-}+{-# SPECIALIZE generalLucas :: Integer -> Integer -> Int -> (Integer, Integer, Integer, Integer) #-}+{-# SPECIALIZE generalLucas :: Natural -> Natural -> Int -> (Natural, Natural, Natural, Natural) #-}++shiftL1 :: Num a => a -> a+shiftL1 n = n + n

@@ -1,5 +1,5 @@ name                : arithmoi-version             : 0.5.0.1+version             : 0.6.0.0 cabal-version       : >= 1.10 author              : Daniel Fischer copyright           : (c) 2011 Daniel Fischer, 2016-2017 Andrew Lelechenko, Carter Schonwald@@ -22,7 +22,7 @@  category            : Math, Algorithms, Number Theory -tested-with         : GHC==7.8.4, GHC==7.10.3, GHC==8.0.2+tested-with         : GHC==7.8.4, GHC==7.10.3, GHC==8.0.2, GHC==8.2.1  extra-source-files  : Changes @@ -50,7 +50,12 @@     exposed-modules     : Math.NumberTheory.ArithmeticFunctions                           Math.NumberTheory.ArithmeticFunctions.Class                           Math.NumberTheory.ArithmeticFunctions.Standard+                          Math.NumberTheory.Curves.Montgomery                           Math.NumberTheory.Moduli+                          Math.NumberTheory.Moduli.Chinese+                          Math.NumberTheory.Moduli.Class+                          Math.NumberTheory.Moduli.Jacobi+                          Math.NumberTheory.Moduli.Sqrt                           Math.NumberTheory.MoebiusInversion                           Math.NumberTheory.MoebiusInversion.Int                           Math.NumberTheory.Recurrencies.Bilinear@@ -75,6 +80,7 @@                           Math.NumberTheory.Primes.Heap                           Math.NumberTheory.UniqueFactorisation                           Math.NumberTheory.Zeta+                          GHC.TypeNats.Compat     other-modules       : Math.NumberTheory.Utils                           Math.NumberTheory.Unsafe                           Math.NumberTheory.Primes.Counting.Impl@@ -123,9 +129,9 @@   ghc-options:          -Wall   main-is:              Test.hs   default-language: Haskell2010-  build-depends:        base >= 4.6 && < 5+  build-depends:        arithmoi+                      , base >= 4.6 && < 5                       , containers >= 0.5 && < 0.6-                      , arithmoi >= 0.5 && < 0.6                       , tasty >= 0.10 && < 0.12                       , tasty-smallcheck >= 0.8 && < 0.9                       , tasty-quickcheck >= 0.9 && < 0.10@@ -136,14 +142,20 @@                       , integer-gmp < 1.1   if impl(ghc < 7.10)     build-depends     : nats >= 1 && <1.2+  if impl(ghc < 8.0)+    build-depends     : semigroups >= 0.8    other-modules :   Math.NumberTheory.ArithmeticFunctionsTests+                  , Math.NumberTheory.CurvesTests                   , Math.NumberTheory.GaussianIntegersTests                   , Math.NumberTheory.GCDTests                   , Math.NumberTheory.GCD.LowLevelTests                   , Math.NumberTheory.Recurrencies.LinearTests                   , Math.NumberTheory.Recurrencies.BilinearTests-                  , Math.NumberTheory.ModuliTests+                  , Math.NumberTheory.Moduli.ChineseTests+                  , Math.NumberTheory.Moduli.ClassTests+                  , Math.NumberTheory.Moduli.JacobiTests+                  , Math.NumberTheory.Moduli.SqrtTests                   , Math.NumberTheory.Powers.CubesTests                   , Math.NumberTheory.MoebiusInversionTests                   , Math.NumberTheory.MoebiusInversion.IntTests

@@ -19,14 +19,37 @@     . mkStdGen     $ salt + bits +-- | bases by Jim Sinclair, https://miller-rabin.appspot.com+fermatBases :: [Integer]+fermatBases = [2, 325, 9375, 28178, 450775, 9780504, 1795265022]++isStrongFermat :: Integer -> Bool+isStrongFermat n = all (isStrongFermatPP n) fermatBases++isFermat :: Integer -> Bool+isFermat n = all (isFermatPP n) fermatBases+ comparePrimalityTests :: Int -> Benchmark comparePrimalityTests bits = bgroup ("primality" ++ show bits)-  [ bench "isPrime"         $ nf (map isPrime)           ns-  , bench "millerRabinV 0"  $ nf (map $ millerRabinV  0) ns-  , bench "millerRabinV 10" $ nf (map $ millerRabinV 10) ns-  , bench "millerRabinV 50" $ nf (map $ millerRabinV 50) ns+  [ bench "isPrime"          $ nf (map isPrime)           ns+  , bench "millerRabinV 0"   $ nf (map $ millerRabinV  0) ns+  , bench "millerRabinV 10"  $ nf (map $ millerRabinV 10) ns+  , bench "millerRabinV 50"  $ nf (map $ millerRabinV 50) ns+  , bench "isStrongFermatPP" $ nf (map isStrongFermat)    ns+  , bench "isFermatPP"       $ nf (map isFermat)          ns   ]   where     ns = take bits [genInteger 0 bits ..] -benchSuite = bgroup "Primes" $ map comparePrimalityTests [50, 100, 200, 500, 1000, 2000]+compareFactorisation :: Int -> Benchmark+compareFactorisation bits =+  bench ("factorise" ++ show bits) $ nf (map factorise) ns+  where+    ns = take (bits `div` 10) [genInteger 0 bits ..]++benchSuite :: Benchmark+benchSuite = bgroup "Primes" $+  map comparePrimalityTests [50, 100, 200, 500, 1000, 2000]+  +++  map compareFactorisation [50, 60, 70, 80, 90, 100]+

@@ -25,6 +25,7 @@ import Data.Foldable #endif +import Data.List (sort) import qualified Data.Set as S import qualified Data.IntSet as IS @@ -49,10 +50,14 @@ divisorsProperty3 :: Natural -> Bool divisorsProperty3 n = all (\d -> n `mod` d == 0) (runFunction divisorsA n) --- | All divisors of n truly divides n.+-- | 'divisorsA' matches 'divisorsSmallA' divisorsProperty4 :: Int -> Bool divisorsProperty4 n = S.toAscList (runFunction divisorsA n) == IS.toAscList (runFunction divisorsSmallA n) +-- | 'divisorsA' matches 'divisorsListA'+divisorsProperty5 :: Int -> Bool+divisorsProperty5 n = S.toAscList (runFunction divisorsA n) == sort (runFunction divisorsListA n)+ -- | tau matches baseline from OEIS. tauOeis :: Assertion tauOeis = oeisAssertion "A000005" tauA@@ -232,6 +237,7 @@     , testSmallAndQuick "sum . divisors = sigma_1" divisorsProperty2     , testSmallAndQuick "matches definition"       divisorsProperty3     , testSmallAndQuick "divisors = divisorsSmall" divisorsProperty4+    , testSmallAndQuick "divisors = divisorsList"  divisorsProperty5     ]   , testGroup "Tau"     [ testCase "OEIS" tauOeis

@@ -0,0 +1,99 @@+-- |+-- Module:      Math.NumberTheory.CurvesTests+-- Copyright:   (c) 2017 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.Curves+--++{-# LANGUAGE CPP                        #-}+{-# LANGUAGE LambdaCase                 #-}+{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.CurvesTests where++import Test.Tasty+import Test.Tasty.QuickCheck as QC hiding (Positive, NonNegative, generate, getNonNegative)++import GHC.TypeNats.Compat++import Math.NumberTheory.Curves.Montgomery+import Math.NumberTheory.TestUtils++#if __GLASGOW_HASKELL__ < 709+import Data.Word+#endif++(==>?) :: Maybe a -> (a -> Property) -> Property+x ==>? f = case x of+  Nothing -> discard+  Just y  -> f y++isValid :: KnownNat n => Point a24 n -> Property+isValid p+  =    counterexample "x is not reduced by modulo"    (x >= 0 && x < n)+  .&&. counterexample "z is not reduced by modulo"    (z >= 0 && z < n)+  where+    n = pointN p+    x = pointX p+    z = pointZ p++isValid' :: KnownNat n => Point a24 n -> Bool+isValid' p+  =  (x >= 0 && x < n)+  && (z >= 0 && z < n)+  where+    n = pointN p+    x = pointX p+    z = pointZ p++newPointRangeProperty :: Shrink2 (Positive Integer) -> Shrink2 (Positive Integer) -> Property+newPointRangeProperty (Shrink2 (Positive s)) (Shrink2 (Positive n)) = newPoint s n ==>? \case+  SomePoint p -> isValid p++multiplyRangeProperty :: Shrink2 (Positive Integer) -> Shrink2 (Positive Integer) -> Shrink2 Word -> Property+multiplyRangeProperty (Shrink2 (Positive s)) (Shrink2 (Positive n)) (Shrink2 k) = newPoint s n ==>? \case+  SomePoint p -> isValid' p ==> isValid (multiply k p)++doubleRangeProperty :: Shrink2 (Positive Integer) -> Shrink2 (Positive Integer) -> Shrink2 Word -> Property+doubleRangeProperty (Shrink2 (Positive s)) (Shrink2 (Positive n)) (Shrink2 k) = newPoint s n ==>? \case+  SomePoint p -> isValid' p ==> isValid' kp ==> isValid (double kp)+    where+      kp = multiply k p++addRangeProperty :: Shrink2 (Positive Integer) -> Shrink2 (Positive Integer) -> Shrink2 Word -> Shrink2 Word -> Property+addRangeProperty (Shrink2 (Positive s)) (Shrink2 (Positive n)) (Shrink2 k) (Shrink2 l) = newPoint s n ==>? \case+  SomePoint p -> isValid' p ==> isValid' kp ==> isValid' lp ==> isValid' klp ==> isValid (add kp lp klp)+    where+      kp  = multiply  k      p+      lp  = multiply      l  p+      klp = multiply (k + l) p++doubleAndMultiplyProperty :: Shrink2 (Positive Integer) -> Shrink2 (Positive Integer) -> Shrink2 Word -> Property+doubleAndMultiplyProperty (Shrink2 (Positive s)) (Shrink2 (Positive n)) (Shrink2 k) = newPoint s n ==>? \case+  SomePoint p+    -> k < maxBound `div` 2 ==> double (multiply k p) === multiply (2 * k) p++addAndMultiplyProperty :: Shrink2 (Positive Integer) -> Shrink2 (Positive Integer) -> Shrink2 Word -> Shrink2 Word -> Property+addAndMultiplyProperty (Shrink2 (Positive s)) (Shrink2 (Positive n)) (Shrink2 k) (Shrink2 l) = newPoint s n ==>? \case+  SomePoint p+    -> k < maxBound `div` 3 && l < maxBound `div` 3 && pointX kp /= 0 && gcd n (pointZ kp) == 1 && gcd n (pointZ lp) == 1 && gcd n (pointZ klp) == 1+    ==> add kp lp klp === k2lp+    where+      kp   = multiply k           p+      lp   = multiply l           p+      klp  = multiply (k + l)     p+      k2lp = multiply (k + 2 * l) p++testSuite :: TestTree+testSuite = localOption (QuickCheckMaxRatio 100) $+  localOption (QuickCheckTests 1000) $ testGroup "Montgomery"+  [ QC.testProperty "range of newPoint"        newPointRangeProperty+  , QC.testProperty "range of double"          doubleRangeProperty+  , QC.testProperty "range of add"             addRangeProperty+  , QC.testProperty "range of multiply"        multiplyRangeProperty+  , QC.testProperty "double matches multiply"  doubleAndMultiplyProperty+  , QC.testProperty "add matches multiply"     addAndMultiplyProperty+  ]

@@ -0,0 +1,48 @@+-- |+-- Module:      Math.NumberTheory.Moduli.ChineseTests+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.Moduli.Chinese+--++{-# LANGUAGE CPP             #-}+{-# LANGUAGE ViewPatterns    #-}++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.Moduli.ChineseTests+  ( testSuite+  ) where++import Test.Tasty++import Control.Arrow+import Data.List (tails)++import Math.NumberTheory.Moduli hiding (invertMod)+import Math.NumberTheory.TestUtils++-- | Check that 'chineseRemainder' is defined iff modulos are coprime.+--   Also check that the result is a solution of input modular equations.+chineseRemainderProperty :: [(Integer, Positive Integer)] -> Bool+chineseRemainderProperty rms' = case chineseRemainder rms of+  Nothing -> not areCoprime+  Just n  -> areCoprime && map (n `mod`) ms == zipWith mod rs ms+  where+    rms = map (second getPositive) rms'+    (rs, ms) = unzip rms+    areCoprime = all (== 1) [ gcd m1 m2 | (m1 : m2s) <- tails ms, m2 <- m2s ]++-- | Check that 'chineseRemainder' matches 'chineseRemainder2'.+chineseRemainder2Property :: Integer -> Positive Integer -> Integer -> Positive Integer -> Bool+chineseRemainder2Property r1 (Positive m1) r2 (Positive m2) = gcd m1 m2 /= 1+  || Just (chineseRemainder2 (r1, m1) (r2, m2)) == chineseRemainder [(r1, m1), (r2, m2)]++testSuite :: TestTree+testSuite = testGroup "Chinese"+  [ testSmallAndQuick "chineseRemainder"  chineseRemainderProperty+  , testSmallAndQuick "chineseRemainder2" chineseRemainder2Property+  ]

@@ -0,0 +1,168 @@+-- |+-- Module:      Math.NumberTheory.Moduli.ClassTests+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.Moduli.Class+--++{-# LANGUAGE CPP             #-}+{-# LANGUAGE ViewPatterns    #-}++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.Moduli.ClassTests+  ( testSuite+  ) where++import Test.Tasty++import Data.Bits+import Data.Maybe+import Numeric.Natural++import Math.NumberTheory.Moduli hiding (invertMod)+import Math.NumberTheory.TestUtils++invertMod :: Integer -> Integer -> Maybe SomeMod+invertMod x m = invertSomeMod (x `modulo` fromInteger m)++powerMod :: Integral a => Integer -> a -> Integer -> SomeMod+powerMod b e m = (b `modulo` fromInteger m) ^ e++-- | Check that 'invertMod' inverts numbers modulo.+invertModProperty :: AnySign Integer -> Positive Integer -> Bool+invertModProperty (AnySign k) (Positive m) = case invertMod k m of+  Nothing            -> k `mod` m == 0 || gcd k m > 1+  Just InfMod{}      -> False+  Just (SomeMod inv) -> gcd k m == 1 && k * getVal inv `mod` m == 1++-- | Check that 'powerMod' is multiplicative by first argument.+powerModProperty2 :: (Integral a, Bits a) => NonNegative a -> AnySign Integer -> AnySign Integer -> Positive Integer -> Bool+powerModProperty2 (NonNegative e) (AnySign b1) (AnySign b2) (Positive m)+  =  e < 0 && (isNothing (invertMod b1 m) || isNothing (invertMod b2 m))+  || pm1 * pm2 == pm12+  where+    pm1  = powerMod b1  e m+    pm2  = powerMod b2  e m+    pm12 = powerMod (b1 * b2) e m++-- | Check that 'powerMod' is additive by second argument.+powerModProperty3 :: (Integral a, Bits a) => NonNegative a -> NonNegative a -> AnySign Integer -> Positive Integer -> Bool+powerModProperty3 (NonNegative e1) (NonNegative 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 == pm12+  where+    pm1  = powerMod b e1 m+    pm2  = powerMod b e2 m+    pm12 = powerMod b (e1 + e2) m++-- | Specialized to trigger 'powerModInteger'.+powerModProperty2_Integer :: NonNegative Integer -> AnySign Integer -> AnySign Integer -> Positive Integer -> Bool+powerModProperty2_Integer = powerModProperty2++-- | Specialized to trigger 'powerModInteger'.+powerModProperty3_Integer :: NonNegative Integer -> NonNegative Integer -> AnySign Integer -> Positive Integer -> Bool+powerModProperty3_Integer = powerModProperty3++someModAddProperty :: Integer -> Positive Natural -> Integer -> Positive Natural -> Bool+someModAddProperty x1 (Positive m1) x2 (Positive m2) = case x1 `modulo` m1 + x2 `modulo` m2 of+  SomeMod z -> getMod z == m3 && getVal z == x3+  InfMod{}  -> False+  where+    m3 = toInteger $ m1 `gcd` m2+    x3 = (x1 + x2) `mod` m3++someModSubProperty :: Integer -> Positive Natural -> Integer -> Positive Natural -> Bool+someModSubProperty x1 (Positive m1) x2 (Positive m2) = case x1 `modulo` m1 - x2 `modulo` m2 of+  SomeMod z -> getMod z == m3 && getVal z == x3+  InfMod{}  -> False+  where+    m3 = toInteger $ m1 `gcd` m2+    x3 = (x1 - x2) `mod` m3++someModMulProperty :: Integer -> Positive Natural -> Integer -> Positive Natural -> Bool+someModMulProperty x1 (Positive m1) x2 (Positive m2) = case (x1 `modulo` m1) * (x2 `modulo` m2) of+  SomeMod z -> getMod z == m3 && getVal z == x3+  InfMod{}  -> False+  where+    m3 = toInteger $ m1 `gcd` m2+    x3 = (x1 * x2) `mod` m3++someModNegProperty :: Integer -> Positive Natural -> Bool+someModNegProperty x1 (Positive m1) = case negate (x1 `modulo` m1) of+  SomeMod z -> getMod z == m3 && getVal z == x3+  InfMod{}  -> False+  where+    m3 = toInteger m1+    x3 = negate x1 `mod` m3++someModAbsSignumProperty :: Integer -> Positive Natural -> Bool+someModAbsSignumProperty x (Positive m) = z == abs z * signum z+  where+    z = x `modulo` m++infModAddProperty :: Integer -> Positive Natural -> Integer -> Bool+infModAddProperty x1 (Positive m1) x2 = case x1 `modulo` m1 + fromInteger x2 of+  SomeMod z -> getMod z == m3 && getVal z == x3+  InfMod{}  -> False+  where+    m3 = toInteger m1+    x3 = (x1 + x2) `mod` m3++infModSubProperty :: Integer -> Positive Natural -> Integer -> Bool+infModSubProperty x1 (Positive m1) x2 = case x1 `modulo` m1 - fromInteger x2 of+  SomeMod z -> getMod z == m3 && getVal z == x3+  InfMod{}  -> False+  where+    m3 = toInteger m1+    x3 = (x1 - x2) `mod` m3++infModMulProperty :: Integer -> Positive Natural -> Integer -> Bool+infModMulProperty x1 (Positive m1) x2 = case x1 `modulo` m1 * fromInteger x2 of+  SomeMod z -> getMod z == m3 && getVal z == x3+  InfMod{}  -> False+  where+    m3 = toInteger m1+    x3 = (x1 * x2) `mod` m3++getValModProperty :: Integer -> Positive Natural -> Bool+getValModProperty x (Positive m) = case z of+  SomeMod t -> z == getVal t `modulo` getNatMod t && z == toInteger (getNatVal t) `modulo` fromInteger (getMod t)+  InfMod{} -> False+  where+    z = x `modulo` m++testSuite :: TestTree+testSuite = testGroup "Class"+  [ testSmallAndQuick "invertMod" invertModProperty+  , testGroup "powerMod"+    [ testGroup "generic"+      [ testIntegralProperty "multiplicative by base"   powerModProperty2+      , testSameIntegralProperty "additive by exponent" powerModProperty3+      ]+    , testGroup "Integer"+      [ testSmallAndQuick "multiplicative by base"  powerModProperty2_Integer+      , testSmallAndQuick "additive by exponent"    powerModProperty3_Integer+      ]+    ]+  , testGroup "SomeMod"+    [ testSmallAndQuick "add" someModAddProperty+    , testSmallAndQuick "sub" someModSubProperty+    , testSmallAndQuick "mul" someModMulProperty+    , testSmallAndQuick "neg" someModNegProperty+    , testSmallAndQuick "abs" someModAbsSignumProperty+    ]+  , testGroup "InfMod"+    [ testSmallAndQuick "add" infModAddProperty+    , testSmallAndQuick "sub" infModSubProperty+    , testSmallAndQuick "mul" infModMulProperty+    ]+  , testSmallAndQuick "getVal/getMod" getValModProperty+  ]

@@ -0,0 +1,66 @@+-- |+-- Module:      Math.NumberTheory.Moduli.JacobiTests+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.Moduli.Jacobi+--++{-# LANGUAGE CPP                        #-}+{-# LANGUAGE ViewPatterns               #-}++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.Moduli.JacobiTests+  ( testSuite+  ) where++import Test.Tasty++import Data.Bits+import Data.Semigroup++import Math.NumberTheory.Moduli hiding (invertMod)+import Math.NumberTheory.TestUtils++-- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 2+jacobiProperty2 :: (Integral a, Bits a) => AnySign a -> (MyCompose Positive Odd) a -> Bool+jacobiProperty2 (AnySign a) (MyCompose (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 -> (MyCompose Positive Odd) a -> Bool+jacobiProperty3 (AnySign a) (MyCompose (Positive (Odd n))) = case jacobi a n of+  MinusOne -> a `gcd` n == 1+  Zero     -> a `gcd` n /= 1+  One      -> a `gcd` n == 1++-- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 4+jacobiProperty4 :: (Integral a, Bits a) => AnySign a -> AnySign a -> (MyCompose Positive Odd) a -> Bool+jacobiProperty4 (AnySign a) (AnySign b) (MyCompose (Positive (Odd n))) = jacobi (a * b) n == jacobi a n <> jacobi b n++jacobiProperty4_Integer :: AnySign Integer -> AnySign Integer -> (MyCompose Positive Odd) Integer -> Bool+jacobiProperty4_Integer = jacobiProperty4++-- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 5+jacobiProperty5 :: (Integral a, Bits a) => AnySign a -> (MyCompose Positive Odd) a -> (MyCompose Positive Odd) a -> Bool+jacobiProperty5 (AnySign a) (MyCompose (Positive (Odd m))) (MyCompose (Positive (Odd n))) = jacobi a (m * n) == jacobi a m <> jacobi a n++jacobiProperty5_Integer :: AnySign Integer -> (MyCompose Positive Odd) Integer -> (MyCompose Positive Odd) Integer -> Bool+jacobiProperty5_Integer = jacobiProperty5++-- https://en.wikipedia.org/wiki/Jacobi_symbol#Properties, item 6+jacobiProperty6 :: (Integral a, Bits a) => (MyCompose Positive Odd) a -> (MyCompose Positive Odd) a -> Bool+jacobiProperty6 (MyCompose (Positive (Odd m))) (MyCompose (Positive (Odd n))) = gcd m n /= 1 || jacobi m n <> jacobi n m == (if m `mod` 4 == 1 || n `mod` 4 == 1 then One else MinusOne)++testSuite :: TestTree+testSuite = testGroup "Jacobi"+  [ testSameIntegralProperty "same modulo n"                jacobiProperty2+  , testSameIntegralProperty "consistent with gcd"          jacobiProperty3+  , testSmallAndQuick        "multiplicative 1"             jacobiProperty4_Integer+  , testSmallAndQuick        "multiplicative 2"             jacobiProperty5_Integer+  , testSameIntegralProperty "law of quadratic reciprocity" jacobiProperty6+  ]

@@ -0,0 +1,121 @@+-- |+-- Module:      Math.NumberTheory.Moduli.SqrtTests+-- Copyright:   (c) 2016 Andrew Lelechenko+-- Licence:     MIT+-- Maintainer:  Andrew Lelechenko <andrew.lelechenko@gmail.com>+-- Stability:   Provisional+--+-- Tests for Math.NumberTheory.Moduli.Sqrt+--++{-# LANGUAGE CPP             #-}+{-# LANGUAGE ViewPatterns    #-}++{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.NumberTheory.Moduli.SqrtTests+  ( testSuite+  ) where++import Test.Tasty+import Test.Tasty.HUnit++import Data.List (nub)++import Math.NumberTheory.Moduli hiding (invertMod)+import Math.NumberTheory.TestUtils++unwrapPP :: (Prime, Power Int) -> (Integer, Int)+unwrapPP (Prime p, Power e) = (p, e)++-- | Check that 'sqrtMod' is defined iff a quadratic residue exists.+--   Also check that the result is a solution of input modular equation.+sqrtModPProperty :: AnySign Integer -> Prime -> Bool+sqrtModPProperty (AnySign n) (Prime p) = case sqrtModP n p of+  Nothing -> jacobi n p == MinusOne+  Just rt -> (p == 2 || jacobi n p /= MinusOne) && rt ^ 2 `mod` p == n `mod` p++sqrtModPListProperty :: AnySign Integer -> Prime -> Bool+sqrtModPListProperty (AnySign n) (Prime p) = all (\rt -> rt ^ 2 `mod` p == n `mod` p) (sqrtModPList n p)++sqrtModP'Property :: Positive Integer -> Prime -> Bool+sqrtModP'Property (Positive n) (Prime p) = (p /= 2 && jacobi n p /= One) || rt ^ 2 `mod` p == n `mod` p+  where+    rt = sqrtModP' n p++tonelliShanksProperty1 :: Positive Integer -> Prime -> Bool+tonelliShanksProperty1 (Positive n) (Prime p) = p `mod` 4 /= 1 || jacobi n p /= One || rt ^ 2 `mod` p == n `mod` p+  where+    rt = tonelliShanks n p++tonelliShanksProperty2 :: Prime -> Bool+tonelliShanksProperty2 (Prime p) = p `mod` 4 /= 1 || rt ^ 2 `mod` p == n `mod` p+  where+    n  = head $ filter (\s -> jacobi s p == One) [2..p-1]+    rt = tonelliShanks n p++tonelliShanksSpecialCases :: Assertion+tonelliShanksSpecialCases =+  assertEqual "OEIS A002224" [6, 32, 219, 439, 1526, 2987, 22193, 11740, 13854, 91168, 326277, 232059, 3230839, 4379725, 11754394, 32020334, 151024619, 345641931, 373671108, 1857111865, 8110112775, 4184367042] rts+  where+    ps = [17, 73, 241, 1009, 2689, 8089, 33049, 53881, 87481, 483289, 515761, 1083289, 3818929, 9257329, 22000801, 48473881, 175244281, 427733329, 898716289, 8114538721, 9176747449, 23616331489]+    rts = map (\p -> tonelliShanks 2 p) ps++sqrtModPPProperty :: AnySign Integer -> (Prime, Power Int) -> Bool+sqrtModPPProperty (AnySign n) (Prime p, Power e) = gcd n p > 1 || case sqrtModPP n (p, e) of+  Nothing -> True+  Just rt -> rt ^ 2 `mod` (p ^ e) == n `mod` (p ^ e)++sqrtModPPBase2Property :: AnySign Integer -> Power Int -> Bool+sqrtModPPBase2Property n e = sqrtModPPProperty n (Prime 2, e)++sqrtModPPSpecialCase1 :: Assertion+sqrtModPPSpecialCase1 =+  assertEqual "sqrtModPP 16 2 2 = 4" (Just 0) (sqrtModPP 16 (2, 2))++sqrtModPPSpecialCase2 :: Assertion+sqrtModPPSpecialCase2 =+  assertEqual "sqrtModPP 16 3 2 = 4" (Just 4) (sqrtModPP 16 (3, 2))++sqrtModPPListProperty :: AnySign Integer -> (Prime, Power Int) -> Bool+sqrtModPPListProperty (AnySign n) (Prime p, Power e) = gcd n p > 1+  || all (\rt -> rt ^ 2 `mod` (p ^ e) == n `mod` (p ^ e)) (sqrtModPPList n (p, e))++sqrtModFProperty :: AnySign Integer -> [(Prime, Power Int)] -> Bool+sqrtModFProperty (AnySign n) (map unwrapPP -> pes) = case sqrtModF n pes of+  Nothing -> True+  Just rt -> all (\(p, e) -> rt ^ 2 `mod` (p ^ e) == n `mod` (p ^ e)) pes++sqrtModFListProperty :: AnySign Integer -> [(Prime, Power Int)] -> Bool+sqrtModFListProperty (AnySign n) (map unwrapPP -> pes)+  = nub ps /= ps || all+    (\rt -> all (\(p, e) -> rt ^ 2 `mod` (p ^ e) == n `mod` (p ^ e)) pes)+    (sqrtModFList n pes)+  where+    ps = map fst pes++sqrtModFListSpecialCase :: Assertion+sqrtModFListSpecialCase =+  assertEqual "sqrtModPPList 0 [(2,1), (3,1), (5,1)]" [0] (sqrtModFList 0 [(2,1), (3,1), (5,1)])++testSuite :: TestTree+testSuite = testGroup "Sqrt"+  [ testSmallAndQuick "sqrtModP"         sqrtModPProperty+  , testSmallAndQuick "sqrtModPList"     sqrtModPListProperty+  , testSmallAndQuick "sqrtModP'"        sqrtModP'Property+  , testGroup "tonelliShanks"+    [ testSmallAndQuick "generic"          tonelliShanksProperty1+    , testSmallAndQuick "smallest residue" tonelliShanksProperty2+    , testCase          "OEIS A002224"     tonelliShanksSpecialCases+    ]+  , testGroup "sqrtModPP"+    [ testSmallAndQuick "generic"        sqrtModPPProperty+    , testSmallAndQuick "_  2 _"         sqrtModPPBase2Property+    , testCase          "16 2 2"         sqrtModPPSpecialCase1+    , testCase          "16 3 2"         sqrtModPPSpecialCase2+    ]+  , testSmallAndQuick "sqrtModPPList"    sqrtModPPListProperty+  , testSmallAndQuick "sqrtModF"         sqrtModFProperty+  , testSmallAndQuick "sqrtModFList"     sqrtModFListProperty+  , testCase          "sqrtModFList 0 [(2,1), (3,1), (5,1)]" sqrtModFListSpecialCase+  ]

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

@@ -28,6 +28,12 @@ totientSumSpecialCase1 :: Assertion totientSumSpecialCase1 = assertEqual "totientSum" 4496 (totientSum 121) +totientSumSpecialCase2 :: Assertion+totientSumSpecialCase2 = assertEqual "totientSum" 0 (totientSum (-9001))++totientSumZero :: Assertion+totientSumZero = assertEqual "totientSum" 0 (totientSum 0)+ generalInversionProperty :: (Int -> Int) -> Positive Int -> Bool generalInversionProperty g (Positive n)   =  g n == sum [f (n `quot` k) | k <- [1 .. n]]@@ -40,6 +46,8 @@   [ testGroup "totientSum"     [ testSmallAndQuick "matches definitions" totientSumProperty     , testCase          "special case 1"      totientSumSpecialCase1+    , testCase          "special case 2"      totientSumSpecialCase2+    , testCase          "zero"                totientSumZero     ]   , QC.testProperty "generalInversion" generalInversionProperty   ]

@@ -28,6 +28,12 @@ totientSumSpecialCase1 :: Assertion totientSumSpecialCase1 = assertEqual "totientSum" 4496 (totientSum 121) +totientSumSpecialCase2 :: Assertion+totientSumSpecialCase2 = assertEqual "totientSum" 0 (totientSum (-9001))++totientSumZero :: Assertion+totientSumZero = assertEqual "totientSum" 0 (totientSum 0)+ generalInversionProperty :: (Int -> Integer) -> Positive Int -> Bool generalInversionProperty g (Positive n)   =  g n == sum [f (n `quot` k) | k <- [1 .. n]]@@ -40,6 +46,8 @@   [ testGroup "totientSum"     [ testSmallAndQuick "matches definitions" totientSumProperty     , testCase          "special case 1"      totientSumSpecialCase1+    , testCase          "special case 2"      totientSumSpecialCase2+    , testCase          "zero"                totientSumZero     ]   , QC.testProperty "generalInversion" generalInversionProperty   ]

@@ -17,10 +17,35 @@ import Test.Tasty import Test.Tasty.HUnit +import Data.List (nub, sort)+ import Math.NumberTheory.Primes.Factorisation import Math.NumberTheory.Primes.Testing import Math.NumberTheory.TestUtils +specialCases :: [(Integer, [(Integer, Int)])]+specialCases =+  [ (4181339589500970917,[(15034813,1),(278110515209,1)])+  , (4181339589500970918,[(2,1),(3,2),(7,1),(2595773,1),(12784336241,1)])+  , (2227144715990344929,[(3,1),(317,1),(17381911,1),(134731889,1)])+  , (10489674846272137811130167281,[(1312601,1),(9555017,1),(836368815445393,1)])+  , (10489674846272137811130167282,[(2,1),(17,1),(577,1),(3863,1),(179347163,1),(771770327021,1)])+  , (10489674846272137811130167283,[(3,1),(7,1),(4634410717,1),(107782489838601619,1)])+  , (10489674846272137811130167287,[(4122913189601,1),(2544238591472087,1)])+  , (6293073306208101456461600748,[(2,2),(3,1),(1613,1),(69973339,1),(4646378436563447,1)])+  , (6293073306208101456461600749,[(7,1),(103,1),(4726591,1),(1846628365511484259,1)])+  , (6293073306208101456461600750,[(2,1),(5,3),(239,1),(34422804769,1),(3059698456333,1)])+  , (6293073306208101456461600751,[(3,1),(13523,1),(1032679,1),(150211485989006401,1)])+  , (6293073306208101456461600753,[(19391,1),(372473053129,1),(871300023127,1)])+  , (6293073306208101456461600754,[(2,1),(3,2),(11,1),(13,1),(71,1),(2311,1),(22859,1),(7798621,1),(83583569,1)])+  , (11999991291828813663324577057,[(14381453,1),(10088205181,1),(82711187849,1)])+  , (11999991291828813663324577062,[(2,1),(3,1),(7,1),(3769,1),(634819511,1),(119413997449529,1)])+  , (16757651897802863152387219654541878160,[(2,4),(5,1),(12323,1),(1424513,1),(6205871923,1),(1922815011093901,1)])+  , (16757651897802863152387219654541878162,[(2,1),(29,1),(78173,1),(401529283,1),(1995634649,1),(4612433663779,1)])+  , (16757651897802863152387219654541878163,[(11,1),(31,1),(112160981904206269,1),(438144115295608147,1)])+  , (16757651897802863152387219654541878166,[(2,1),(23,1),(277,1),(505353699591289,1),(2602436338718275457,1)])+  ]+ factoriseProperty1 :: Assertion factoriseProperty1 = assertEqual "0" [] (factorise 1) @@ -31,14 +56,24 @@ factoriseProperty3 (Positive n) = all (isPrime . fst) (factorise n)  factoriseProperty4 :: Positive Integer -> Bool-factoriseProperty4 (Positive n) = product (map (uncurry (^)) (factorise n)) == n+factoriseProperty4 (Positive n) = bases == nub (sort bases)+  where+    bases = map fst $ factorise n +factoriseProperty5 :: Positive Integer -> Bool+factoriseProperty5 (Positive n) = product (map (uncurry (^)) (factorise n)) == n++factoriseProperty6 :: (Integer, [(Integer, Int)]) -> Assertion+factoriseProperty6 (n, fs) = assertEqual (show n) fs (factorise n)+ testSuite :: TestTree testSuite = testGroup "Factorisation"-  [ testGroup "factorise"-    [ testCase          "0"                factoriseProperty1-    , testSmallAndQuick "negate"                  factoriseProperty2-    , testSmallAndQuick          "bases are prime" factoriseProperty3-    , testSmallAndQuick          "factorback" factoriseProperty4-    ]+  [ testGroup "factorise" $+    [ testCase          "0"                              factoriseProperty1+    , testSmallAndQuick "negate"                         factoriseProperty2+    , testSmallAndQuick "bases are prime"                factoriseProperty3+    , testSmallAndQuick "bases are ordered and distinct" factoriseProperty4+    , testSmallAndQuick "factorback"                     factoriseProperty5+    ] +++    map (\x -> testCase ("special case " ++ show (fst x)) (factoriseProperty6 x)) specialCases   ]

@@ -35,7 +35,7 @@ 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)+import Math.NumberTheory.Primes (isPrime, nthPrime)  ------------------------------------------------------------------------------- -- AnySign@@ -171,10 +171,12 @@   deriving (Eq, Ord, Show)  instance Arbitrary Prime where-  arbitrary = Prime <$> arbitrary `suchThat` (\p -> p > 0 && isPrime p)+  arbitrary = do+    n <- arbitrary+    return $ Prime $ head $ filter isPrime [abs n ..]  instance Monad m => Serial m Prime where-  series = Prime <$> series `suchThatSerial` (\p -> p > 0 && isPrime p)+  series = Prime . nthPrime <$> series `suchThatSerial` (> 0)  ------------------------------------------------------------------------------- -- Utils

@@ -6,7 +6,10 @@ import qualified Math.NumberTheory.Recurrencies.BilinearTests as RecurrenciesBilinear import qualified Math.NumberTheory.Recurrencies.LinearTests as RecurrenciesLinear -import qualified Math.NumberTheory.ModuliTests as Moduli+import qualified Math.NumberTheory.Moduli.ChineseTests as ModuliChinese+import qualified Math.NumberTheory.Moduli.ClassTests as ModuliClass+import qualified Math.NumberTheory.Moduli.JacobiTests as ModuliJacobi+import qualified Math.NumberTheory.Moduli.SqrtTests as ModuliSqrt  import qualified Math.NumberTheory.MoebiusInversionTests as MoebiusInversion import qualified Math.NumberTheory.MoebiusInversion.IntTests as MoebiusInversionInt@@ -28,6 +31,7 @@ import qualified Math.NumberTheory.ArithmeticFunctionsTests as ArithmeticFunctions import qualified Math.NumberTheory.UniqueFactorisationTests as UniqueFactorisation import qualified Math.NumberTheory.ZetaTests as Zeta+import qualified Math.NumberTheory.CurvesTests as Curves  main :: IO () main = defaultMain tests@@ -49,7 +53,10 @@     , RecurrenciesBilinear.testSuite     ]   , testGroup "Moduli"-    [ Moduli.testSuite+    [ ModuliChinese.testSuite+    , ModuliClass.testSuite+    , ModuliJacobi.testSuite+    , ModuliSqrt.testSuite     ]   , testGroup "MoebiusInversion"     [ MoebiusInversion.testSuite@@ -74,5 +81,8 @@     ]   , testGroup "Zeta"     [ Zeta.testSuite+    ]+  , testGroup "Curves"+    [ Curves.testSuite     ]   ]