#!/usr/bin/env python
# coding: utf-8

# # Index calculus algorithm
# 
# 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.

# In[1]:


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$.

# In[3]:


base = [-1] + prime_range(50)
base


# 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.

# In[4]:


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)


# 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.

# In[5]:


t = M^-1 * Matrix(zip(v)) % (p-1)
t


# Let us verify that the table is correct.

# In[6]:


[pow(a, i, p) for i, in t]


# We can now find a power of $b_i$ that can be factorized with our base. We will use the following function.

# In[7]:


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.

# In[8]:


i1, f1 = findPower(b1, p, base)
i1, f1


# 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.

# In[9]:


g1 = gcd(i1, p-1)
g1


# We can thus obtain a solution modulo $p-1$ divided by this GCD.

# In[10]:


m1 = ((p-1)/g1)
(x1,), = Matrix([f1])*t / i1 % m1
x1


# Let us now put this into a function which will also verify the potential solutions.

# In[11]:


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)


# We have thus computed $\log_a b_1$. Let us now compute $\log_a b_2$.

# In[12]:


i2, f2 = findPower(b2, p, base)
i2, f2


# In[13]:


checkSolutions(a, b2, p, t, i2, f2)


# ## Running time comparison
# 
# 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.

# In[14]:


from algorithms.discreteLogarithm import logarithmTable, indexCalculus

aa = 47
bb = 191


# In[15]:


pp = 100000000003
table = get_ipython().run_line_magic('time', 'logarithmTable(aa, pp, base, trace = True)')
print(table)
get_ipython().run_line_magic('time', 'indexCalculus(aa, bb, pp, base, table, trace = True)')


# In[ ]:


pp = 10000000000000000051
base = [-1] + prime_range(5000)
table = get_ipython().run_line_magic('time', 'logarithmTable(aa, pp, base, trace = True)')
get_ipython().run_line_magic('time', 'indexCalculus(aa, bb, pp, base, table, trace = True)')