Discrete log

Math/NumberTheory/GaussianIntegers.hs

@@ -40,7 +40,7 @@ import Data.Maybe (fromMaybe)
 import Data.Ord (comparing)
 import GHC.Generics
 
-import qualified Math.NumberTheory.Moduli as Moduli
+import qualified Math.NumberTheory.Moduli.Sqrt as Moduli
 import Math.NumberTheory.Moduli.Sqrt (FieldCharacteristic(..))
 import Math.NumberTheory.Powers (integerSquareRoot)
 import Math.NumberTheory.Primes.Types (PrimeNat(..))

Math/NumberTheory/Moduli.hs

@@ -12,11 +12,15 @@
 module Math.NumberTheory.Moduli
   ( module Math.NumberTheory.Moduli.Class
   , module Math.NumberTheory.Moduli.Chinese
+  , module Math.NumberTheory.Moduli.DiscreteLogarithm
   , module Math.NumberTheory.Moduli.Jacobi
+  , module Math.NumberTheory.Moduli.PrimitiveRoot
   , module Math.NumberTheory.Moduli.Sqrt
   ) where
 
 import Math.NumberTheory.Moduli.Chinese
 import Math.NumberTheory.Moduli.Class
+import Math.NumberTheory.Moduli.DiscreteLogarithm
 import Math.NumberTheory.Moduli.Jacobi
+import Math.NumberTheory.Moduli.PrimitiveRoot
 import Math.NumberTheory.Moduli.Sqrt

Math/NumberTheory/Moduli/Class.hs

@@ -29,6 +29,11 @@ module Math.NumberTheory.Moduli.Class
   , invertMod
   , powMod
   , (^%)
+  -- * Multiplicative group
+  , MultMod
+  , multElement
+  , isMultElement
+  , invertGroup
   -- * Unknown modulo
   , SomeMod(..)
   , modulo
@@ -38,8 +43,10 @@ module Math.NumberTheory.Moduli.Class
   , KnownNat
   ) where
 
+
 import Data.Proxy
 import Data.Ratio
+import Data.Semigroup
 import Data.Type.Equality
 import GHC.Integer.GMP.Internals
 import GHC.TypeNats.Compat
@@ -180,6 +187,39 @@ infixr 8 ^%
 -- of type classes in Core.
 -- {-# RULES "^%Mod" forall (x :: KnownNat m => Mod m) p. x ^ p = x ^% p #-}
 
+-- | This type represents elements of the multiplicative group mod m, i.e.
+-- those elements which are coprime to m. Use @toMultElement@ to construct.
+newtype MultMod m = MultMod { multElement :: Mod m }
+  deriving (Eq, Ord, Show)
+
+instance KnownNat m => Semigroup (MultMod m) where
+  MultMod a <> MultMod b = MultMod (a * b)
+  stimes k a@(MultMod a')
+    | k >= 0 = MultMod (powMod a' k)
+    | otherwise = invertGroup $ stimes (-k) a
+  -- ^ This Semigroup is in fact a group, so @stimes@ can be called with a negative first argument.
+
+instance KnownNat m => Monoid (MultMod m) where
+  mempty = MultMod 1
+  mappend = (<>)
+
+instance KnownNat m => Bounded (MultMod m) where
+  minBound = MultMod 1
+  maxBound = MultMod (-1)
+
+-- | Attempt to construct a multiplicative group element.
+isMultElement :: KnownNat m => Mod m -> Maybe (MultMod m)
+isMultElement a = if getVal a `gcd` getMod a == 1
+                     then Just $ MultMod a
+                     else Nothing
+
+-- | For elements of the multiplicative group, we can safely perform the inverse
+-- without needing to worry about failure.
+invertGroup :: KnownNat m => MultMod m -> MultMod m
+invertGroup (MultMod a) = case invertMod a of
+                            Just b -> MultMod b
+                            Nothing -> error "Math.NumberTheory.Moduli.invertGroup: failed to invert element"
+
 -- | 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:

Math/NumberTheory/Moduli/DiscreteLogarithm.hs

@@ -0,0 +1,41 @@
+-- |
+-- Module:       Math.NumberTheory.Moduli.DiscreteLogarithm
+-- Copyright:    (c) 2018 Bhavik Mehta
+-- License:      MIT
+-- Maintainer:   Andrew Lelechenko <[email protected]>
+-- Stability:    Provisional
+-- Portability:  Non-portable
+--
+
+module Math.NumberTheory.Moduli.DiscreteLogarithm where
+
+import Data.Semigroup
+import Data.Maybe
+-- import Data.List
+import qualified Data.IntMap.Strict as M
+import Numeric.Natural
+
+import Math.NumberTheory.Moduli.Class
+import Math.NumberTheory.Moduli.PrimitiveRoot
+import Math.NumberTheory.Prefactored
+import Math.NumberTheory.Powers.Squares
+
+-- | Computes the discrete logarithm. Currently uses a naive search.
+discreteLogarithm
+  :: KnownNat m
+  => PrimitiveRoot m
+  -> MultMod m
+  -> Natural
+-- discreteLogarithm a b = let n = prefValue . groupSize . getGroup $ a
+--                             a' = unPrimitiveRoot a
+--                             vals = genericTake n $ iterate (<> a') mempty
+--                          in fromIntegral $ fromJust $ elemIndex b vals
+
+discreteLogarithm a b = let n = prefValue . groupSize . getGroup $ a
+                            a' = unPrimitiveRoot a
+                            m = integerSquareRoot (n - 1) + 1 -- simple way of ceiling . sqrt
+                            babies = fromInteger . getVal . multElement <$> iterate (<> a') mempty
+                            table = M.fromList $ zip babies [0..(m-1)]
+                            bigGiant = stimes (- toInteger m) a'
+                            giants = fromInteger . getVal . multElement <$> iterate (<> bigGiant) b
+                         in head [i*m + j | (v,i) <- zip giants [0..(m-1)], j <- maybeToList $ M.lookup v table]

Math/NumberTheory/Moduli/PrimitiveRoot.hs

@@ -18,29 +18,36 @@
 {-# LANGUAGE UndecidableInstances #-}
 
 module Math.NumberTheory.Moduli.PrimitiveRoot
-  ( isPrimitiveRoot
-    -- * Cyclic groups
-  , CyclicGroup(..)
+  ( -- * Cyclic groups
+    CyclicGroup(..)
   , cyclicGroupFromModulo
   , cyclicGroupToModulo
+  , groupSize
+    -- * Primitive roots
+  , PrimitiveRoot
+  , unPrimitiveRoot
+  , getGroup
+  , isPrimitiveRoot
   , isPrimitiveRoot'
   ) where
 
-import Control.DeepSeq
 #if __GLASGOW_HASKELL__ < 803
 import Data.Semigroup
 #endif
 
 import Math.NumberTheory.ArithmeticFunctions (totient)
 import Math.NumberTheory.GCD as Coprimes
-import Math.NumberTheory.Moduli (Mod, getNatMod, getNatVal, KnownNat)
+import Math.NumberTheory.Moduli.Class (getNatMod, getNatVal, KnownNat, Mod, MultMod, isMultElement)
 import Math.NumberTheory.Powers.General (highestPower)
 import Math.NumberTheory.Powers.Modular
 import Math.NumberTheory.Prefactored
 import Math.NumberTheory.UniqueFactorisation
 import Math.NumberTheory.Utils.FromIntegral
 
+import Control.DeepSeq
+import Control.Monad (guard)
 import GHC.Generics
+import Numeric.Natural
 
 -- | A multiplicative group of residues is called cyclic,
 -- if there is a primitive root @g@,
@@ -112,6 +119,13 @@ cyclicGroupToModulo = fromFactors . \case
   CGOddPrimePower p k       -> Coprimes.singleton (unPrime p) k
   CGDoubleOddPrimePower p k -> Coprimes.singleton 2 1 <> Coprimes.singleton (unPrime p) k
 
+-- | 'PrimitiveRoot m' is a type which is only inhabited by primitive roots of n.
+data PrimitiveRoot m =
+  PrimitiveRoot { unPrimitiveRoot :: MultMod m -- ^ Extract primitive root value.
+                , getGroup        :: CyclicGroup Natural -- ^ Get cyclic group structure.
+                }
+  deriving (Eq, Show)
+
 -- | 'isPrimitiveRoot'' @cg@ @a@ checks whether @a@ is
 -- a <https://en.wikipedia.org/wiki/Primitive_root_modulo_n primitive root>
 -- of a given cyclic multiplicative group of residues @cg@.
@@ -153,15 +167,21 @@ isPrimitiveRoot' cg r =
 --
 -- Here is how to list all primitive roots:
 --
--- >>> filter isPrimitiveRoot [minBound .. maxBound] :: [Mod 13]
+-- >>> mapMaybe isPrimitiveRoot [minBound .. maxBound] :: [Mod 13]
 -- [(2 `modulo` 13), (6 `modulo` 13), (7 `modulo` 13), (11 `modulo` 13)]
 --
 -- This function is a convenient wrapper around 'isPrimitiveRoot''. The latter
 -- provides better control and performance, if you need them.
 isPrimitiveRoot
   :: KnownNat n
   => Mod n
-  -> Bool
-isPrimitiveRoot r = case cyclicGroupFromModulo (getNatMod r) of
-  Nothing -> False
-  Just cg -> isPrimitiveRoot' cg (getNatVal r)
+  -> Maybe (PrimitiveRoot n)
+isPrimitiveRoot r = do
+  r' <- isMultElement r
+  cg <- cyclicGroupFromModulo (getNatMod r)
+  guard $ isPrimitiveRoot' cg (getNatVal r)
+  return $ PrimitiveRoot r' cg
+
+-- | Calculate the size of a given cyclic group.
+groupSize :: (Integral a, UniqueFactorisation a) => CyclicGroup a -> Prefactored a
+groupSize = totient . cyclicGroupToModulo

arithmoi.cabal

@@ -64,6 +64,7 @@ library
     Math.NumberTheory.Moduli
     Math.NumberTheory.Moduli.Chinese
     Math.NumberTheory.Moduli.Class
+    Math.NumberTheory.Moduli.DiscreteLogarithm
     Math.NumberTheory.Moduli.Equations
     Math.NumberTheory.Moduli.Jacobi
     Math.NumberTheory.Moduli.PrimitiveRoot
@@ -142,6 +143,7 @@ test-suite spec
     Math.NumberTheory.GaussianIntegersTests
     Math.NumberTheory.GCDTests
     Math.NumberTheory.Moduli.ChineseTests
+    Math.NumberTheory.Moduli.DiscreteLogarithmTests
     Math.NumberTheory.Moduli.ClassTests
     Math.NumberTheory.Moduli.EquationsTests
     Math.NumberTheory.Moduli.JacobiTests
@@ -191,6 +193,7 @@ benchmark criterion
       semigroups >=0.8
   other-modules:
     Math.NumberTheory.ArithmeticFunctionsBench
+    Math.NumberTheory.DiscreteLogarithmBench
     Math.NumberTheory.GaussianIntegersBench
     Math.NumberTheory.GCDBench
     Math.NumberTheory.JacobiBench

benchmark/Bench.hs

@@ -3,6 +3,7 @@ module Main where
 import Gauge.Main
 
 import Math.NumberTheory.ArithmeticFunctionsBench as ArithmeticFunctions
+import Math.NumberTheory.DiscreteLogarithmBench as DiscreteLogarithm
 import Math.NumberTheory.GaussianIntegersBench as Gaussian
 import Math.NumberTheory.GCDBench as GCD
 import Math.NumberTheory.JacobiBench as Jacobi
@@ -17,6 +18,7 @@ import Math.NumberTheory.SmoothNumbersBench as SmoothNumbers
 main :: IO ()
 main = defaultMain
   [ ArithmeticFunctions.benchSuite
+  , DiscreteLogarithm.benchSuite
   , Gaussian.benchSuite
   , GCD.benchSuite
   , Jacobi.benchSuite

benchmark/Math/NumberTheory/DiscreteLogarithmBench.hs

@@ -0,0 +1,22 @@
+{-# LANGUAGE DataKinds #-}
+module Math.NumberTheory.DiscreteLogarithmBench 
+  ( benchSuite
+  ) where
+
+import Gauge.Main
+import Data.Maybe
+
+-- import Math.NumberTheory.Moduli (Mod, discreteLogarithm, PrimitiveRoot, isPrimitiveRoot, MultMod, isMultElement)
+import Math.NumberTheory.Moduli.Class (isMultElement)
+import Math.NumberTheory.Moduli.PrimitiveRoot (PrimitiveRoot, isPrimitiveRoot)
+import Math.NumberTheory.Moduli.DiscreteLogarithm (discreteLogarithm)
+
+type Modulus = 1000000007
+
+root :: PrimitiveRoot Modulus
+root = fromJust $ isPrimitiveRoot 5
+
+benchSuite :: Benchmark
+benchSuite = bgroup "Discrete logarithm"
+  [ bench "5^x = 8 mod 10^9+7" $ nf (uncurry discreteLogarithm) (root,fromJust $ isMultElement 8)
+  ]

test-suite/Math/NumberTheory/Moduli/DiscreteLogarithmTests.hs

@@ -0,0 +1,52 @@
+{-# LANGUAGE ScopedTypeVariables #-}
+{-# OPTIONS_GHC -fno-warn-type-defaults #-}
+module Math.NumberTheory.Moduli.DiscreteLogarithmTests
+  ( testSuite
+  ) where
+
+import Data.Maybe
+import Numeric.Natural
+import Test.Tasty
+import Data.Semigroup
+import Data.Proxy
+import GHC.TypeNats.Compat
+
+import Math.NumberTheory.Moduli.Class
+import Math.NumberTheory.Moduli.PrimitiveRoot
+import Math.NumberTheory.Moduli.DiscreteLogarithm
+import Math.NumberTheory.ArithmeticFunctions (totient)
+import Math.NumberTheory.TestUtils
+
+-- | Ensure 'discreteLogarithm' returns in the appropriate range.
+discreteLogRange :: Positive Natural -> Integer -> Integer -> Bool
+discreteLogRange (Positive m) a b =
+  case someNatVal m of
+    SomeNat (_ :: Proxy m) -> fromMaybe True $ do
+      a' <- isPrimitiveRoot (fromInteger a :: Mod m)
+      b' <- isMultElement (fromInteger b)
+      return $ discreteLogarithm a' b' < totient m
+
+-- | Check that 'discreteLogarithm' inverts exponentiation.
+discreteLogarithmProperty :: Positive Natural -> Integer -> Integer -> Bool
+discreteLogarithmProperty (Positive m) a b =
+  case someNatVal m of
+    SomeNat (_ :: Proxy m) -> fromMaybe True $ do
+      a' <- isPrimitiveRoot (fromInteger a :: Mod m)
+      b' <- isMultElement (fromInteger b)
+      return $ discreteLogarithm a' b' `stimes` unPrimitiveRoot a' == b'
+
+-- | Check that 'discreteLogarithm' inverts exponentiation in the other direction.
+discreteLogarithmProperty' :: Positive Natural -> Integer -> Natural -> Bool
+discreteLogarithmProperty' (Positive m) a k =
+  case someNatVal m of
+    SomeNat (_ :: Proxy m) -> fromMaybe True $ do
+      a'' <- isPrimitiveRoot (fromInteger a :: Mod m)
+      let a' = unPrimitiveRoot a''
+      return $ discreteLogarithm a'' (k `stimes` a') == k `mod` totient m
+
+testSuite :: TestTree
+testSuite = testGroup "Discrete logarithm"
+  [ testSmallAndQuick "output is correct range" discreteLogRange
+  , testSmallAndQuick "a^(log_a b) == b"        discreteLogarithmProperty
+  , testSmallAndQuick "log_a a^k == k"          discreteLogarithmProperty'
+  ]

test-suite/Math/NumberTheory/Moduli/PrimitiveRootTests.hs

@@ -21,7 +21,7 @@ import Test.Tasty
 
 import qualified Data.Set as S
 import Data.List (genericTake, genericLength)
-import Data.Maybe (isJust, isNothing)
+import Data.Maybe (isJust, isNothing, mapMaybe)
 import Control.Arrow (first)
 import Numeric.Natural
 import Data.Proxy
@@ -73,34 +73,34 @@ isPrimitiveRootProperty1 :: AnySign Integer -> Positive Natural -> Bool
 isPrimitiveRootProperty1 (AnySign n) (Positive m)
   = case n `modulo` m of
     SomeMod n' -> gcd n (toInteger m) == 1
-               || not (isPrimitiveRoot n')
+               || isNothing (isPrimitiveRoot n')
     InfMod{}   -> False
 
 isPrimitiveRootProperty2 :: Positive Natural -> Bool
 isPrimitiveRootProperty2 (Positive m)
   = isNothing (cyclicGroupFromModulo m)
   || case someNatVal m of
-    SomeNat (_ :: Proxy t) -> any isPrimitiveRoot [(minBound :: Mod t) .. maxBound]
+    SomeNat (_ :: Proxy t) -> not $ null $ mapMaybe isPrimitiveRoot [(minBound :: Mod t) .. maxBound]
 
 isPrimitiveRootProperty3 :: AnySign Integer -> Positive Natural -> Bool
 isPrimitiveRootProperty3 (AnySign n) (Positive m)
   = case n `modulo` m of
-    SomeMod n' -> not (isPrimitiveRoot n')
+    SomeMod n' -> isNothing (isPrimitiveRoot n')
                || allUnique (genericTake (totient m - 1) (iterate (* n') 1))
     InfMod{}   -> False
 
 isPrimitiveRootProperty4 :: AnySign Integer -> Positive Natural -> Bool
 isPrimitiveRootProperty4 (AnySign n) (Positive m)
   = isJust (cyclicGroupFromModulo m)
   || case n `modulo` m of
-    SomeMod n' -> not (isPrimitiveRoot n')
+    SomeMod n' -> isNothing (isPrimitiveRoot n')
     InfMod{}   -> False
 
 isPrimitiveRootProperty5 :: Positive Natural -> Bool
 isPrimitiveRootProperty5 (Positive m)
   = isNothing (cyclicGroupFromModulo m)
   || case someNatVal m of
-       SomeNat (_ :: Proxy t) -> genericLength (filter isPrimitiveRoot [(minBound :: Mod t) .. maxBound]) == totient (totient m)
+       SomeNat (_ :: Proxy t) -> genericLength (mapMaybe isPrimitiveRoot [(minBound :: Mod t) .. maxBound]) == totient (totient m)
 
 testSuite :: TestTree
 testSuite = testGroup "Primitive root"

test-suite/Test.hs

@@ -8,6 +8,7 @@ import qualified Math.NumberTheory.Recurrencies.LinearTests as RecurrenciesLinea
 
 import qualified Math.NumberTheory.Moduli.ChineseTests as ModuliChinese
 import qualified Math.NumberTheory.Moduli.ClassTests as ModuliClass
+import qualified Math.NumberTheory.Moduli.DiscreteLogarithmTests as ModuliDiscreteLogarithm
 import qualified Math.NumberTheory.Moduli.EquationsTests as ModuliEquations
 import qualified Math.NumberTheory.Moduli.JacobiTests as ModuliJacobi
 import qualified Math.NumberTheory.Moduli.PrimitiveRootTests as ModuliPrimitiveRoot
@@ -63,6 +64,7 @@ tests = testGroup "All"
   , testGroup "Moduli"
     [ ModuliChinese.testSuite
     , ModuliClass.testSuite
+    , ModuliDiscreteLogarithm.testSuite
     , ModuliEquations.testSuite
     , ModuliJacobi.testSuite
     , ModuliPrimitiveRoot.testSuite