Skip to content

Commit

Permalink
Merge pull request #180 from b-mehta/dirichlet-characters
Browse files Browse the repository at this point in the history 
Dirichlet characters
  • Loading branch information
@Bodigrim
Bodigrim authored Jan 13, 2020
2 parents 0592774 + d71843c commit 54a4448
Show file tree
Hide file tree
Showing 8 changed files with 1,015 additions and 104 deletions.
614 changes: 614 additions & 0 deletions Math/NumberTheory/DirichletCharacters.hs
View file

Large diffs are not rendered by default.

126 changes: 126 additions & 0 deletions Math/NumberTheory/Moduli/Internal.hs
View file
Original file line number Diff line number Diff line change
@@ -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
13 changes: 13 additions & 0 deletions Math/NumberTheory/Moduli/JacobiSymbol.hs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,6 +17,7 @@
module Math.NumberTheory.Moduli.JacobiSymbol
( JacobiSymbol(..)
, jacobi
, symbolToNum
) where

import Data.Bits
Expand Down Expand Up @@ -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.
Expand Down
106 changes: 2 additions & 104 deletions Math/NumberTheory/Moduli/Multiplicative.hs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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.
Expand Down Expand Up @@ -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>.
--
Expand Down Expand Up @@ -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
1 change: 1 addition & 0 deletions Math/NumberTheory/Moduli/Sqrt.hs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@ module Math.NumberTheory.Moduli.Sqrt
-- * Jacobi symbol
, JacobiSymbol(..)
, jacobi
, symbolToNum
) where

import Control.Monad (liftM2)
Expand Down
3 changes: 3 additions & 0 deletions arithmoi.cabal
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -151,6 +153,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
Expand Down
Loading

0 comments on commit 54a4448

Please sign in to comment.