#from math import *

# from stackexchange
def egcd(a, b):
    if a == 0:
        return (b, 0, 1)
    else:
        g, y, x = egcd(b % a, a)
        return (g, x - (b // a) * y, y)

def modinv(a, m):
    a %= m
    if a == 0:
        import pdb; pdb.set_trace()
    g, x, y = egcd(a, m)
    if g != 1:
        raise Exception('modular inverse for %d does not exist mod %d' % (a, m))
    else:
        return x % m

# Jacobi and msqrt copied from python-ecdsa
def jacobi(a, n):
    assert n >= 3
    assert n%2 == 1
    a = a % n
    if a == 0:
        return 0
    if a == 1:
        return 1
    a1, e = a, 0
    while a1%2 == 0:
        a1, e = a1//2, e+1
    if e%2 == 0 or n%8 == 1 or n%8 == 7:
        s = 1
    else:
        s = -1
    if a1 == 1:
        return s
    if n%4 == 3 and a1%4 == 3:
        s = -s
    return s * jacobi(n % a1, a1)

def msqrt(a, p):
    a %= p
    if a == 0:
        return 0
    if p == 2:
        return a
    jac = jacobi(a, p)
    if jac == -1:
        raise Exception("%d has no square root modulo %d" % (a, p))
    if p % 4 == 3:
        return pow(a, (p+1)//4, p)
    if p % 8 == 5:
        d = pow(a, (p-1)//4, p)
    if d == 1:
        return pow(a, (p+3)//8, p)
    if d == p-1:
        return (2 * a * pow(4*a, (p-5)//8, p)) % p
    raise Exception("Shouldn't get here.")

# This is likely an error on Wikipedia
def findBsq(m, n):
    return modinv(1-m, n)

# This is likely an error on Wikipedia
def findBinvSq(m, n):
    return (1-m) % n

def findDN(cn, sn, m, n):
    # This matches other sources.  Wikipedia is wrong here.
    return msqrt(cn**2 + sn**2*findBinvSq(m, n), n)

class Point:

    def __init__(self, x, y, cn, sn, dn, m, n):
        self.x = x
        self.y = y
        self.cn = cn
        self.sn = sn
        self.dn = dn
        self.m = m
        self.n = n
        self.verifyPoint()

    def verifyPoint(self):
        x, y, cn, sn, dn, m, n = self.x, self.y, self.cn, self.sn, self.dn, self.m, self.n
        binvSq = findBinvSq(m, n)
        if(x != sn):
            import pdb; pdb.set_trace()
            raise Exception("Edwards and ellips points are different!")
        if ((cn**2 + sn**2) % n != 1 or
            (cn**2 + sn**2*binvSq) % n != dn**2 % n or
            (x**2 + y**2) % n != (1 + m*x**2*y**2) % n):
            raise Exception("Invalid point %s" % self)

    def __str__(self):
        return "(x:%d, y:%d) (sn:%d, cn:%d dn:%d) m:%d n:%d" % (
            self.x, self.y, self.sn, self.cn, self.dn, self.m, self.n)

def pointsEqual(a, b):
    return a.x == b.x and a.y == b.y

def genPointFromX(x, upper, m, n):
    y = msqrt((1 - x**2)*modinv(1 - m*x**2, n), n)
    if upper != (y >= n//2):
        y = -y % n
    sn = x
    cn = msqrt(1 - sn**2, n)
    if upper != (cn >= n//2):
        cn = -cn % n
    dn = findDN(cn, sn, m, n)
    return Point(x, y, cn, sn, dn, m, n)

def ecAdd(a, b):
    c1, s1, d1, x1, y1 = a.cn, a.sn, a.dn, a.x, a.y
    c2, s2, d2, x2, y2 = b.cn, b.sn, b.dn, b.x, b.y
    if a.m != b.m or a.n != b.n:
        raise Exception("Points are on different curves")
    m, n = a.m, a.n

    # Add the points on the ellipse
    Ninv = modinv(1 - m*s1**2*s2**2, n)
    c3 = (c1*c2 - s1*s2*d1*d2) * Ninv % n
    s3 = (s1*c2*d2 + s2*c1*d1) * Ninv % n
    d3 = (d1*d2 - m*s1*s2*c1*c2) * Ninv % n

    # Add the points on the Edwards curve
    x3 = (x1*y2 + x2*y1)*modinv(1 + m*x1*x2*y1*y2, n) % n
    y3 = (y1*y2 - x1*x2)*modinv(1 - m*x1*x2*y1*y2, n) % n
    return Point(x3, y3, c3, s3, d3, m, n)

def doubleEC(p):
    return addEC(p, p)

def mulEC(m, p):
    r = Point(0, 1, 1, 0, 1, m, p.n)
    while m != 0:
        if m & 1:
            r = addEC(r, p)
        m >>= 1
        p = doubleEC(p)
    return r

#n = 29
#x = 8
#m = 9*modinv(10, n) % n
n = 7883
x = 17
m = 11*modinv(12, n) % n
print "Bsq", findBsq(m, n)
g = genPointFromX(x, True, m, n)
# Create the origin point
origin = Point(0, 1, 1, 0, 1, m, n)
p = g
i = 1
while not pointsEqual(p, origin):
    print i, p
    p = ecAdd(p, g)
    i += 1