tests = [
  'firstmatrices',
  '328',
  '328radius',
  '328wradius',
  '60321097',
  '22293',
  'F',
  'pos34',
  'pos14',
  'abseigenvalues',
  'notabseigenvalues',
  'PG',
  '5467345',
  '1287513',
  '2718281',
  'matnorm0',
  'matnorm',
  'intvecnorm',
  'scalevecnorm',
  'scalematnorm',
  'matvecnorm',
  'matmatnorm',
  'normisatmost',
  'notnormisatmost',
  '1sqrt2',
  'beta',
  'gamma',
  'biggamma',
  'main0',
  'main',
  'Rnorm',
  'wnorm',
  '14410',
]
coverage = {t:0 for t in tests}

M2Q = MatrixSpace(QQ,2)

def M(e,q):
  return M2Q((0,1/2^e,-1/2^e,q/2^(2*e)))

def posM(e,q):
  return M2Q((0,1/2^e,1/2^e,q/2^(2*e)))

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 gcdalgo(R0,R1):
  R = [R0,R1]
  e = [R[0].ord(2)]
  P = []
  while R[-1] != 0:
    Ri,Ri1 = R[-1],R[-2]
    ei,ei1 = Ri.ord(2),Ri1.ord(2)
    qi = 2^ei*div2(ZZ(Ri1/2^ei1),Ri)
    assert ei1 == e[-1]
    e += [ei]
    P += [M(ei,qi)]
    R += [-mod2(ZZ(Ri1/2^ei1),Ri)/2^ei]

  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

  return sum(e),P

for e in range(1,8):
  for q in range(1-2^(e+1),2^(e+1),2):
    assert all(abs(v) == 1/2^e for v in M(e,q).eigenvalues())
    coverage['abseigenvalues'] += 1
  for q in [-1-2^(e+1),1+2^(e+1)]:
    assert not all(abs(v) == 1/2^e for v in M(e,q).eigenvalues())
    coverage['notabseigenvalues'] += 1

assert M(1,1) == M2Q((0,1/2,-1/2,1/4))
assert M(1,3) == M2Q((0,1/2,-1/2,3/4))
assert M(2,1) == M2Q((0,1/4,-1/4,1/16))
assert M(2,3) == M2Q((0,1/4,-1/4,3/16))
assert M(2,5) == M2Q((0,1/4,-1/4,5/16))
assert M(2,7) == M2Q((0,1/4,-1/4,7/16))
coverage['firstmatrices'] += 1

P = (1/4^7)*M2Q((328,-380,-158,233))
assert P == M(2,1)*M(1,3)^2*M(1,1)*M(1,3)^2
L = [M(2,1),M(1,3),M(1,3),M(1,1),M(1,3),M(1,3)]
assert P == prod(L)
coverage['328'] += 1

r = min((P^3).eigenvalues())
assert r > 1/2^27.2
assert r < 1/2^27.1
assert r > 1/2^(21*1.30)
assert r < 1/2^(21*1.29)
r = max((P^3).eigenvalues())
assert r > 1/2^14.9
assert r < 1/2^14.8
assert r > 1/2^(21*0.708)
assert r < 1/2^(21*0.707)
assert 133569/2^31 > 2^(-13.973)
assert 133569/2^31 < 2^(-13.972)

assert P^(-3) == M2Q((60321097,113368060,47137246,88663112))
r = max(abs(v) for v in (P^(-3)).eigenvalues())
assert r > 2^27.1
assert r < 2^27.2

w,P = gcdalgo(60321097,47137246)
assert w == 21
assert len(P) == 18
assert P == list(reversed(L+L+L))
coverage['60321097'] += 1

w,P = gcdalgo(22293,-128330)
assert w == 23
assert len(P) == 19
coverage['22293'] += 1

F = M2Q((0,1,1,-1))
r = max(F.eigenvalues())
assert r == (sqrt(5)-1)/2
assert r > 1/2^0.70
assert r < 1/2^0.69
r = min(F.eigenvalues())
assert r == -(sqrt(5)+1)/2
assert r > -2^0.70
assert r < -2^0.69
assert abs(r) > 1
r = max(abs(v) for v in (F^(-1)).eigenvalues())
assert r == (sqrt(5)+1)/2
assert r < 2^0.70
assert r > 2^0.69
coverage['F'] += 1

for rotation in range(6):
  P = prod(L[rotation:]+L[:rotation])
  assert max(abs(v) for v in P.eigenvalues()) == (561+sqrt(249185))/2^15
  assert (561+sqrt(249185))/2^15 > 0.03235425826
  assert (561+sqrt(249185))/2^15 < 0.03235425827
  assert (561+sqrt(249185))/2^15 > 0.61253803919^7
  assert (561+sqrt(249185))/2^15 < 0.61253803920^7
  assert (561+sqrt(249185))/2^15 > 1/2^4.949900587
  assert (561+sqrt(249185))/2^15 < 1/2^4.949900586
  assert (561+sqrt(249185))/2^15 > 1/2^(7*0.7071286553)
  assert (561+sqrt(249185))/2^15 < 1/2^(7*0.7071286552)
  for w in range(7,35,7):
    w7 = ZZ(w/7)
    r = max(abs(v) for v in (P^w7).eigenvalues())
    assert r == ((561+sqrt(249185))/2^15)^w7
    assert r > 1/2^(w*0.7071286553)
    assert r < 1/2^(w*0.7071286552)
    coverage['328wradius'] += 1
  assert 7/(15-log(561+sqrt(249185),2)) > 1.414169815
  assert 7/(15-log(561+sqrt(249185),2)) < 1.414169816
  assert 14/(15-log(561+sqrt(249185),2)) > 2.828339631
  assert 14/(15-log(561+sqrt(249185),2)) < 2.828339632
  coverage['328radius'] += 1

P = posM(1,3)
assert P == M2Q((0,1/2,1/2,3/4))
assert max(abs(v) for v in P.eigenvalues()) == 1
assert vector((1,2)) == P * vector((1,2))
coverage['pos34'] += 1

P = posM(1,1)
assert P == M2Q((0,1/2,1/2,1/4))
r = max(abs(v) for v in P.eigenvalues())
assert r == (1+sqrt(17))/8
assert r > 0.6403882032
assert r < 0.6403882033
assert 2/(log(sqrt(17)-1,2)-1) > 3.110510062
assert 2/(log(sqrt(17)-1,2)-1) < 3.110510063
r = min(P.eigenvalues())
assert r == -2/(sqrt(17)+1)
assert r > -0.39039
assert r < -0.39038
assert log(2)/log((sqrt(17)+1)/2) > 0.73690
assert log(2)/log((sqrt(17)+1)/2) < 0.73691
assert 2*log(2)/log((sqrt(17)+1)/2) > 1.4738
assert 2*log(2)/log((sqrt(17)+1)/2) < 1.4739
coverage['pos14'] += 1

def G(n):
  assert n >= 0
  if n == 0: return 0
  if n == 1: return 1
  return -G(n-1)+4*G(n-2)

def szsteps(R0,R1):
  result = 0
  a,b = R0,R1
  while b != 0:
    result += 1
    q = ZZ(Integers(2^(b.ord(2)-a.ord(2)+1))(-a/2^a.ord(2))/ZZ(b/2^
b.ord(2)))
    if q > 2^(b.ord(2)-a.ord(2)): q -= 2^(1+b.ord(2)-a.ord(2))
    assert q&1
    assert abs(q) < 2^(b.ord(2)-a.ord(2))
    r = ZZ(a+q*b/2^(b.ord(2)-a.ord(2)))
    assert r.ord(2) >= b.ord(2)+1
    a,b = b,r

  G = gcd(R0,R1)
  assert abs(a)/2^a.ord(2) == G/2^G.ord(2)
  return result

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 cdivsteps(delta,f,g):
  result = 0
  while g != 0:
    result += 1
    delta,f,g = cdivstep(delta,f,g)
  return result

for n in range(10):
  assert P^(-(n+1)) == M2Q((G(n+2),2*G(n+1),2*G(n+1),4*G(n)))
  coverage['PG'] += 1

assert G(17) == 2135149
assert G(18) == -5467345
assert szsteps(-5467345,2*2135149) == 17
assert cdivsteps(1,-5467345,2135149) == 34
coverage['5467345'] += 1

assert szsteps(-1287513,2*622123) == 24
assert cdivsteps(1,-1287513,622123) == 58
coverage['1287513'] += 1

assert cdivsteps(1,-2718281,3141592) == 45
coverage['2718281'] += 1

def vecnorm(v):
  assert len(v) == 2
  return sqrt(v[0]^2+v[1]^2)

def matnorm(P):
  (a,b),(c,d) = P.transpose() * P
  return sqrt((a+d+sqrt((a-d)^2+4*b^2))/2)

assert matnorm(M2Q(0)) == 0
coverage['matnorm0'] += 1

for loop in range(100):
  P = M2Q([randrange(-100,100) for i in range(4)])
  Pnorm = matnorm(P)

  def f(x):
    return vecnorm(P * vector((cos(x),sin(x))))
  m,x = find_local_maximum(f,0,3.2)
  assert m >= 0.999999999 * Pnorm
  assert Pnorm >= 0.999999999 * m
  coverage['matnorm'] += 1

  v = vector((randrange(-5,5),randrange(-5,5)))
  if v != 0:
    assert vecnorm(v) >= 1
    coverage['intvecnorm'] += 1

  v = vector((randrange(-100,100),randrange(-100,100)))
  r = randrange(-100,100)
  assert vecnorm(r*v) == abs(r) * vecnorm(v)
  coverage['scalevecnorm'] += 1

  assert matnorm(r*P) == abs(r) * Pnorm
  coverage['scalematnorm'] += 1

  assert vecnorm(P*v) <= Pnorm * vecnorm(v)
  coverage['matvecnorm'] += 1

  Q = M2Q([randrange(-100,100) for i in range(4)])
  assert matnorm(P*Q) <= Pnorm * matnorm(Q)
  coverage['matmatnorm'] += 1

  def normisatmost(P,N):
    (a,b),(c,d) = P.transpose() * P
    return N>=0 and 2*N^2 >= a+d and (2*N^2-a-d)^2 >= (a-d)^2 + 4*b^2

  for N in [Pnorm,(1000/999)*Pnorm,(999/1000)*Pnorm,
            0,1/1000,-1/1000,
            -Pnorm,-(1000/999)*Pnorm,-(999/1000)*Pnorm]:
    if matnorm(P) <= N:
      assert normisatmost(P,N)
      coverage['normisatmost'] += 1
    else:
      assert not normisatmost(P,N)
      coverage['notnormisatmost'] += 1

for e in range(1,8):
  for q in range(1-2^(e+1),2^(e+1),2):
    assert matnorm(M(e,q)) < (1+sqrt(2))/2^e
    coverage['1sqrt2'] += 1

earlybounds = { 0:1, 1:1, 2:689491/2^20, 3:779411/2^21,
  4:880833/2^22, 5:165219/2^20, 6:97723/2^20, 7:882313/2^24,
  8:306733/2^23, 9:92045/2^22, 10:439213/2^25, 11:281681/2^25,
  12:689007/2^27, 13:824303/2^28, 14:257817/2^27, 15:634229/2^29,
  16:386245/2^29, 17:942951/2^31, 18:583433/2^31, 19:713653/2^32,
  20:432891/2^32, 21:133569/2^31, 22:328293/2^33, 23:800421/2^35,
  24:489233/2^35, 25:604059/2^36, 26:738889/2^37, 27:112215/2^35,
  28:276775/2^37, 29:84973/2^36, 30:829297/2^40, 31:253443/2^39,
  32:625405/2^41, 33:95625/2^39, 34:465055/2^42, 35:286567/2^42,
  36:175951/2^42, 37:858637/2^45, 38:65647/2^42, 39:40469/2^42,
  40:24751/2^42, 41:240917/2^46, 42:593411/2^48, 43:364337/2^48,
  44:889015/2^50, 45:543791/2^50, 46:41899/2^47, 47:205005/2^50,
  48:997791/2^53, 49:307191/2^52, 50:754423/2^54, 51:57527/2^51,
  52:281515/2^54, 53:694073/2^56, 54:212249/2^55, 55:258273/2^56,
  56:636093/2^58, 57:781081/2^59, 58:952959/2^60, 59:291475/2^59,
  60:718599/2^61, 61:878997/2^62, 62:534821/2^62, 63:329285/2^62,
  64:404341/2^63, 65:986633/2^65, 66:603553/2^65,
}

def alpha(w):
  if w >= 67: return (633/1024)^w
  return earlybounds[w]

assert all(alpha(w)^49<2^(-(34*w-23)) for w in range(31,100))
assert min((633/1024)^w/alpha(w) for w in range(68))==633^5/(2^30*165219)

betacache = {}
def beta(w):
  if w in betacache: return betacache[w]
  if w >= 67:
    result = (633/1024)^w*633^5/(2^30*165219)
  else:
    result = min(alpha(w+j)/alpha(j) for j in range(68))
  betacache[w] = result
  return result

for w in range(100):
  assert beta(w) == min(alpha(w+j)/alpha(j) for j in range(100))
  coverage['beta'] += 1

gammacache = {}
def gamma(w,e):
  assert w >= 0
  assert e >= 1
  if (w,e) in gammacache: return gammacache[w,e]
  result = min(beta(w+j)*2^j*70/169 for j in range(e,e+68))
  gammacache[w,e] = result
  return result

for loop in range(1000):
  w = randrange(100)
  e = randrange(1,100)
  assert gamma(w,e) == min(beta(w+j)*2^j*70/169 for j in range(e,e+100))
  coverage['gamma'] += 1
  if w + e >= 67:
    assert gamma(w,e) == 2^e*(633/1024)^(w+e)*(70/169)*633^5/(2^30*165219)
    coverage['biggamma'] += 1

assert matnorm(M2Q(1)) <= alpha(0)
coverage['main0'] += 1

for loop in range(1000):
  j = randrange(1,10)
  e = [max(1,randrange(4)) for i in range(j)]
  q = [randrange(1,2^(e[i]+1),2) for i in range(j)]
  assert matnorm(prod(M(e[i],q[i]) for i in range(j))) <= alpha(sum(e))
  coverage['main'] += 1

for loop in range(100):
  bits = randrange(100)
  R0 = 1+2*randrange(-2^bits,2^bits)
  R1 = 2*randrange(-2^bits,2^bits)
  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

  startnorm = vecnorm((R0,R1))
  for j in range(t):
    w = sum(e[1:j+1])
    assert vecnorm((R[j]/2^e[j],R[j+1])) <= alpha(w) * startnorm
    coverage['Rnorm'] += 1
    if w >= 67:
      assert w <= log(startnorm,1024/633)
      coverage['wnorm'] += 1

  b = log(max(abs(R0),abs(R1)),2)
  assert abs(R0) <= 2^b
  assert abs(R1) <= 2^b
  assert startnorm <= 2^(b+1/2)
  assert log(startnorm,1024/633) <= (b+1/2)/log(1024/633,2)
  assert 1/log(1024/633,2) > 1.4410
  assert 1/log(1024/633,2) < 1.4411
  coverage['14410'] += 1

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