tests = [
  'divstep2e',
  'delta',
  'divstep2w',
  'gcd',
  '4917',
  'W',
]
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)

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 manydivstep(n,delta,f,g):
  while n>0:
    delta,f,g = divstep(delta,f,g)
    n -= 1
  return delta,f,g

for loop in range(10000):
  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)

  q = ZZ(2^e*div2(f,g))
  r = mod2(f,g)

  assert q&1
  assert q >= 1
  assert q <= 2^(e+1)-1
  assert r == f-q*ZZ(g/2^e)
  assert r.ord(2) >= e+1

  delta,fn,gn = 1,f,ZZ(g/2)
  for n in range(2*e):
    delta,fn,gn = divstep(delta,fn,gn)

  assert (delta,fn,gn) == manydivstep(2*e,1,f,ZZ(g/2))

  assert delta == 1
  assert fn == g/2^e
  assert gn == (q*g/2^e-f)/2^(e+1)
  assert gn == -r/2^(e+1)
  coverage['divstep2e'] += 1

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

      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)
      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]
        if R[i] == 0:
          assert i == len(R)-1

      t = len(R)-2
      assert R[t+1] == 0
      assert abs(R[t]/2^e[t]) == G

      w = sum(e[1:t+1])
      for n in [2*w,2*w+1,2*w+2]:
        assert manydivstep(n,1,R0,ZZ(R1/2)) == (1+n-2*w,R[t]/2^e[t],0)
        coverage['delta'] += 1

      for j in range(t):
        w = sum(e[1:j+1])
        assert manydivstep(2*w,1,R0,ZZ(R1/2)) == (1,R[j]/2^e[j],R[j+1]/2)
        coverage['divstep2w'] += 1

      delta,fn,gn = 1,R0,ZZ(R1/2)
      while gn != 0: delta,fn,gn = divstep(delta,fn,gn)
      assert fn in [G,-G]
      delta,fn,gn = divstep(delta,fn,gn)
      assert fn in [G,-G]
      delta,fn,gn = divstep(delta,fn,gn)
      assert fn in [G,-G]
      coverage['gcd'] += 1

      b = log(sqrt(R0^2+R1^2),2)
      if b <= 21: n = ZZ(floor(19*b/7))
      if 21 < b and b <= 46: n = floor((49*b+23)/17)
      if b > 46: n = floor(49*b/17)
      delta,fn,gn = manydivstep(n,1,R0,ZZ(R1/2))
      assert gn == 0
      assert fn in [G,-G]
      assert n <= floor((49*b+23)/17)
      assert 19/7 < 49/17
      coverage['4917'] += 1

def steps(R0,R1):
  delta,f,g = 1,R0,ZZ(R1/2)
  steps = 0
  while g != 0:
    delta,f,g = divstep(delta,f,g)
    steps += 1
  return steps

def W(s):
  H = 1
  result = None
  while not result:
    for R0 in range(-H,H+1):
      if R0&1:
        for R1 in range(-H,H+1):
          if R1&1: continue
          if R0^2+R1^2 > H^2: continue
          if steps(R0,R1) >= s:
            if not result: result = R0^2+R1^2
            if result > R0^2+R1^2: result = R0^2+R1^2
    H *= 2
  return result

for s in range(11):
  Ws = W(s)
  assert Ws == [1,5,5,5,17,17,65,65,157,181,421][s]
  R0,R1 = [
    (1,0),(1,-2),(1,-2),(1,-2),(1,-4),(1,-4),(1,-8),(1,-8),(11,-6),(9,-10),(15,-14)
  ][s]
  assert steps(R0,R1) >= s
  assert R0^2+R1^2 == Ws
coverage['W'] += 1

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