Dirichlet characters

Math/NumberTheory/Moduli/Internal.hs

@@ -0,0 +1,126 @@
+-- |
+-- Module:      Math.NumberTheory.Moduli.Internal
+-- Copyright:   (c) 2020 Bhavik Mehta
+-- Licence:     MIT
+-- Maintainer:  Bhavik Mehta <[email protected]>
+--
+-- Multiplicative groups of integers modulo m.
+--
+
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE ScopedTypeVariables #-}
+
+module Math.NumberTheory.Moduli.Internal
+  ( isPrimitiveRoot'
+  , discreteLogarithmPP
+  ) where
+
+import qualified Data.Map as M
+import Data.Maybe
+import Data.Mod
+import Data.Proxy
+import GHC.Integer.GMP.Internals
+import GHC.TypeNats.Compat
+import Numeric.Natural
+
+import Math.NumberTheory.ArithmeticFunctions
+import Math.NumberTheory.Moduli.Chinese
+import Math.NumberTheory.Moduli.Equations
+import Math.NumberTheory.Moduli.Singleton
+import Math.NumberTheory.Primes
+import Math.NumberTheory.Powers.Modular
+import Math.NumberTheory.Roots
+
+-- https://en.wikipedia.org/wiki/Primitive_root_modulo_n#Finding_primitive_roots
+isPrimitiveRoot'
+  :: (Integral a, UniqueFactorisation a)
+  => CyclicGroup a m
+  -> a
+  -> Bool
+isPrimitiveRoot' cg r =
+  case cg of
+    CG2                       -> r == 1
+    CG4                       -> r == 3
+    CGOddPrimePower p k       -> oddPrimePowerTest (unPrime p) k r
+    CGDoubleOddPrimePower p k -> doubleOddPrimePowerTest (unPrime p) k r
+  where
+    oddPrimeTest p g              = let phi  = totient p
+                                        pows = map (\pk -> phi `quot` unPrime (fst pk)) (factorise phi)
+                                        exps = map (\x -> powMod g x p) pows
+                                     in g /= 0 && gcd g p == 1 && notElem 1 exps
+    oddPrimePowerTest p 1 g       = oddPrimeTest p (g `mod` p)
+    oddPrimePowerTest p _ g       = oddPrimeTest p (g `mod` p) && powMod g (p-1) (p*p) /= 1
+    doubleOddPrimePowerTest p k g = odd g && oddPrimePowerTest p k g
+
+-- Implementation of Bach reduction (https://www2.eecs.berkeley.edu/Pubs/TechRpts/1984/CSD-84-186.pdf)
+{-# INLINE discreteLogarithmPP #-}
+discreteLogarithmPP :: Integer -> Word -> Integer -> Integer -> Natural
+discreteLogarithmPP p 1 a b = discreteLogarithmPrime p a b
+discreteLogarithmPP p k a b = fromInteger $ if result < 0 then result + pkMinusPk1 else result
+  where
+    baseSol    = toInteger $ discreteLogarithmPrime p (a `rem` p) (b `rem` p)
+    thetaA     = theta p pkMinusOne a
+    thetaB     = theta p pkMinusOne b
+    pkMinusOne = p^(k-1)
+    pkMinusPk1 = pkMinusOne * (p - 1)
+    c          = (recipModInteger thetaA pkMinusOne * thetaB) `rem` pkMinusOne
+    result     = fromJust $ chineseCoprime (baseSol, p-1) (c, pkMinusOne)
+
+-- compute the homomorphism theta given in https://math.stackexchange.com/a/1864495/418148
+{-# INLINE theta #-}
+theta :: Integer -> Integer -> Integer -> Integer
+theta p pkMinusOne a = (numerator `quot` pk) `rem` pkMinusOne
+  where
+    pk           = pkMinusOne * p
+    p2kMinusOne  = pkMinusOne * pk
+    numerator    = (powModInteger a (pk - pkMinusOne) p2kMinusOne - 1) `rem` p2kMinusOne
+
+-- TODO: Use Pollig-Hellman to reduce the problem further into groups of prime order.
+-- While Bach reduction simplifies the problem into groups of the form (Z/pZ)*, these
+-- have non-prime order, and the Pollig-Hellman algorithm can reduce the problem into
+-- smaller groups of prime order.
+-- In addition, the gcd check before solveLinear is applied in Pollard below will be
+-- made redundant, since n would be prime.
+discreteLogarithmPrime :: Integer -> Integer -> Integer -> Natural
+discreteLogarithmPrime p a b
+  | p < 100000000 = fromIntegral $ discreteLogarithmPrimeBSGS (fromInteger p) (fromInteger a) (fromInteger b)
+  | otherwise     = discreteLogarithmPrimePollard p a b
+
+discreteLogarithmPrimeBSGS :: Int -> Int -> Int -> Int
+discreteLogarithmPrimeBSGS p a b = head [i*m + j | (v,i) <- zip giants [0..m-1], j <- maybeToList (M.lookup v table)]
+  where
+    m        = integerSquareRoot (p - 2) + 1 -- simple way of ceiling (sqrt (p-1))
+    babies   = iterate (.* a) 1
+    table    = M.fromList (zip babies [0..m-1])
+    aInv     = recipModInteger (toInteger a) (toInteger p)
+    bigGiant = fromInteger $ powModInteger aInv (toInteger m) (toInteger p)
+    giants   = iterate (.* bigGiant) b
+    x .* y   = x * y `rem` p
+
+-- TODO: Use more advanced walks, in order to reduce divisions, cf
+-- https://maths-people.anu.edu.au/~brent/pd/rpb231.pdf
+-- This will slightly improve the expected time to collision, and can reduce the
+-- number of divisions performed.
+discreteLogarithmPrimePollard :: Integer -> Integer -> Integer -> Natural
+discreteLogarithmPrimePollard p a b =
+  case concatMap runPollard [(x,y) | x <- [0..n], y <- [0..n]] of
+    (t:_)  -> fromInteger t
+    []     -> error ("discreteLogarithm: pollard's rho failed, please report this as a bug. inputs " ++ show [p,a,b])
+  where
+    n                 = p-1 -- order of the cyclic group
+    halfN             = n `quot` 2
+    mul2 m            = if m < halfN then m * 2 else m * 2 - n
+    sqrtN             = integerSquareRoot n
+    step (xi,!ai,!bi) = case xi `rem` 3 of
+                          0 -> (xi*xi `rem` p, mul2 ai, mul2 bi)
+                          1 -> ( a*xi `rem` p,    ai+1,      bi)
+                          _ -> ( b*xi `rem` p,      ai,    bi+1)
+    initialise (x,y)  = (powModInteger a x n * powModInteger b y n `rem` n, x, y)
+    begin t           = go (step t) (step (step t))
+    check t           = powModInteger a t p == b
+    go tort@(xi,ai,bi) hare@(x2i,a2i,b2i)
+      | xi == x2i, gcd (bi - b2i) n < sqrtN = case someNatVal (fromInteger n) of
+        SomeNat (Proxy :: Proxy n) -> map (toInteger . unMod) $ solveLinear (fromInteger (bi - b2i) :: Mod n) (fromInteger (ai - a2i))
+      | xi == x2i                           = []
+      | otherwise                           = go (step tort) (step (step hare))
+    runPollard        = filter check . begin . initialise

Math/NumberTheory/Moduli/JacobiSymbol.hs

@@ -17,6 +17,7 @@
 module Math.NumberTheory.Moduli.JacobiSymbol
   ( JacobiSymbol(..)
   , jacobi
+  , symbolToNum
   ) where
 
 import Data.Bits
@@ -48,6 +49,18 @@ negJS = \case
   Zero     -> Zero
   One      -> MinusOne
 
+{-# SPECIALISE symbolToNum :: JacobiSymbol -> Integer,
+                              JacobiSymbol -> Int,
+                              JacobiSymbol -> Word,
+                              JacobiSymbol -> Natural
+  #-}
+-- | Convenience function to convert out of a Jacobi symbol
+symbolToNum :: Num a => JacobiSymbol -> a
+symbolToNum = \case
+  Zero -> 0
+  One -> 1
+  MinusOne -> -1
+
 -- | <https://en.wikipedia.org/wiki/Jacobi_symbol Jacobi symbol> of two arguments.
 -- The lower argument (\"denominator\") must be odd and positive,
 -- this condition is checked.

Math/NumberTheory/Moduli/Multiplicative.hs

@@ -7,9 +7,9 @@
 -- Multiplicative groups of integers modulo m.
 --
 
-{-# LANGUAGE BangPatterns #-}
 {-# LANGUAGE ScopedTypeVariables #-}
 {-# LANGUAGE ViewPatterns #-}
+{-# LANGUAGE PatternSynonyms #-}
 
 module Math.NumberTheory.Moduli.Multiplicative
   ( -- * Multiplicative group
@@ -26,22 +26,14 @@ module Math.NumberTheory.Moduli.Multiplicative
 
 import Control.Monad
 import Data.Constraint
-import qualified Data.Map as M
-import Data.Maybe
 import Data.Mod
-import Data.Proxy
 import Data.Semigroup
-import GHC.Integer.GMP.Internals
 import GHC.TypeNats.Compat
 import Numeric.Natural
 
-import Math.NumberTheory.ArithmeticFunctions
-import Math.NumberTheory.Moduli.Chinese
-import Math.NumberTheory.Moduli.Equations
+import Math.NumberTheory.Moduli.Internal
 import Math.NumberTheory.Moduli.Singleton
 import Math.NumberTheory.Primes
-import Math.NumberTheory.Powers.Modular
-import Math.NumberTheory.Roots
 
 -- | This type represents elements of the multiplicative group mod m, i.e.
 -- those elements which are coprime to m. Use @toMultElement@ to construct.
@@ -84,27 +76,6 @@ newtype PrimitiveRoot m = PrimitiveRoot
   }
   deriving (Eq, Show)
 
--- https://en.wikipedia.org/wiki/Primitive_root_modulo_n#Finding_primitive_roots
-isPrimitiveRoot'
-  :: (Integral a, UniqueFactorisation a)
-  => CyclicGroup a m
-  -> a
-  -> Bool
-isPrimitiveRoot' cg r =
-  case cg of
-    CG2                       -> r == 1
-    CG4                       -> r == 3
-    CGOddPrimePower p k       -> oddPrimePowerTest (unPrime p) k r
-    CGDoubleOddPrimePower p k -> doubleOddPrimePowerTest (unPrime p) k r
-  where
-    oddPrimeTest p g              = let phi  = totient p
-                                        pows = map (\pk -> phi `quot` unPrime (fst pk)) (factorise phi)
-                                        exps = map (\x -> powMod g x p) pows
-                                     in g /= 0 && gcd g p == 1 && all (/= 1) exps
-    oddPrimePowerTest p 1 g       = oddPrimeTest p (g `mod` p)
-    oddPrimePowerTest p _ g       = oddPrimeTest p (g `mod` p) && powMod g (p-1) (p*p) /= 1
-    doubleOddPrimePowerTest p k g = odd g && oddPrimePowerTest p k g
-
 -- | Check whether a given modular residue is
 -- a <https://en.wikipedia.org/wiki/Primitive_root_modulo_n primitive root>.
 --
@@ -148,76 +119,3 @@ discreteLogarithm cg (multElement . unPrimitiveRoot -> a) (multElement -> b) = c
   CGDoubleOddPrimePower (unPrime -> p) k
     -> discreteLogarithmPP p k (toInteger (unMod a) `rem` p^k) (toInteger (unMod b) `rem` p^k)
     -- we have the isomorphism t -> t `rem` p^k from (Z/2p^kZ)* -> (Z/p^kZ)*
-
--- Implementation of Bach reduction (https://www2.eecs.berkeley.edu/Pubs/TechRpts/1984/CSD-84-186.pdf)
-{-# INLINE discreteLogarithmPP #-}
-discreteLogarithmPP :: Integer -> Word -> Integer -> Integer -> Natural
-discreteLogarithmPP p 1 a b = discreteLogarithmPrime p a b
-discreteLogarithmPP p k a b = fromInteger $ if result < 0 then result + pkMinusPk1 else result
-  where
-    baseSol    = toInteger $ discreteLogarithmPrime p (a `rem` p) (b `rem` p)
-    thetaA     = theta p pkMinusOne a
-    thetaB     = theta p pkMinusOne b
-    pkMinusOne = p^(k-1)
-    pkMinusPk1 = pkMinusOne * (p - 1)
-    c          = (recipModInteger thetaA pkMinusOne * thetaB) `rem` pkMinusOne
-    result     = fromJust $ chineseCoprime (baseSol, p-1) (c, pkMinusOne)
-
--- compute the homomorphism theta given in https://math.stackexchange.com/a/1864495/418148
-{-# INLINE theta #-}
-theta :: Integer -> Integer -> Integer -> Integer
-theta p pkMinusOne a = (numerator `quot` pk) `rem` pkMinusOne
-  where
-    pk           = pkMinusOne * p
-    p2kMinusOne  = pkMinusOne * pk
-    numerator    = (powModInteger a (pk - pkMinusOne) p2kMinusOne - 1) `rem` p2kMinusOne
-
--- TODO: Use Pollig-Hellman to reduce the problem further into groups of prime order.
--- While Bach reduction simplifies the problem into groups of the form (Z/pZ)*, these
--- have non-prime order, and the Pollig-Hellman algorithm can reduce the problem into
--- smaller groups of prime order.
--- In addition, the gcd check before solveLinear is applied in Pollard below will be
--- made redundant, since n would be prime.
-discreteLogarithmPrime :: Integer -> Integer -> Integer -> Natural
-discreteLogarithmPrime p a b
-  | p < 100000000 = fromIntegral $ discreteLogarithmPrimeBSGS (fromInteger p) (fromInteger a) (fromInteger b)
-  | otherwise     = discreteLogarithmPrimePollard p a b
-
-discreteLogarithmPrimeBSGS :: Int -> Int -> Int -> Int
-discreteLogarithmPrimeBSGS p a b = head [i*m + j | (v,i) <- zip giants [0..m-1], j <- maybeToList (M.lookup v table)]
-  where
-    m        = integerSquareRoot (p - 2) + 1 -- simple way of ceiling (sqrt (p-1))
-    babies   = iterate (.* a) 1
-    table    = M.fromList (zip babies [0..m-1])
-    aInv     = recipModInteger (toInteger a) (toInteger p)
-    bigGiant = fromInteger $ powModInteger aInv (toInteger m) (toInteger p)
-    giants   = iterate (.* bigGiant) b
-    x .* y   = x * y `rem` p
-
--- TODO: Use more advanced walks, in order to reduce divisions, cf
--- https://maths-people.anu.edu.au/~brent/pd/rpb231.pdf
--- This will slightly improve the expected time to collision, and can reduce the
--- number of divisions performed.
-discreteLogarithmPrimePollard :: Integer -> Integer -> Integer -> Natural
-discreteLogarithmPrimePollard p a b =
-  case concatMap runPollard [(x,y) | x <- [0..n], y <- [0..n]] of
-    (t:_)  -> fromInteger t
-    []     -> error ("discreteLogarithm: pollard's rho failed, please report this as a bug. inputs " ++ show [p,a,b])
-  where
-    n                 = p-1 -- order of the cyclic group
-    halfN             = n `quot` 2
-    mul2 m            = if m < halfN then m * 2 else m * 2 - n
-    sqrtN             = integerSquareRoot n
-    step (xi,!ai,!bi) = case xi `rem` 3 of
-                          0 -> (xi*xi `rem` p, mul2 ai, mul2 bi)
-                          1 -> ( a*xi `rem` p,    ai+1,      bi)
-                          _ -> ( b*xi `rem` p,      ai,    bi+1)
-    initialise (x,y)  = (powModInteger a x n * powModInteger b y n `rem` n, x, y)
-    begin t           = go (step t) (step (step t))
-    check t           = powModInteger a t p == b
-    go tort@(xi,ai,bi) hare@(x2i,a2i,b2i)
-      | xi == x2i, gcd (bi - b2i) n < sqrtN = case someNatVal (fromInteger n) of
-        SomeNat (Proxy :: Proxy n) -> map (toInteger . unMod) $ solveLinear (fromInteger (bi - b2i) :: Mod n) (fromInteger (ai - a2i))
-      | xi == x2i                           = []
-      | otherwise                           = go (step tort) (step (step hare))
-    runPollard        = filter check . begin . initialise

Math/NumberTheory/Moduli/Sqrt.hs

@@ -21,6 +21,7 @@ module Math.NumberTheory.Moduli.Sqrt
     -- * Jacobi symbol
   , JacobiSymbol(..)
   , jacobi
+  , symbolToNum
   ) where
 
 import Control.Monad (liftM2)

arithmoi.cabal

@@ -59,6 +59,7 @@ library
     Math.NumberTheory.ArithmeticFunctions.Moebius
     Math.NumberTheory.ArithmeticFunctions.SieveBlock
     Math.NumberTheory.Curves.Montgomery
+    Math.NumberTheory.DirichletCharacters
     Math.NumberTheory.Euclidean
     Math.NumberTheory.Euclidean.Coprimes
     Math.NumberTheory.Moduli
@@ -94,6 +95,7 @@ library
   other-modules:
     Math.NumberTheory.ArithmeticFunctions.Class
     Math.NumberTheory.ArithmeticFunctions.Standard
+    Math.NumberTheory.Moduli.Internal
     Math.NumberTheory.Moduli.JacobiSymbol
     Math.NumberTheory.Moduli.SomeMod
     Math.NumberTheory.Primes.Counting.Approximate
@@ -150,6 +152,7 @@ test-suite spec
     Math.NumberTheory.ArithmeticFunctions.MertensTests
     Math.NumberTheory.ArithmeticFunctions.SieveBlockTests
     Math.NumberTheory.CurvesTests
+    Math.NumberTheory.DirichletCharactersTests
     Math.NumberTheory.EisensteinIntegersTests
     Math.NumberTheory.GaussianIntegersTests
     Math.NumberTheory.EuclideanTests