import sys
import random

print('version 2020.03.17')

def multfactorialmod(alpha,zeta,l,n):
  Zn = Integers(n)
  ZnX.<X> = Zn[]
  b = floor(sqrt(l))
  if b == 0:
    I = []
    J = []
  else:
    I = range(1,b+1)
    J = range(0,b^2,b)
    assert min(I) == 1
    assert max(I) == b
    assert len(I) == b
    assert min(J) == 0
    assert max(J) == b^2-b
    assert len(J) == b
    assert all(j%b == 0 for j in J)
  hI = ZnX(prod(X-zeta^i for i in I))
  hJ = ZnX(prod(X-alpha/zeta^j for j in J))
  HJ = ZnX(prod(zeta^j*X-alpha for j in J))
  res = hJ.resultant(hI)
  Res = HJ.resultant(hI)
  assert res == prod(hI(alpha/zeta^j) for j in J)
  assert Res == prod(zeta^(j*b)*hI(alpha/zeta^j) for j in J)
  assert all(zeta^(j*b)*hI(alpha/zeta^j) == prod(alpha-zeta^(i+j) for i in I) for j in J)
  assert Res == prod(alpha-zeta^(i+j) for i in I for j in J)
  Res *= prod(alpha-zeta^s for s in range(b^2+1,l+1))
  assert Res == prod(alpha-zeta^s for s in range(1,l+1))
  return Res

for n in range(2,100):
  print('try','multfactorialmod',n)
  sys.stdout.flush()

  alpha = Integers(n)(randrange(n))
  while True:
    zeta = Integers(n)(randrange(n))
    if gcd(n,zeta) == 1: break
  for l in range(n+1):
    multfactorialmod(alpha,zeta,l,n)
    qfactorial = prod(sum(zeta^j for j in range(i)) for i in range(2,l+1))
    assert multfactorialmod(1,zeta,l,n) == (1-zeta)^l*qfactorial

for loop in range(100):
  q = random.choice([3,5,7,9,11,13,17,19,23,25,27,29])
  k.<fieldgen> = GF(q)
  kZ.<Z> = k[]

  while True:
    zeta = k.random_element()
    if zeta.is_unit(): break

  I = set()
  J = set()
  K = set()
  S = set()
  for ijloop in range(loop):
    if randrange(2):
      inew = randrange(-loop*loop,loop*loop+1)
      if inew not in I:
        if all(inew+j not in S for j in J):
          I.add(inew)
          for j in J:
            S.add(inew+j)
    else:
      jnew = randrange(loop*loop)
      if jnew not in J:
        if all(i+jnew not in S for i in I):
          J.add(jnew)
          for i in I:
            S.add(i+jnew)

  for ijloop in range(loop):
    snew = randrange(-loop*loop,loop*loop+1)
    if snew not in S:
      K.add(snew)
      S.add(snew)

  print('try','res',loop,q,len(S),len(I),len(J),len(K))
  sys.stdout.flush()

  alpha = k.random_element()

  hI = kZ(prod(Z-zeta^i for i in I))
  HJ = kZ(prod(alpha-zeta^j*Z for j in J))
  hK = k(prod(alpha-zeta^s for s in K))
  hS = k(prod(alpha-zeta^s for s in S))

  assert hS == hI.resultant(HJ)*hK

for n in range(1000):
  print('try','I+J',n)
  sys.stdout.flush()

  m = randrange(1000)
  r = randrange(1,1000)
  S = set(range(m,m+n*r,r))
  assert len(S) == n
  b = floor(sqrt(n))
  I = set(i*r for i in range(b))
  J = set(m+j*b*r for j in range(b))
  IJ = set(i+j for i in I for j in J)
  assert IJ == set(m+k*r for k in range(b^2))
  K = set(m+k*r for k in range(b^2,n))
  assert IJ.isdisjoint(K)
  assert S == IJ.union(K)
  assert len(I) == b
  assert len(J) == b
  assert len(K) == n-b^2
  assert max(len(I),len(J),len(K)) >= sqrt(len(S)+1/4)-1/2
  assert n-b^2 <= 2*b