diff --git a/Changes b/Changes --- a/Changes +++ b/Changes @@ -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. diff --git a/GHC/TypeNats/Compat.hs b/GHC/TypeNats/Compat.hs new file mode 100644 --- /dev/null +++ b/GHC/TypeNats/Compat.hs @@ -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 diff --git a/Math/NumberTheory/ArithmeticFunctions/Standard.hs b/Math/NumberTheory/ArithmeticFunctions/Standard.hs --- a/Math/NumberTheory/ArithmeticFunctions/Standard.hs +++ b/Math/NumberTheory/ArithmeticFunctions/Standard.hs @@ -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 } diff --git a/Math/NumberTheory/Curves/Montgomery.hs b/Math/NumberTheory/Curves/Montgomery.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Curves/Montgomery.hs @@ -0,0 +1,172 @@ +-- | +-- Module: Math.NumberTheory.Curves.Montgomery +-- Copyright: (c) 2017 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- 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 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 +-- 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 +-- 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 +-- 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) diff --git a/Math/NumberTheory/Moduli.hs b/Math/NumberTheory/Moduli.hs --- a/Math/NumberTheory/Moduli.hs +++ b/Math/NumberTheory/Moduli.hs @@ -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 ±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 diff --git a/Math/NumberTheory/Moduli/Chinese.hs b/Math/NumberTheory/Moduli/Chinese.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Moduli/Chinese.hs @@ -0,0 +1,59 @@ +-- | +-- Module: Math.NumberTheory.Moduli.Chinese +-- Copyright: (c) 2011 Daniel Fischer +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- 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 diff --git a/Math/NumberTheory/Moduli/Class.hs b/Math/NumberTheory/Moduli/Class.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Moduli/Class.hs @@ -0,0 +1,306 @@ +-- | +-- Module: Math.NumberTheory.Moduli.Class +-- Copyright: (c) 2017 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- 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 #-} diff --git a/Math/NumberTheory/Moduli/Jacobi.hs b/Math/NumberTheory/Moduli/Jacobi.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Moduli/Jacobi.hs @@ -0,0 +1,150 @@ +-- | +-- Module: Math.NumberTheory.Moduli.Jacobi +-- Copyright: (c) 2011 Daniel Fischer, 2017 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- 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 ±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)] diff --git a/Math/NumberTheory/Moduli/Sqrt.hs b/Math/NumberTheory/Moduli/Sqrt.hs new file mode 100644 --- /dev/null +++ b/Math/NumberTheory/Moduli/Sqrt.hs @@ -0,0 +1,229 @@ +-- | +-- Module: Math.NumberTheory.Moduli.Sqrt +-- Copyright: (c) 2011 Daniel Fischer +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- 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 diff --git a/Math/NumberTheory/MoebiusInversion.hs b/Math/NumberTheory/MoebiusInversion.hs --- a/Math/NumberTheory/MoebiusInversion.hs +++ b/Math/NumberTheory/MoebiusInversion.hs @@ -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 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@. diff --git a/Math/NumberTheory/MoebiusInversion/Int.hs b/Math/NumberTheory/MoebiusInversion/Int.hs --- a/Math/NumberTheory/MoebiusInversion/Int.hs +++ b/Math/NumberTheory/MoebiusInversion/Int.hs @@ -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 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@. diff --git a/Math/NumberTheory/Powers.hs b/Math/NumberTheory/Powers.hs --- a/Math/NumberTheory/Powers.hs +++ b/Math/NumberTheory/Powers.hs @@ -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 diff --git a/Math/NumberTheory/Primes/Factorisation/Certified.hs b/Math/NumberTheory/Primes/Factorisation/Certified.hs --- a/Math/NumberTheory/Primes/Factorisation/Certified.hs +++ b/Math/NumberTheory/Primes/Factorisation/Certified.hs @@ -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 diff --git a/Math/NumberTheory/Primes/Factorisation/Montgomery.hs b/Math/NumberTheory/Primes/Factorisation/Montgomery.hs --- a/Math/NumberTheory/Primes/Factorisation/Montgomery.hs +++ b/Math/NumberTheory/Primes/Factorisation/Montgomery.hs @@ -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 +-- 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) diff --git a/Math/NumberTheory/Primes/Testing/Certificates/Internal.hs b/Math/NumberTheory/Primes/Testing/Certificates/Internal.hs --- a/Math/NumberTheory/Primes/Testing/Certificates/Internal.hs +++ b/Math/NumberTheory/Primes/Testing/Certificates/Internal.hs @@ -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)) diff --git a/Math/NumberTheory/Primes/Testing/Probabilistic.hs b/Math/NumberTheory/Primes/Testing/Probabilistic.hs --- a/Math/NumberTheory/Primes/Testing/Probabilistic.hs +++ b/Math/NumberTheory/Primes/Testing/Probabilistic.hs @@ -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 diff --git a/Math/NumberTheory/Recurrencies/Linear.hs b/Math/NumberTheory/Recurrencies/Linear.hs --- a/Math/NumberTheory/Recurrencies/Linear.hs +++ b/Math/NumberTheory/Recurrencies/Linear.hs @@ -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 diff --git a/arithmoi.cabal b/arithmoi.cabal --- a/arithmoi.cabal +++ b/arithmoi.cabal @@ -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 diff --git a/benchmark/Math/NumberTheory/PrimesBench.hs b/benchmark/Math/NumberTheory/PrimesBench.hs --- a/benchmark/Math/NumberTheory/PrimesBench.hs +++ b/benchmark/Math/NumberTheory/PrimesBench.hs @@ -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] + diff --git a/test-suite/Math/NumberTheory/ArithmeticFunctionsTests.hs b/test-suite/Math/NumberTheory/ArithmeticFunctionsTests.hs --- a/test-suite/Math/NumberTheory/ArithmeticFunctionsTests.hs +++ b/test-suite/Math/NumberTheory/ArithmeticFunctionsTests.hs @@ -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 diff --git a/test-suite/Math/NumberTheory/CurvesTests.hs b/test-suite/Math/NumberTheory/CurvesTests.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/NumberTheory/CurvesTests.hs @@ -0,0 +1,99 @@ +-- | +-- Module: Math.NumberTheory.CurvesTests +-- Copyright: (c) 2017 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- 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 + ] diff --git a/test-suite/Math/NumberTheory/Moduli/ChineseTests.hs b/test-suite/Math/NumberTheory/Moduli/ChineseTests.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/NumberTheory/Moduli/ChineseTests.hs @@ -0,0 +1,48 @@ +-- | +-- Module: Math.NumberTheory.Moduli.ChineseTests +-- Copyright: (c) 2016 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- 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 + ] diff --git a/test-suite/Math/NumberTheory/Moduli/ClassTests.hs b/test-suite/Math/NumberTheory/Moduli/ClassTests.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/NumberTheory/Moduli/ClassTests.hs @@ -0,0 +1,168 @@ +-- | +-- Module: Math.NumberTheory.Moduli.ClassTests +-- Copyright: (c) 2016 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- 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 + ] diff --git a/test-suite/Math/NumberTheory/Moduli/JacobiTests.hs b/test-suite/Math/NumberTheory/Moduli/JacobiTests.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/NumberTheory/Moduli/JacobiTests.hs @@ -0,0 +1,66 @@ +-- | +-- Module: Math.NumberTheory.Moduli.JacobiTests +-- Copyright: (c) 2016 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- 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 + ] diff --git a/test-suite/Math/NumberTheory/Moduli/SqrtTests.hs b/test-suite/Math/NumberTheory/Moduli/SqrtTests.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/NumberTheory/Moduli/SqrtTests.hs @@ -0,0 +1,121 @@ +-- | +-- Module: Math.NumberTheory.Moduli.SqrtTests +-- Copyright: (c) 2016 Andrew Lelechenko +-- Licence: MIT +-- Maintainer: Andrew Lelechenko +-- 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 + ] diff --git a/test-suite/Math/NumberTheory/ModuliTests.hs b/test-suite/Math/NumberTheory/ModuliTests.hs deleted file mode 100644 --- a/test-suite/Math/NumberTheory/ModuliTests.hs +++ /dev/null @@ -1,225 +0,0 @@ --- | --- Module: Math.NumberTheory.ModuliTests --- Copyright: (c) 2016 Andrew Lelechenko --- Licence: MIT --- Maintainer: Andrew Lelechenko --- 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 - ] - ] diff --git a/test-suite/Math/NumberTheory/MoebiusInversion/IntTests.hs b/test-suite/Math/NumberTheory/MoebiusInversion/IntTests.hs --- a/test-suite/Math/NumberTheory/MoebiusInversion/IntTests.hs +++ b/test-suite/Math/NumberTheory/MoebiusInversion/IntTests.hs @@ -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 ] diff --git a/test-suite/Math/NumberTheory/MoebiusInversionTests.hs b/test-suite/Math/NumberTheory/MoebiusInversionTests.hs --- a/test-suite/Math/NumberTheory/MoebiusInversionTests.hs +++ b/test-suite/Math/NumberTheory/MoebiusInversionTests.hs @@ -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 ] diff --git a/test-suite/Math/NumberTheory/Primes/FactorisationTests.hs b/test-suite/Math/NumberTheory/Primes/FactorisationTests.hs --- a/test-suite/Math/NumberTheory/Primes/FactorisationTests.hs +++ b/test-suite/Math/NumberTheory/Primes/FactorisationTests.hs @@ -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 ] diff --git a/test-suite/Math/NumberTheory/TestUtils/Wrappers.hs b/test-suite/Math/NumberTheory/TestUtils/Wrappers.hs --- a/test-suite/Math/NumberTheory/TestUtils/Wrappers.hs +++ b/test-suite/Math/NumberTheory/TestUtils/Wrappers.hs @@ -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 diff --git a/test-suite/Test.hs b/test-suite/Test.hs --- a/test-suite/Test.hs +++ b/test-suite/Test.hs @@ -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 ] ]