"""
(axy + bx + cy + d)^n

Reduction rules:

x^2 = (x + 3)
y^2 = 2y
constants a, b, c, and d are computed mod p

(a x y + b x + c y + d) (e x y + f x + g y + h)
= dexy + bex^2y + cexy^2 + aex^2y^2 + dfx + bfx^2 + cfxy + afx^2y + dgy + bgxy + cgy^2 + agy^2x + dh + bhx + chy + ahxy
= dexy + be(x+3)y + 2cexy + 2ae(x+3)y + dfx + bf(x+3) + cfxy + af(x+3)y + dgy + bgxy + 2cgy + 2agyx + dh + bhx + chy + ahxy
= dexy + bexy + be3y + 2cexy + 2aexy + 2ae3y + dfx + bfx + bf3 + cfxy + afxy + af3y + dgy + bgxy + 2cgy + 2agyx + dh + bhx + chy + ahxy
= dexy + bexy + 3bey + 2cexy + 2aexy + 6aey + dfx + bfx + 3bf + cfxy + afxy + 3afy + dgy + bgxy + 2cgy + 2agyx + dh + bhx + chy + ahxy
= xy(de + be + 2ce + 2ae + cf + af + bg + 2ag + ah) + x(df + bf + bh) + y(3be + 6ae + 3af + dg + 2cg + ch) + (3bf + dh)
"""

def mul(m, n, p):
    a, b, c, d = m[0], m[1], m[2], m[3]
    e, f, g, h = n[0], n[1], n[2], n[3]
    #xy(de + be + 2ce + 2ae + cf + af + bg + 2ag + ah)
    A = (d*e + b*e + 2*c*e + 2*a*e + c*f + a*f + b*g + 2*a*g + a*h) % p
    #x(df + bf + bh)
    B = (d*f + b*f + b*h) % p
    #y(3be + 6ae + 3af + dg + 2cg + ch)
    C = (3*b*e + 6*a*e + 3*a*f + d*g + 2*c*g + c*h) % p
    #(3bf + dh)
    D = (3*b*f + d*h) % p
    return (A, B, C, D)

def exp(g, n, p):
    r = (0, 0, 0, 1)
    while n != 0:
        if n & 1:
            r = mul(r, g, p)
        n >>= 1
        g = mul(g, g, p)
    return r

def testBase(g, p):
    s = set()
    r1 = (0, 0, 0, 1)
    r2 = r1
    i = 0
    while True:
        e = exp(g, i, p)
        #print i, ":  ", r1, "=", r2, "=", e
        if r1 != r2 or r1 != e:
            raise Exception("Mismatch")
        if r1 in s:
            print "length", i
            return i
        s.add(r1)
        r1 = mul(r1, g, p)
        r2 = mul(g, r2, p)
        i += 1

p = 11

for i in range(1, 100):
    g = (i, 0, 0, 0)
    testBase(g, p)
