# pollard's rho algorithm for discrete logarithms
 
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 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 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]) * inverse(j[1]-k[1],p-1)) % (p-1)
    m, l = 0, ((k[2]-j[2]) * inverse(j[1]-k[1],(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

IyBwb2xsYXJkJ3MgcmhvIGFsZ29yaXRobSBmb3IgZGlzY3JldGUgbG9nYXJpdGhtcwoKZGVmIGRsb2coZyx0LHApOgogICAgIyBsIHN1Y2ggdGhhdCBnKipsID09IHQgKG1vZCBwKSwgd2l0aCBwIHByaW1lCiAgICAjIGFsZ29yaXRobSBkdWUgdG8gQ3JhbmRhbGwvUG9tZXJhbmNlICJQcmltZSBOdW1iZXJzIiBzZWMgNS4yLjIKICAgIGZyb20gZnJhY3Rpb25zIGltcG9ydCBnY2QKICAgIGRlZiBpbnZlcnNlKHgsIG0pOgogICAgICAgIGEsIGIsIHUgPSAwLCBtLCAxCiAgICAgICAgd2hpbGUgeCA+IDA6CiAgICAgICAgICAgIHgsIGEsIGIsIHUgPSBiICUgeCwgdSwgeCwgYSAtIGIgLy8geCAqIHUKICAgICAgICBpZiBiID09IDE6IHJldHVybiBhICUgbQogICAgICAgIHJldHVybiAwICMgbXVzdCBiZSBjb3ByaW1lCiAgICBkZWYgZih4YWIpOgogICAgICAgIHgsIGEsIGIgPSB4YWJbMF0sIHhhYlsxXSwgeGFiWzJdCiAgICAgICAgaWYgeCA8IHAvMzoKICAgICAgICAgICAgcmV0dXJuIFsodCp4KSVwLCAoYSsxKSUocC0xKSwgYl0KICAgICAgICBpZiAyKnAvMyA8IHg6CiAgICAgICAgICAgIHJldHVybiBbKGcqeCklcCwgYSwgKGIrMSklKHAtMSldCiAgICAgICAgcmV0dXJuIFsoeCp4KSVwLCAoMiphKSUocC0xKSwgKDIqYiklKHAtMSldCiAgICBpLCBqLCBrID0gMSwgWzEsMCwwXSwgZihbMSwwLDBdKQogICAgd2hpbGUgalswXSA8PiBrWzBdOgogICAgICAgIHByaW50IGksIGosIGsKICAgICAgICBpLCBqLCBrID0gaSsxLCBmKGopLCBmKGYoaykpCiAgICBwcmludCBpLCBqLCBrCiAgICBkID0gZ2NkKGpbMV0gLSBrWzFdLCBwIC0gMSkKICAgIGlmIGQgPT0gMTogcmV0dXJuICgoa1syXS1qWzJdKSAqIGludmVyc2UoalsxXS1rWzFdLHAtMSkpICUgKHAtMSkKICAgIG0sIGwgPSAwLCAoKGtbMl0talsyXSkgKiBpbnZlcnNlKGpbMV0ta1sxXSwocC0xKS9kKSkgJSAoKHAtMSkvZCkKICAgIHdoaWxlIG0gPD0gZDoKICAgICAgICBwcmludCBtLCBsCiAgICAgICAgaWYgcG93KGcsbCxwKSA9PSB0OiByZXR1cm4gbAogICAgICAgIG0sIGwgPSBtKzEsIChsKygocC0xKS9kKSklKHAtMSkKICAgIHJldHVybiBGYWxzZQoKcHJpbnQgImRsb2coODMsNTU1LDk5NykgPT4gMTI5IgpwcmludCBkbG9nKDgzLDU1NSw5OTcpICMgMTI5CnByaW50CnByaW50ICJkbG9nKDgzLDU2Niw5OTcpID0+IDE0NyIKcHJpbnQgZGxvZyg4Myw1NjYsOTk3KSAjIDE0Nw==