Merge pull request #26488 from teschlg/teschlg-patch-2

Implementation of Index Calculus for discrete_log

sympy/ntheory/residue_ntheory.py

@@ -6,6 +6,7 @@
 from sympy.polys.domains import ZZ
 from sympy.polys.galoistools import gf_crt1, gf_crt2, linear_congruence, gf_csolve
 from .primetest import isprime
+from .generate import primerange
 from .factor_ import factorint, _perfect_power
 from .modular import crt
 from sympy.utilities.decorator import deprecated
@@ -1432,6 +1433,117 @@ def _discrete_log_pollard_rho(n, a, b, order=None, retries=10, rseed=None):
     raise ValueError("Pollard's Rho failed to find logarithm")
 
 
+def _discrete_log_is_smooth(n: int, factorbase: list):
+    """Try to factor n with respect to a given factorbase.
+    Upon success a list of exponents with repect to the factorbase is returned.
+    Otherwise None."""
+    factors = [0]*len(factorbase)
+    for i, p in enumerate(factorbase):
+        while n % p == 0: # divide by p as many times as possible
+            factors[i] += 1
+            n = n // p
+    if n != 1:
+        return None # the number factors if at the end nothing is left
+    return factors
+
+
+def _discrete_log_index_calculus(n, a, b, order, rseed=None):
+    """
+    Index Calculus algorithm for computing the discrete logarithm of ``a`` to
+    the base ``b`` modulo ``n``.
+
+    The group order must be given and prime. It is not suitable for small orders
+    and the algorithm might fail to find a solution in such situations.
+
+    Examples
+    ========
+
+    >>> from sympy.ntheory.residue_ntheory import _discrete_log_index_calculus
+    >>> _discrete_log_index_calculus(24570203447, 23859756228, 2, 12285101723)
+    4519867240
+
+    See Also
+    ========
+
+    discrete_log
+
+    References
+    ==========
+
+    .. [1] "Handbook of applied cryptography", Menezes, A. J., Van, O. P. C., &
+        Vanstone, S. A. (1997).
+    """
+    randint = _randint(rseed)
+    from math import sqrt, exp, log
+    a %= n
+    b %= n
+    # assert isprime(order), "The order of the base must be prime."
+    # First choose a heuristic the bound B for the factorbase.
+    # We have added an extra term to the asymptotic value which
+    # is closer to the theoretical optimum for n up to 2^70.
+    B = int(exp(0.5 * sqrt( log(n) * log(log(n)) )*( 1 + 1/log(log(n)) )))
+    max = 5 * B * B  # expected number of trys to find a relation
+    factorbase = list(primerange(B)) # compute the factorbase
+    lf = len(factorbase) # length of the factorbase
+    ordermo = order-1
+    abx = a
+    for x in range(order):
+        if abx == 1:
+            return (order - x) % order
+        relationa = _discrete_log_is_smooth(abx, factorbase)
+        if relationa:
+            relationa = [r % order for r in relationa] + [x]
+            break
+        abx = abx * b % n # abx = a*pow(b, x, n) % n
+
+    else:
+        raise ValueError("Index Calculus failed")
+
+    relations = [None] * lf
+    k = 1  # number of relations found
+    kk = 0
+    while k < 3 * lf and kk < max:  # find relations for all primes in our factor base
+        x = randint(1,ordermo)
+        relation = _discrete_log_is_smooth(pow(b,x,n), factorbase)
+        if relation is None:
+            kk += 1
+            continue
+        k += 1
+        kk = 0
+        relation += [ x ]
+        index = lf  # determine the index of the first nonzero entry
+        for i in range(lf):
+            ri = relation[i] % order
+            if ri> 0 and relations[i] is not None:  # make this entry zero if we can
+                for j in range(lf+1):
+                    relation[j] = (relation[j] - ri*relations[i][j]) % order
+            else:
+                relation[i] = ri
+            if relation[i] > 0 and index == lf:  # is this the index of the first nonzero entry?
+                index = i
+        if index == lf or relations[index] is not None:  # the relation contains no new information
+            continue
+        # the relation contains new information
+        rinv = pow(relation[index],-1,order)  # normalize the first nonzero entry
+        for j in range(index,lf+1):
+            relation[j] = rinv * relation[j] % order
+        relations[index] = relation
+        for i in range(lf):  # subtract the new relation from the one for a
+            if relationa[i] > 0 and relations[i] is not None:
+                rbi = relationa[i]
+                for j in range(lf+1):
+                    relationa[j] = (relationa[j] - rbi*relations[i][j]) % order
+            if relationa[i] > 0:  # the index of the first nonzero entry
+                break  # we do not need to reduce further at this point
+        else:  # all unkowns are gone
+            #print(f"Success after {k} relations out of {lf}")
+            x = (order -relationa[lf]) % order
+            if pow(b,x,n) == a:
+                return x
+            raise ValueError("Index Calculus failed")
+    raise ValueError("Index Calculus failed")
+
+
 def _discrete_log_pohlig_hellman(n, a, b, order=None, order_factors=None):
     """
     Pohlig-Hellman algorithm for computing the discrete logarithm of ``a`` to
@@ -1494,6 +1606,7 @@ def discrete_log(n, a, b, order=None, prime_order=None):
         * Trial multiplication
         * Baby-step giant-step
         * Pollard's Rho
+        * Index Calculus
         * Pohlig-Hellman
 
     Examples
@@ -1511,6 +1624,7 @@ def discrete_log(n, a, b, order=None, prime_order=None):
         Vanstone, S. A. (1997).
 
     """
+    from math import sqrt, log
     n, a, b = as_int(n), as_int(a), as_int(b)
     if order is None:
         # Compute the order and its factoring in one pass
@@ -1550,11 +1664,16 @@ def discrete_log(n, a, b, order=None, prime_order=None):
     if order < 1000:
         return _discrete_log_trial_mul(n, a, b, order)
     elif prime_order:
-        if order < 1000000000000:
+        # Shanks and Pollard rho are O(sqrt(order)) while index calculus is O(exp(2*sqrt(log(n)log(log(n)))))
+        # we compare the expected running times to determine the algorithmus which is expected to be faster
+        if 4*sqrt(log(n)*log(log(n))) < log(order) - 10:  # the number 10 was determined experimental
+            return _discrete_log_index_calculus(n, a, b, order)
+        elif order < 1000000000000:
+            # Shanks seems typically faster, but uses O(sqrt(order)) memory
             return _discrete_log_shanks_steps(n, a, b, order)
-        return _discrete_log_pollard_rho(n, a, b, order, order_factors)
+        return _discrete_log_pollard_rho(n, a, b, order)
 
-    return _discrete_log_pohlig_hellman(n, a, b, order)
+    return _discrete_log_pohlig_hellman(n, a, b, order, order_factors)
 
 
 

sympy/ntheory/tests/test_residue.py

@@ -11,7 +11,7 @@
 from sympy.ntheory.residue_ntheory import _primitive_root_prime_iter, \
     _primitive_root_prime_power_iter, _primitive_root_prime_power2_iter, \
     _nthroot_mod_prime_power, _discrete_log_trial_mul, _discrete_log_shanks_steps, \
-    _discrete_log_pollard_rho, _discrete_log_pohlig_hellman, \
+    _discrete_log_pollard_rho, _discrete_log_index_calculus, _discrete_log_pohlig_hellman, \
     _binomial_mod_prime_power, binomial_mod
 from sympy.polys.domains import ZZ
 from sympy.testing.pytest import raises
@@ -256,7 +256,11 @@ def test_residue():
     assert _discrete_log_pollard_rho(24567899, 3**333, 3, rseed=0) == 333
     raises(ValueError, lambda: _discrete_log_pollard_rho(11, 7, 31, rseed=0))
     raises(ValueError, lambda: _discrete_log_pollard_rho(227, 3**7, 5, rseed=0))
-
+    assert _discrete_log_index_calculus(983, 948, 2, 491) == 183
+    assert _discrete_log_index_calculus(633383, 21794, 2, 316691) == 68048
+    assert _discrete_log_index_calculus(941762639, 68822582, 2, 470881319) == 338029275
+    assert _discrete_log_index_calculus(999231337607, 888188918786, 2, 499615668803) == 142811376514
+    assert _discrete_log_index_calculus(47747730623, 19410045286, 43425105668, 645239603) == 590504662
     assert _discrete_log_pohlig_hellman(98376431, 11**9, 11) == 9
     assert _discrete_log_pohlig_hellman(78723213, 11**31, 11) == 31
     assert _discrete_log_pohlig_hellman(32942478, 11**98, 11) == 98
@@ -265,6 +269,10 @@ def test_residue():
     assert discrete_log(2456747, 3**51, 3) == 51
     assert discrete_log(32942478, 11**127, 11) == 127
     assert discrete_log(432751500361, 7**324, 7) == 324
+    assert discrete_log(265390227570863,184500076053622, 2) == 17835221372061
+    assert discrete_log(22708823198678103974314518195029102158525052496759285596453269189798311427475159776411276642277139650833937,
+                        17463946429475485293747680247507700244427944625055089103624311227422110546803452417458985046168310373075327,
+                        123456) == 2068031853682195777930683306640554533145512201725884603914601918777510185469769997054750835368413389728895
     args = 5779, 3528, 6215
     assert discrete_log(*args) == 687
     assert discrete_log(*Tuple(*args)) == 687