import random

tests = [
  'truncate0',
  'truncate',
]
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 jumpdivsteps2(n,t,delta,f,g):
  assert t >= n and n >= 0
  if n <= 1: return divsteps2(n,t,delta,f,g)

  # j = n//2
  # "j can be taken anywhere between 1 and n-1"
  j = randrange(1,n)

  delta,f1,g1,P1 = jumpdivsteps2(j,j,delta,f,g)
  f,g = P1*vector((f,g))
  f,g = truncate(ZZ(f),t-j),truncate(ZZ(g),t-j)

  delta,f2,g2,P2 = jumpdivsteps2(n-j,n-j,delta,f,g)
  f,g = P2*vector((f,g))
  f,g = truncate(ZZ(f),t-n+1),truncate(ZZ(g),t-n)

  return delta,f,g,P2*P1

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 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))

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

  deltan = {}
  fn = {}
  gn = {}
  Tn = {}
  Sn = {}

  for n in range(20):
    deltan[n] = delta
    fn[n] = f
    gn[n] = g
    Tn[n] = T(delta,f,g)
    Sn[n] = S(delta,f,g)

    delta,f,g = divstep(delta,f,g)

  M2Q = MatrixSpace(QQ,2)

  for n in range(20):
    t = randrange(40)
    if n <= t:
      delta,f,g,P = divsteps2(n,t,deltan[0],fn[0],gn[0])
      assert delta == deltan[n]
      if n == 0:
        assert f == truncate(fn[n],t)
        coverage['truncate0'] += 1
      else:
        assert f == truncate(fn[n],t-(n-1))
        coverage['truncate'] += 1
      assert g == truncate(gn[n],t-n)
      assert P == M2Q(prod(Tn[i] for i in reversed(range(0,n))))
      assert (delta,f,g,P) == jumpdivsteps2(n,t,deltan[0],fn[0],gn[0])

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