实现指数演算求解DLP

Question

所以我需要实现求解DLP的索引演算算法。主SolveCalculation函数应该只取a、b和素数p并返回x，这样a^x≡ b模p。现在我只有这个问题：

def gen_dlp_prime(nbits):
    p_0 = random_prime(pow(2, nbits), False, pow(2, nbits-1))
    a = 2
    p = a*p_0 + 1
    while p not in Primes():
        a += 1
        p = a * p_0 + 1
        return (p, p_0, a)
def rnd_dlp(nbits):
    p, p_0, a = gen_dlp_prime(nbits)
    g = 2
    while pow(g, p_0, p) != 1:
        g+=1
        x = randrange(2, p)
        return (p, g, x, pow(g, x, p))

def dlp_set_generator(beg, end, step=10):
    print("----- a = b^x (mod p) -----")
    for nbits in range(beg, end, step):
        p, g, x, a = rnd_dlp(nbits)
        print("Nbits = {}\n\tp = {}\n\tg = {}\n\ta = {}\n\tx = {}".format(
            nbits, p,g,a,x))
dlp_set_generator(10, 160)

其中p是素数，g是F_p*的生成元，x是从2到p-1的随机数，因此a=g^x mod p。我的问题是如何实现这个算法？我在网上搜索了一些关于它的书，但我越来越困惑它是如何运作的。也许有人用python或sage实现了？我将非常感激