h:= proc(l) local n; n:=nops(l); add(i, i=l)!/mul(mul(1+l[i]-j
+add(`if`(l[k]>=j, 1, 0), k=i+1..n), j=1..l[i]), i=1..n)
end:
g:= proc(n, i, l)
`if`(n=0 or i=1, h([l[], 1$n])^2, `if`(i<1, 0,
add(g(n-i*j, i-1, [l[], i$j]), j=0..n/i)))
end:
a:= n-> n! -g(n, 3, []):
# Alternative:
a:= proc(n) option remember; `if`(n<3, 0, `if`(n=4, 1,
((13-11*n-40*n^2+10*n^3+n^4)*a(n-1) -(10*n^2-9*n-31)*(n-1)^2*a(n-2)
+9*(n-1)^2*(n-2)^2*a(n-3)) / ((n-4)*(n+2)^2)))
end:
