# (Python)  
# Returns n's prime factors in increasing order
def factors(n):
out,k=[],2
while n>=k*k:
  while n%k==0:
    out.append(k)
    n//=k
  k+=1+(k&1)
return(out if n<2 else out+[n])
# Find k1...km such that k1!k2....km!=p, ki<ki+1, sum(ki)=k1 and having factors l
def exki(k,m,prod,l,lfac):
for k1 in range(l[-1],min(k,m)+1):
  if prod%lfac[k1]: continue
  f2=lfac[k1]
  l2=list(l)
  for i in range(len(l)-1,-1,-1):
    if f2%l2[i]==0: f2//=l2.pop(i)
  if f2>1: continue # f2 has factors not in prod
  if lfac[k1]==prod: return([k1]+[1]*(k-k1))
  s=exki(k-k1,k1,prod//lfac[k1],l2,lfac)
  if s: return([k1]+s)
# Finds smallest k such that k=sum(k1,... km)
# and k!= n x (k1!k2!...km!)
def imult(n,lfac):
if n==1: return([1])
bg=n<<1
if (lfac[-1]>bg):
   inf=2
   sup=len(lfac)-1
   while(sup-inf>1):
     mid=(inf+sup)>>1
     if lfac[mid]>bg: sup=mid
     else: inf=mid
   k0=inf+(lfac[inf]<bg)
else:
   while(lfac[-1]<bg):
    lfac.append(lfac[-1]*(lfac[-1]+1))
   k0=len(lfac)-1
if lfac[k0-1]==n: return([1]*(k0-1))
maxk=.5+(n+n+.25)**.5
if not maxk.is_integer():
  maxk=.5+(n+.25)**.5
maxk+=maxk.is_integer()
for k in range(k0,int(maxk)): # n's factors must be in k
  if (k>=len(lfac)): lfac.append(lfac[-1]*(k+1))
  if lfac[k]%n==0: #  k! so that k!/n integer
   prod=lfac[k]//n  # Prod=k1!k2!...kp! w p<k
   res=exki(k,k,prod,factors(prod),lfac)
   if res: return(res)
return([n-1,1])
# Print the sequence
order=1
nseq=10011 # Max number
seedf=[1,1]+[0]*(nseq*2-1)
for p in range(2,nseq*2+1): seedf[p]=seedf[p-1]*p
for test in range(1,nseq):
    r=imult(test,seedf)
    if r[0]<test-1:
     print("a(%d) = %d"%(order,test))
     order+=1