# pollard's rho algorithm for discrete logarithms
 
def inverse(x, m):
    a, b, u = 0, m, 1
    while x > 0:
        x, a, b, u = b % x, u, x, a - b // x * u
    if b == 1: return a % m
    return 0 # must be coprime
 
def dlog(g,t,p):
    # l such that g**l == t (mod p), with p prime
    # algorithm due to Crandall/Pomerance "Prime Numbers" sec 5.2.2
    from fractions import gcd
    def f(xab):
        x, a, b = xab[0], xab[1], xab[2]
        if x < p/3:
            return [(t*x)%p, (a+1)%(p-1), b]
        if 2*p/3 < x:
            return [(g*x)%p, a, (b+1)%(p-1)]
        return [(x*x)%p, (2*a)%(p-1), (2*b)%(p-1)]
    i, j, k = 1, [1,0,0], f([1,0,0])
    while j[0] <> k[0]:
        print i, j, k
        i, j, k = i+1, f(j), f(f(k))
    print i, j, k
    d = gcd(j[1] - k[1], p - 1)
    if d == 1: return ((k[2]-j[2])%(p-1) * inverse((j[1]-k[1])%(p-1),p-1)) % (p-1)
    m, l = 0, ((k[2]-j[2])%((p-1)/d) * inverse((j[1]-k[1])%((p-1)/d),(p-1)/d)) % ((p-1)/d)
    while m <= d:
        print m, l
        if pow(g,l,p) == t: return l
        m, l = m+1, (l+((p-1)/d))%(p-1)
    return False
 
print "dlog(83,555,997) => 129"
print dlog(83,555,997) # 129
print
print "dlog(83,566,997) => 147"
print dlog(83,566,997) # 147

IyBwb2xsYXJkJ3MgcmhvIGFsZ29yaXRobSBmb3IgZGlzY3JldGUgbG9nYXJpdGhtcwoKZGVmIGludmVyc2UoeCwgbSk6CiAgICBhLCBiLCB1ID0gMCwgbSwgMQogICAgd2hpbGUgeCA+IDA6CiAgICAgICAgeCwgYSwgYiwgdSA9IGIgJSB4LCB1LCB4LCBhIC0gYiAvLyB4ICogdQogICAgaWYgYiA9PSAxOiByZXR1cm4gYSAlIG0KICAgIHJldHVybiAwICMgbXVzdCBiZSBjb3ByaW1lCgpkZWYgZGxvZyhnLHQscCk6CiAgICAjIGwgc3VjaCB0aGF0IGcqKmwgPT0gdCAobW9kIHApLCB3aXRoIHAgcHJpbWUKICAgICMgYWxnb3JpdGhtIGR1ZSB0byBDcmFuZGFsbC9Qb21lcmFuY2UgIlByaW1lIE51bWJlcnMiIHNlYyA1LjIuMgogICAgZnJvbSBmcmFjdGlvbnMgaW1wb3J0IGdjZAogICAgZGVmIGYoeGFiKToKICAgICAgICB4LCBhLCBiID0geGFiWzBdLCB4YWJbMV0sIHhhYlsyXQogICAgICAgIGlmIHggPCBwLzM6CiAgICAgICAgICAgIHJldHVybiBbKHQqeCklcCwgKGErMSklKHAtMSksIGJdCiAgICAgICAgaWYgMipwLzMgPCB4OgogICAgICAgICAgICByZXR1cm4gWyhnKngpJXAsIGEsIChiKzEpJShwLTEpXQogICAgICAgIHJldHVybiBbKHgqeCklcCwgKDIqYSklKHAtMSksICgyKmIpJShwLTEpXQogICAgaSwgaiwgayA9IDEsIFsxLDAsMF0sIGYoWzEsMCwwXSkKICAgIHdoaWxlIGpbMF0gPD4ga1swXToKICAgICAgICBwcmludCBpLCBqLCBrCiAgICAgICAgaSwgaiwgayA9IGkrMSwgZihqKSwgZihmKGspKQogICAgcHJpbnQgaSwgaiwgawogICAgZCA9IGdjZChqWzFdIC0ga1sxXSwgcCAtIDEpCiAgICBpZiBkID09IDE6IHJldHVybiAoKGtbMl0talsyXSklKHAtMSkgKiBpbnZlcnNlKChqWzFdLWtbMV0pJShwLTEpLHAtMSkpICUgKHAtMSkKICAgIG0sIGwgPSAwLCAoKGtbMl0talsyXSklKChwLTEpL2QpICogaW52ZXJzZSgoalsxXS1rWzFdKSUoKHAtMSkvZCksKHAtMSkvZCkpICUgKChwLTEpL2QpCiAgICB3aGlsZSBtIDw9IGQ6CiAgICAgICAgcHJpbnQgbSwgbAogICAgICAgIGlmIHBvdyhnLGwscCkgPT0gdDogcmV0dXJuIGwKICAgICAgICBtLCBsID0gbSsxLCAobCsoKHAtMSkvZCkpJShwLTEpCiAgICByZXR1cm4gRmFsc2UKCnByaW50ICJkbG9nKDgzLDU1NSw5OTcpID0+IDEyOSIKcHJpbnQgZGxvZyg4Myw1NTUsOTk3KSAjIDEyOQpwcmludApwcmludCAiZGxvZyg4Myw1NjYsOTk3KSA9PiAxNDciCnByaW50IGRsb2coODMsNTY2LDk5NykgIyAxNDc=