We will perform the Index calculus algorithm to compute $\log_{a} b_i$ ($i = 1, 2$) in $\mathbf{Z}_{p}^*$, where $a = 5$, $b_1 = 4389733$, $b_2 = 1234567$, and $p = 9330887$ is a prime number.
from algorithms.euclidean import gcd
from algorithms.factorization import factorizeByBase
p = 9330887
a = 5
b1 = 4389733
b2 = 1234567
We will use a base consisting of $-1$ and all primes less than $50$.
base = [-1] + prime_range(50)
base
[-1, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
Let us now try to find the logarithms table. We construct a matrix by trying random powers of $a$ and factoring them with the numbers in our base. We stop when the matrix has full rank.
set_random_seed(0)
v = []
s = set()
r, l = 0, len(base)
M = Matrix(nrows=0, ncols=l)
while r < l:
i = randint(1, p-1)
if i in s:
continue
s.add(i)
x = pow(a, i, p)
f = factorizeByBase(Integer(x), base, p)
if not f:
continue
MM = Matrix(list(M) + [f])
rr = MM.rank()
if rr > r:
M = MM
r = rr
v.append(i)
print(M)
print(v)
[1 0 2 1 0 0 1 0 0 0 0 0 0 0 2 0] [1 0 1 0 0 0 0 0 0 0 0 1 3 0 0 0] [0 3 8 0 1 0 0 0 0 1 0 0 0 0 0 0] [0 1 2 1 2 0 0 0 0 0 1 0 1 0 0 0] [1 0 0 0 1 3 1 0 0 0 0 0 0 1 0 0] [1 4 0 0 1 0 0 1 0 0 1 0 0 0 0 0] [1 2 0 3 3 0 0 0 0 1 0 0 0 0 0 0] [1 2 2 0 0 1 1 1 0 1 0 0 0 0 0 0] [0 8 1 0 0 2 0 0 0 0 0 0 0 0 1 0] [0 3 1 0 0 1 0 0 1 1 0 0 0 0 0 0] [0 5 0 1 0 0 0 1 0 0 1 0 1 0 0 0] [0 0 0 1 1 0 0 0 1 0 0 0 0 2 0 0] [0 0 0 0 2 0 0 0 1 0 1 0 0 0 0 1] [0 1 3 0 4 0 1 0 0 0 0 0 0 0 0 0] [1 4 1 5 0 0 0 0 0 1 0 0 0 0 0 0] [1 1 1 2 0 2 1 0 1 0 0 0 0 0 0 0] [1342485, 669118, 7237936, 7922551, 9149713, 4288717, 563606, 4482702, 7606125, 6167680, 7916700, 6825736, 4712055, 7120183, 7905940, 4584527]
We now solve the system of equations we have obtained. The solution represents our table of logarithms and can be used to find logarithms of multiple numbers.
t = M^-1 * Matrix(zip(v)) % (p-1)
t
[4665443] [6670912] [2030334] [ 1] [5786904] [1078197] [8534197] [7606749] [8519903] [2519168] [6200403] [9068634] [7409417] [5590350] [6037417] [6410599]
Let us verify that the table is correct.
[pow(a, i, p) for i, in t]
[9330886, 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47]
We can now find a power of $b_i$ that can be factorized with our base. We will use the following function.
def findPower(b, p, base):
while True:
i = randint(1, p-1)
x = pow(b, i, p)
f = factorizeByBase(Integer(x), base, p)
if not f:
continue
return (i, f)
Let us find such a power of $b_1$ and its factorization.
i1, f1 = findPower(b1, p, base)
i1, f1
(3157805, [0, 0, 3, 1, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0])
We can now use the logarithm table to obtain a congruence in $\log_a b_1$. However, there may be multiple solutions, and we need to check which one is our answer. The number of solutions is given by the greatest common divisor of $p-1$ and the obtained exponent.
g1 = gcd(i1, p-1)
g1
1
We can thus obtain a solution modulo $p-1$ divided by this GCD.
m1 = ((p-1)/g1)
(x1,), = Matrix([f1])*t / i1 % m1
x1
5753305
Let us now put this into a function which will also verify the potential solutions.
def checkSolutions(a, b, p, t, i, f):
g = gcd(i, p-1)
m = (p-1)/g
(x,), = Matrix([f])*t / i % m
for j in range(g):
y = pow(a, x, p)
print("%d^%d mod %d = %d" % (a, x, p, y))
if y == b:
return x
x += m
checkSolutions(a, b1, p, t, i1, f1)
5^5753305 mod 9330887 = 4389733
5753305
We have thus computed $\log_a b_1$. Let us now compute $\log_a b_2$.
i2, f2 = findPower(b2, p, base)
i2, f2
(5353502, [0, 0, 2, 1, 0, 0, 2, 1, 0, 0, 0, 0, 0, 1, 0, 0])
checkSolutions(a, b2, p, t, i2, f2)
5^2608331 mod 9330887 = 8096320 5^7273774 mod 9330887 = 1234567
7273774
We will now use the function logarithmTable
to compute tables of logarithms, which will then be used to compute some more discrete logarithms with the function indexCalculus
. We will measure the evaluation times.
from algorithms.discreteLogarithm import logarithmTable, indexCalculus
aa = 47
bb = 191
pp = 100000000003
table = %time logarithmTable(aa, pp, base, trace = True)
print(table)
%time indexCalculus(aa, bb, pp, base, table, trace = True)
found factorization 47^8920294660 = -1^1 * 2^2 * 3^2 * 7^1 * 11^1 * 13^1 * 19^1 * 41^1 * 43^1 * 47^1 (mod 100000000003) found factorization 47^86731650251 = -1^1 * 2^2 * 3^2 * 7^1 * 11^3 * 13^2 * 19^1 * 31^1 (mod 100000000003) found factorization 47^39752784926 = 2^7 * 3^1 * 5^1 * 11^3 * 19^2 * 23^1 (mod 100000000003) found factorization 47^69419634802 = 2^1 * 3^4 * 5^7 * 19^1 * 37^1 (mod 100000000003) found factorization 47^75320357659 = -1^1 * 3^1 * 5^3 * 7^1 * 11^2 * 13^2 * 31^1 * 47^1 (mod 100000000003) found factorization 47^40532049132 = -1^1 * 2^4 * 5^1 * 11^1 * 19^1 * 31^1 * 37^1 * 41^1 * 47^1 (mod 100000000003) found factorization 47^32937591819 = 3^3 * 5^1 * 11^4 * 13^1 * 31^2 (mod 100000000003) found factorization 47^6602367151 = -1^1 * 2^1 * 7^2 * 17^1 * 29^4 * 43^1 (mod 100000000003) found factorization 47^70383598779 = -1^1 * 2^1 * 11^1 * 17^1 * 31^2 * 37^1 * 41^1 * 43^1 (mod 100000000003) found factorization 47^82842048017 = 5^1 * 13^3 * 17^1 * 29^1 (mod 100000000003) found factorization 47^93452946324 = 2^3 * 3^1 * 5^1 * 11^1 * 19^1 * 23^1 * 29^1 * 31^1 * 41^1 (mod 100000000003) found factorization 47^85301660587 = 3^1 * 11^2 * 19^4 * 23^1 * 47^1 (mod 100000000003) found factorization 47^95880313549 = -1^1 * 2^1 * 3^4 * 13^4 * 23^3 (mod 100000000003) found factorization 47^32201305943 = 2^9 * 11^3 * 17^1 * 31^1 * 41^1 (mod 100000000003) found factorization 47^26135740991 = -1^1 * 2^2 * 3^8 * 5^1 * 7^1 * 11^2 * 13^1 * 37^1 (mod 100000000003) found factorization 47^47177975199 = -1^1 * 2^5 * 7^2 * 11^1 * 13^1 * 23^1 * 31^1 * 47^1 (mod 100000000003) CPU times: user 2.47 s, sys: 10.5 ms, total: 2.48 s Wall time: 2.51 s [50000000001, 88948660783, 74403602789, 76906212201, 2995382539, 11520765452, 33861633403, 95303711677, 78308239391, 74623569337, 9047223934, 26437938024, 66204838085, 75359412859, 30170333882, 1] found factorization 191^92538805852 = -1^1 * 2^2 * 3^3 * 23^2 * 29^1 * 47^2 (mod 100000000003) checking 2 solutions for 47^x = 191, x = 6935101882 (mod 50000000001) 47^6935101882 mod 100000000003 = 191 CPU times: user 86.4 ms, sys: 0 ns, total: 86.4 ms Wall time: 89.4 ms
6935101882
pp = 10000000000000000051
base = [-1] + prime_range(5000)
table = %time logarithmTable(aa, pp, base, trace = True)
%time indexCalculus(aa, bb, pp, base, table, trace = True)