tests = [
  'q',
  'ord2',
  'div2',
  'mod2',
  'e',
  'R',
  'R0',
  'matrix',
  'gcd',
  'algo',
  'it',
  'double',
  'pos',
  'sz',
  'szscaled',
  'szord',
  'szunique',
  'szmatrix',
]
coverage = {t:0 for t in tests}

def div2(f,g):
  assert f&1
  e = g.ord(2)
  q = ZZ(Integers(2^(e+1))(f)/ZZ(g/2^e))
  return q/2^e

def mod2(f,g):
  assert f&1
  e = g.ord(2)
  q = ZZ(Integers(2^(e+1))(f)/ZZ(g/2^e))
  return ZZ(f-q*g/2^e)

for loop in range(5000):
  e = randrange(1,10)
  f = 1+2*randrange(-2^50,2^50)
  g = 2^e*(1+2*randrange(-2^50,2^50))
  f,g = ZZ(f),ZZ(g)

  numq = 0
  for q in range(1,2^(e+1),2):
    r = q*ZZ(g/2^e)-f
    if r.ord(2) >= e+1:
      numq += 1
  assert numq == 1
  coverage['q'] += 1

  assert e == g.ord(2)
  coverage['ord2'] += 1

  q = ZZ(2^e*div2(f,g))
  r = q*ZZ(g/2^e)-f
  assert r.ord(2) >= e+1
  coverage['div2'] += 1

  assert f == div2(f,g)*g + mod2(f,g)
  coverage['mod2'] += 1

  bitsd = randrange(20)
  d = ZZ(randrange(1-2^bitsd,2^bitsd))
  bits0 = randrange(100)
  bits1 = randrange(100)
  R0 = d*ZZ(randrange(1-2^bits0,2^bits0))
  R1 = 2*d*ZZ(randrange(1-2^bits1,2^bits1))
  G = gcd(R0,R1)

  if R0:
    R = [R0,R1]
    e = [R0.ord(2)]
    while R[-1] != 0:
      e += [R[-1].ord(2)]
      R += [-mod2(ZZ(R[-2]/2^e[-2]),R[-1])/2^e[-1]]
  
    assert len(e) == len(R)-1
    for i in range(len(R)):
      if R[i] != 0:
        assert e[i] == R[i].ord(2)
        coverage['e'] += 1
    for i in range(1,len(R)):
      if R[i] != 0:
        assert R[i+1].ord(2) >= 1
        assert R[i+1] == -mod2(ZZ(R[i-1]/2^e[i-1]),R[i])/2^e[i]
        coverage['R'] += 1
      if R[i] == 0:
        assert i == len(R)-1
        coverage['R0'] += 1
  
    M2Q = MatrixSpace(QQ,2)
    for i in range(1,len(R)-1):
      M = M2Q((0,1/2^e[i],-1/2^e[i],div2(ZZ(R[i-1]/2^e[i-1]),R[i])/2^e[i]))
      assert vector((R[i]/2^e[i],R[i+1])) == M * vector((R[i-1]/2^e[i-1],R[i]))
      coverage['matrix'] += 1

  if R0&1:
    t = len(R)-2
    assert R[t+1] == 0
    assert abs(R[t]/2^e[t]) == G
    coverage['gcd'] += 1

    for i in range(1,len(R)):
      if R[i] != 0:
        assert e[i] == R[i].ord(2)
        q = ZZ((2^e[i]*R[i+1]+R[i-1]/2^e[i-1])/(R[i]/2^e[i]))
        assert q&1
        assert q>0
        assert q<2^(e[i]+1)
        assert (q*R[i]/2^e[i]-R[i-1]/2^e[i-1]).ord(2) >= e[i]+1
        assert R[i+1] == (q*R[i]/2^e[i]-R[i-1]/2^e[i-1])/2^e[i]
        assert R[i+1].ord(2) >= 1
        coverage['algo'] += 1
      else:
        assert i == t+1
        coverage['it'] += 1

  if sign(R0) == sign(R1) and R0 != 0:
    R = [R0,R1]
    e = [R0.ord(2)]
    while len(R) < 100:
      assert R[-1] != 0 # this algorithm never terminates
      e += [R[-1].ord(2)]
      R += [-mod2(ZZ(-R[-2]/2^e[-2]),R[-1])/2^e[-1]]

    for i in range(1,len(e)):
      assert e[i] == R[i].ord(2)
      q = ZZ((2^e[i]*R[i+1]-R[i-1]/2^e[i-1])/(R[i]/2^e[i]))
      assert q&1
      assert q>0
      assert q<2^(e[i]+1)
      assert (q*R[i]/2^e[i]+R[i-1]/2^e[i-1]).ord(2) >= e[i]+1
      assert R[i+1] == (q*R[i]/2^e[i]+R[i-1]/2^e[i-1])/2^e[i]
      assert R[i+1].ord(2) >= 1
      coverage['pos'] += 1

  if R0.ord(2) < R1.ord(2):
    a,b = R0,R1
    while b != 0:
      q = ZZ(Integers(2^(b.ord(2)-a.ord(2)+1))(-a/2^a.ord(2))/ZZ(b/2^b.ord(2)))
      if q > 2^(b.ord(2)-a.ord(2)): q -= 2^(1+b.ord(2)-a.ord(2))
      assert q&1
      assert abs(q) < 2^(b.ord(2)-a.ord(2))
      r = ZZ(a+q*b/2^(b.ord(2)-a.ord(2)))
      assert r.ord(2) >= b.ord(2)+1
      a,b = b,r

    assert abs(a)/2^a.ord(2) == G/2^G.ord(2)
    coverage['sz'] += 1

    a,b = R0,R1
    ab = []
    R = [a]
    while b != 0:
      ab += [(a,b)]
      R += [b]
      q = ZZ(Integers(2^(b.ord(2)-a.ord(2)+1))(-a/2^a.ord(2))/ZZ(b/2^b.ord(2)))
      if q > 2^(b.ord(2)-a.ord(2)): q -= 2^(1+b.ord(2)-a.ord(2))
      assert q&1
      assert abs(q) < 2^(b.ord(2)-a.ord(2))
      r = ZZ(a+q*b/2^(b.ord(2)-a.ord(2)))
      assert r.ord(2) >= b.ord(2)+1
      a,b = ZZ(b/2^b.ord(2)),ZZ(r/2^b.ord(2))

    ab += [(a,b)]
    R += [b]
    assert abs(a)/2^a.ord(2) == G/2^G.ord(2)
    coverage['szscaled'] += 1

    for i in range(len(ab)):
      if R0&1: assert ab[i][0]&1
      if i>0: assert ab[i][0]&1
      if i < len(ab)-1: assert ab[i][1] != 0
      if i == len(ab)-1: assert ab[i][1] == 0
      if i < len(ab)-1:
        e = ab[i][1].ord(2)
        assert ab[i+1][0] == ab[i][1]/2^e
        if R0&1 and e<15:
          numq = 0
          for q in range(1-2^e,2^e,2):
            if ab[i+1][1] == (ab[i][0]+q*ab[i][1]/2^e)/2^e:
              numq += 1
          assert numq == 1
          coverage['szunique'] += 1
      coverage['szord'] += 1

    if R0&1:
      assert R[0] == R0
      assert R[1] == R1
      for i in range(1,len(R)-1):
        e = R[i].ord(2)
        assert e >= 1
        q = ZZ((2^e*R[i+1]-R[i-1]/2^(R[i-1].ord(2)))/(R[i]/2^(R[i].ord(2))))
        assert q&1
        assert q>=1-2^e
        assert q<=2^e-1
        M = M2Q((0,1/2^e,1/2^e,q/4^e))
        assert vector((R[i]/2^R[i].ord(2),R[i+1])) == M * vector((R[i-1]/2^R[i-1].ord(2),R[i]))
        coverage['szmatrix'] += 1

  if R0&1:
    R1 = ZZ(randrange(1-2^bits1,2^bits1))
    assert gcd(R0,R1) == gcd(R0,2*R1)
    coverage['double'] += 1

for t in tests:
  print coverage[t],t