import random

tests = [
  'posdivstep',
  'pmdivstep',
  'pmdivstep1',
  'pmdivstep2',
  'pmdivstep3',
  'T',
  'S',
  'decomposition',
  'range',
  'range2',
  'range3',
  '3divstep',
  '3cdivstep',
]
coverage = {t:0 for t in tests}

def divstep(delta,f,g):
  if delta>0 and g&1:
    return 1-delta,g,ZZ((g-f)/2)
  return 1+delta,f,ZZ((g+(g&1)*f)/2)

def posdivstep(delta,f,g):
  if delta>0 and g&1:
    return 1-delta,g,ZZ((g+f)/2)
  return 1+delta,f,ZZ((g+(g&1)*f)/2)

def cdivstep(delta,f,g):
  if delta>0 and g&1:
    return 1-delta,g,ZZ((f-g)/2)
  if delta == 0:
    return 1+delta,f,ZZ((g+(g&1)*f)/2)
  return 1+delta,f,ZZ((g-(g&1)*f)/2)

def pmdivstep(delta,f,g):
  if delta>=0 and (g-f)%4 == 0: return -delta,g,ZZ((f-g)/2)
  if delta>=0 and (g+f)%4 == 0: return -delta,g,ZZ((f+g)/2)
  if delta<0 and (g-f)%4 == 0: return delta,f,ZZ((g-f)/2)
  if delta<0 and (g+f)%4 == 0: return delta,f,ZZ((g+f)/2)
  return 1+delta,f,ZZ(g/2)

def condswap(delta,f,g):
  if delta>0 and g&1: return -delta,g,-f
  return delta,f,g

def elimination(delta,f,g):
  return 1+delta,f,ZZ((g+(g&1)*f)/2)

def T(delta,f,g):
  M2Q = MatrixSpace(QQ,2)
  if delta>0 and g&1:
    return M2Q((0,1,-1/2,1/2))
  return M2Q((1,0,(g&1)/2,1/2))

def S(delta,f,g):
  M2ZZ = MatrixSpace(ZZ,2)
  if delta>0 and g&1:
    return M2ZZ((1,0,1,-1))
  return M2ZZ((1,0,1,1))

delta,f,g = 1,1,1
delta,f,g = posdivstep(delta,f,g)
delta,f,g = posdivstep(delta,f,g)
assert (delta,f,g) == (1,1,1)
coverage['posdivstep'] += 1

for b in range(2,20):
  delta,f,g = 1,1,3*2^(b-2)
  iterations = 0
  while g != 0:
    delta,f,g = pmdivstep(delta,f,g)
    iterations += 1
    if iterations <= b-2:
      assert (delta,f,g) == (iterations+1,1,3*2^(b-2-iterations))
      coverage['pmdivstep1'] += 1
    elif iterations == 3*b-1:
      assert (delta,f,g) == (-1,1,0)
      coverage['pmdivstep2'] += 1
    else:
      assert delta == 2-b+(iterations-b)//2
      assert f == 3 if iterations<3*b-3 else 1
      assert g == 2 if (iterations-b)&1 else 1
      coverage['pmdivstep3'] += 1
  assert iterations == 3*b-1
  coverage['pmdivstep'] += 1

for loop in range(500000):
  bits1 = randrange(300)
  bits2 = randrange(300)
  d = randrange(2**bits1)
  if not d&1: continue
  f = d*randrange(2**bits2)
  if randrange(2): f = -f
  if not f&1: continue
  g = d*randrange(2**bits2)
  if randrange(2): g = -g

  delta = randrange(-9,10)
  f,g = ZZ(f),ZZ(g)

  delta1,f1,g1 = divstep(delta,f,g)

  assert vector((f1,g1)) == T(delta,f,g)*vector((f,g))
  coverage['T'] += 1

  assert vector((1,delta1)) == S(delta,f,g)*vector((1,delta))
  coverage['S'] += 1

  deltac,fc,gc = condswap(delta,f,g)
  assert elimination(deltac,fc,gc) == (delta1,f1,g1)
  coverage['decomposition'] += 1

  assert max(abs(f1),abs(g1)) <= max(abs(f),abs(g))
  coverage['range'] += 1

  n = max(f.nbits()-1,g.nbits()-1)
  while True:
    if -2^n <= f and f < 2^n and -2^n <= g and g < 2^n: break
    n += 1
  assert -2^n <= f1 and f1 < 2^n and -2^n <= g1 and g1 < 2^n
  coverage['range2'] += 1

  n = max(f.nbits()-1,g.nbits()-1)
  while True:
    if -2^n < f and f < 2^n and -2^n < g and g < 2^n: break
    n += 1
  assert -2^n < f1 and f1 < 2^n and -2^n < g1 and g1 < 2^n
  coverage['range3'] += 1

  delta3,f3,g3 = -2,f,g
  delta3,f3,g3 = divstep(delta3,f3,g3)
  delta3,f3,g3 = divstep(delta3,f3,g3)
  delta3,f3,g3 = divstep(delta3,f3,g3)
  assert delta3 == 1
  assert f3 == f
  assert 8*g3 in [g,g+f,g+2*f,g+3*f,g+4*f,g+5*f,g+6*f,g+7*f]
  coverage['3divstep'] += 1

  delta3,f3,g3 = -2,f,g
  delta3,f3,g3 = cdivstep(delta3,f3,g3)
  delta3,f3,g3 = cdivstep(delta3,f3,g3)
  delta3,f3,g3 = cdivstep(delta3,f3,g3)
  assert delta3 == 1
  assert f3 == f
  assert 8*g3 in [g-3*f,g-2*f,g-f,g,g+f,g+2*f,g+3*f,g+4*f]
  coverage['3cdivstep'] += 1

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