import random

tests = [
  'STprod',
  'STprod1',
  'truncate',
]
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 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(10000):
  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)
  M2Z = MatrixSpace(ZZ,2)

  for n in range(-10,10):
    for m in range(-10,10):
      if n >= m and m >= 0:
        Tprod = M2Q(prod(Tn[i] for i in reversed(range(m,n))))
        Sprod = M2Z(prod(Sn[i] for i in reversed(range(m,n))))
        assert vector((fn[n],gn[n])) == Tprod * vector((fn[m],gn[m]))
        assert vector((1,deltan[n])) == Sprod * vector((1,deltan[m]))
        if n == m:
          assert Tprod == 1
          assert Sprod == 1
          coverage['STprod1'] += 1
        if n > m:
          ZZ(2^(n-m-1)*Tprod[0][0])
          ZZ(2^(n-m-1)*Tprod[0][1])
          ZZ(2^(n-m)*Tprod[1][0])
          ZZ(2^(n-m)*Tprod[1][1])
          assert Sprod[0][0] == 1
          assert Sprod[0][1] == 0
          assert Sprod[1][0] >= 2-(n-m)
          assert Sprod[1][0] <= n-m
          assert (Sprod[1][0] - (n-m)) % 2 == 0
          assert Sprod[1][1] in [1,-1]
          coverage['STprod'] += 1

  primedeltan = {}
  primefn = {}
  primegn = {}
  primeTn = {}
  primeSn = {}

  for m in range(20):
    t = randrange(20)
    delta = deltan[m]
    f = fn[m] % 2^t
    g = gn[m] % 2^t
    for n in range(m,20):
      primedeltan[n] = delta
      primefn[n] = f
      primegn[n] = g
      primeTn[n] = T(delta,f,g)
      primeSn[n] = S(delta,f,g)

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

    for n in range(20):
      if m < n and n <= m + t:
        assert primeTn[n-1] == Tn[n-1]
        assert primeSn[n-1] == Sn[n-1]
        assert primedeltan[n] == deltan[n]
        assert primefn[n] % 2^(t-(n-m-1)) == fn[n] % 2^(t-(n-m-1))
        assert primegn[n] % 2^(t-(n-m)) == gn[n] % 2^(t-(n-m))
        coverage['truncate'] += 1

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