import sys

print('version 2020.03.17')

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

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

  p = min(p for p,e in factor(n))
  for l in range(p):
    assert gcd(n,factorialmod(l,n)) == 1
  for l in range(p+1,n+1):
    assert gcd(n,factorialmod(l,n)) > 1