(PARI) lpf(n)=my(f=factor(n)[, 1]); f[1]
factorial_lval(n, p)={
my(s);
while(n\=p, s+=n);
s
};
isfactorial(n)={
if(n<3, return(n>0));
my(v2=valuation(n, 2), mn=v2+1, mx=mn+log(v2+1.5)\log(2), t, c);
while (mx - mn > 1,
t = mn + (mx - mn)\2;
c = factorial_lval(t, 2);
if (c < v2,
mn = t+1
,
if (c > v2,
mx = t-1
,
mx = bitor(t, 1);
mn = max(mn, mx-1)
)
)
);
if (mn < mx,
my(p=lpf(mx));
t = valuation(n, p);
c = factorial_lval(mx, p);
if (t > c, return(0));
if (t == c,
mn = mx
)
);
n==mn!
};
is(n)=isfactorial(sumdigits(n!))