import random

tests = [
  'gcd2',
  'recip2',
  'divsteps',
]
coverage = {t:0 for t in tests}

def truncate(f,t):
  if t == 0: return 0
  twot = 1<<(t-1)
  return ((f+twot)&(2*twot-1))-twot

def divsteps2(n,t,delta,f,g):
  assert t >= n and n >= 0
  f,g = truncate(f,t),truncate(g,t)
  u,v,q,r = 1,0,0,1

  while n > 0:
    f = truncate(f,t)
    if delta > 0 and g&1: delta,f,g,u,v,q,r = -delta,g,-f,q,r,-u,-v
    g0 = g&1
    delta,g,q,r = 1+delta,(g+g0*f)/2,(q+g0*u)/2,(r+g0*v)/2
    n,t = n-1,t-1
    g = truncate(ZZ(g),t)

  M2Q = MatrixSpace(QQ,2)
  return delta,f,g,M2Q((u,v,q,r))

def iterations(d):
  return (49*d+80)//17 if d<46 else (49*d+57)//17

for d in range(1000):
  assert iterations(d) <= 4+3*d

def gcd2(f,g):
  assert f & 1
  d = max(f.nbits(),g.nbits())
  m = iterations(d)
  delta,fm,gm,P = divsteps2(m,m+d,1,f,g)
  assert gm == 0
  return abs(fm)

def recip2(f,g):
  assert f & 1
  d = max(f.nbits(),g.nbits())
  m = iterations(d)
  precomp = Integers(f)((f+1)/2)^(m-1)
  delta,fm,gm,P = divsteps2(m,m+1,1,f,g)
  assert gm == 0
  V = sign(fm)*ZZ(P[0][1]*2^(m-1))
  return ZZ(V*precomp)

for loop in range(20000):
  if randrange(2):
    bits1 = 1
    d = 1
  else:
    bits1 = randrange(150)
    d = randrange(2**bits1)
    if not d&1: continue
  bits2 = randrange(150)
  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)

  assert gcd2(f,g) == gcd(f,g)
  coverage['gcd2'] += 1

  if gcd(f,g) == 1 and f > 1:
    assert (recip2(f,g)*g)%abs(f) == 1
    coverage['recip2'] += 1

  d = 0
  while f^2 + 4*g^2 > 5*4^d:
    d += RR.random_element(0,1)

  m = floor((49*d+80)/17) if d<46 else floor((49*d+57)/17)
  m = ZZ(m)
  assert m >= 1
  m += randrange(3)
  assert m >= 1

  deltam,fm,gm,P = divsteps2(m,ZZ(m+bits1+bits2+100),1,f,g)
  assert gm == 0
  assert fm in [gcd(f,g),-gcd(f,g)]
  v = P[0][1]
  ZZ(2^(m-1)*v)
  assert (ZZ(2^(m-1)*v)*g)%f == (2^(m-1)*fm)%f
  coverage['divsteps'] += 1

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