import random

tests = [
  'divstepR2',
  'parta',
  'partb',
  'partc',
]
coverage = {t:0 for t in tests}

def divstep(delta,f,g):
  kx = parent(f)
  x = kx.gen()
  if delta>0 and g[0]!=0:
    return 1-delta,g,kx((g[0]*f-f[0]*g)/x)
  return 1+delta,f,kx((f[0]*g-g[0]*f)/x)

def T(delta,f,g):
  kx = parent(f)
  x = kx.gen()
  M2kx = MatrixSpace(kx.fraction_field(),2)
  if delta>0 and g[0]!=0:
    return M2kx((0,1,g[0]/x,-f[0]/x))
  return M2kx((1,0,-g[0]/x,f[0]/x))

def bezout(R0,R1):
  # sage gcd and xgcd fail to return monic for, e.g., (2*x,0)
  if R0 == 0 and R1 == 0: return 0,0,0
  if R1 == 0:
    return R0/R0.leading_coefficient(),1/R0.leading_coefficient(),0
  if R0 == 0:
    return R1/R1.leading_coefficient(),0,1/R1.leading_coefficient()
  return xgcd(R0,R1)

def sanedeg(f):
  if f == 0: return -Infinity
  return f.degree()

for loop in range(200):
  q = random.choice([2,3,5,7,11,13,17,19,23,29])
  k = GF(q)
  kx.<x> = k[]
  coeffs1 = randrange(20)
  coeffs2 = randrange(1,100)
  coeffs3 = randrange(coeffs2)
  h = kx(sum(k.random_element()*x^i for i in range(coeffs1)))
  if h[0] == 0: continue
  R0 = h*kx(sum(k.random_element()*x^i for i in range(coeffs2)))
  R1 = h*kx(sum(k.random_element()*x^i for i in range(coeffs3)))

  G,U,V = bezout(R0,R1)
  assert G == U*R0+V*R1

  d0 = sanedeg(R0)
  d1 = sanedeg(R1)
  if d0 > d1 and d1 >= 0:

    l0 = R0[d0]
    l1 = R1[d1]
    R2 = R0%R1

    delta = 1
    f = kx(x^d0*R0(1/x))
    g = kx(x^(d0-1)*R1(1/x))

    for n in range(2*d0-2*d1):
      delta,f,g = divstep(delta,f,g)

    assert delta == 1
    assert f == l0^(d0-d1-1)*x^d1*R1(1/x)
    assert g == (l0^(d0-d1-1)*l1)^(d0-d1+1)*x^(d1-1)*R2(1/x)
    coverage['divstepR2'] += 1

  # importing B.1 and B.2...

  R = [R0,R1]
  while R[-1] != 0:
    R += [R[-2] % R[-1]]
  r = len(R)-1
  assert R[r] == 0

  for i in range(1,r):
    assert R[i] != 0
    assert R[i+1] == R[i-1]%R[i]

  d = [sanedeg(Ri) for Ri in R]
  l = [Ri.leading_coefficient() for Ri in R]

  for i in range(1,r): assert d[i] > d[i+1]
  assert d[r] == -Infinity

  for i in range(1,r): assert l[i] != 0

  if l[r-1] != 0:
    assert G == R[r-1]/l[r-1]

  if l[r-1] == 0:
    assert R0 == 0
    assert R1 == 0
    assert G == 0
    assert r == 1

  U = [kx(1),kx(0)]
  V = [kx(0),kx(1)]
  for i in range(1,r):
    U += [U[-2]-(R[i-1]//R[i])*U[-1]]
    V += [V[-2]-(R[i-1]//R[i])*V[-1]]

  assert U[0] == 1
  assert V[0] == 0
  assert U[1] == 0
  assert V[1] == 1
  for i in range(1,r):
    assert U[i+1] == U[i-1] - (R[i-1]//R[i])*U[i]
    assert V[i+1] == V[i-1] - (R[i-1]//R[i])*V[i]

  for i in range(r+1):
    assert R[i] == U[i]*R0+V[i]*R1

  for i in range(r):
    assert U[i]*V[i+1]-U[i+1]*V[i] == (-1)^i

  for i in range(r):
    assert V[i+1]*R[i]-V[i]*R[i+1] == (-1)^i*R0

  if d[0] > d[1]:
    for i in range(r):
      assert sanedeg(V[i+1]) == d[0] - d[i]

  # ... and now C.2

  M2kx = MatrixSpace(kx.fraction_field(),2)

  assert d[0] == d0
  assert d[1] == d1
  if d0 > d1:
    delta = 1
    f = kx(x^d0*R0(1/x))
    g = kx(x^(d0-1)*R1(1/x))
    a = 1
    b = 1
    P = M2kx(1)

    for n in range(2*d0+10):
      if n >= 2*d[0]-2*d[r-1]:
        coverage['partc'] += 1
        assert delta == 1+n-(2*d[0]-2*d[r-1])
        assert delta > 0
        assert f == a*x^d[r-1]*R[r-1](1/x)
        assert g == 0
        assert P[0][0] == a*(1/x)^(d[0]-d[r-1])*U[r-1](1/x)
        assert P[0][1] == a*(1/x)^(d[0]-d[r-1]-1)*V[r-1](1/x)
        assert P[1][0] == b*(1/x)^(n-(d[0]-d[r-1])+1)*U[r](1/x)
        assert P[1][1] == b*(1/x)^(n-(d[0]-d[r-1]))*V[r](1/x)
      else:
        oki = 0
        for i in range(r-1):
          if 2*d[0]-2*d[i] <= n and n < 2*d[0]-d[i]-d[i+1]:
            oki += 1
            coverage['parta'] += 1
            assert delta == 1+n-(2*d[0]-2*d[i])
            assert delta > 0
            assert f == a*x^d[i]*R[i](1/x)
            assert g == b*x^(2*d[0]-d[i]-n-1)*R[i+1](1/x)
            assert P[0][0] == a*(1/x)^(d[0]-d[i])*U[i](1/x)
            assert P[0][1] == a*(1/x)^(d[0]-d[i]-1)*V[i](1/x)
            assert P[1][0] == b*(1/x)^(n-(d[0]-d[i])+1)*U[i+1](1/x)
            assert P[1][1] == b*(1/x)^(n-(d[0]-d[i]))*V[i+1](1/x)
          if 2*d[0]-d[i]-d[i+1] <= n and n < 2*d[0]-2*d[i+1]:
            oki += 1
            coverage['partb'] += 1
            assert delta == 1+n-(2*d[0]-2*d[i+1])
            assert delta <= 0
            assert f == a*x^d[i+1]*R[i+1](1/x)
            Z = R[i] % (x^(2*d[0]-2*d[i+1]-n)*R[i+1])
            assert g == b*x^(2*d[0]-d[i+1]-n-1)*Z(1/x)
            X = U[i]-(R[i]//(x^(2*d[0]-2*d[i+1]-n)*R[i+1]))*x^(2*d[0]-2*d[i+1]-n)*U[i+1]
            Y = V[i]-(R[i]//(x^(2*d[0]-2*d[i+1]-n)*R[i+1]))*x^(2*d[0]-2*d[i+1]-n)*V[i+1]
            assert P[0][0] == a*(1/x)^(d[0]-d[i+1])*U[i+1](1/x)
            assert P[0][1] == a*(1/x)^(d[0]-d[i+1]-1)*V[i+1](1/x)
            assert P[1][0] == b*(1/x)^(n-(d[0]-d[i+1])+1)*X(1/x)
            assert P[1][1] == b*(1/x)^(n-(d[0]-d[i+1]))*Y(1/x)

	assert oki == 1

      if delta>0 and g[0]:
        a,b = b,g[0]*a
      else:
        a,b = a,f[0]*b

      P = T(delta,f,g)*P
      delta,f,g = divstep(delta,f,g)

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