A131685 - OEIS
List:=NULL: for n from 1 to 1000 do m:=1: #running till n=50 will last 2 min.
for i from 1 to numtheory[pi](n) do div:=ithprime(i): d:=1: e:=0: oldmini:=-1:mini:=0:
while oldmini<>mini do e:=e+1: #the last time consuming loop could be skipped by proving e<=floor(ln(n)/ln(div)):
d:=d*div; for x from 0 to d-1 do [seq((x &^k mod d)+k mod d, k=1..n)]:contrib[d, x]:=nops(select(has, %, 0)): od:
L:=seq(add(contrib[div^j, x mod div^j], j=1..e), x=0..div^e-1); oldmini:=mini: mini:=min(L): od:
if mini<padic[ordp](n!, div) then m:=m*div^(padic[ordp](n!, div)-mini) fi; od: print(n, m); List:=List, m: od:
[List];