fac(n) computes factorial, facmod(n, m) is faster in reducing all recursive results "mod m"

fac.bc

define fac(n) {
  if (n==0)  return 1;
  return n*fac(n-1)
}

define facmod(n, m) {
  if (n==0)  return 1;
  return n*fac(n-1)%m
}

Wilson's theorem:
(n-1)! ≡ -1 (mod n) ⇔ n is prime
https://en.wikipedia.org/wiki/Wilson%27s_theorem

pi@pi400-64:~ $ time (echo "facmod(9972, 9973)" | bc -q fac.bc )
9972

real	0m6.119s
user	0m5.550s
sys	0m0.050s
pi@pi400-64:~ $
pi@pi400-64:~ $ bc -q fac.bc
fac(0)
1
fac(10)
3628800

Only for n starting at 30,000 "fac(n)%m" becomes slower than "facmod(n, m)":

pi@pi400-64:~ $ time (echo "fac(39988)%39989" | bc -q fac.bc )
39988

real	0m55.654s
user	0m54.578s
sys	0m0.959s
pi@pi400-64:~ $ time (echo "facmod(39988, 39989)" | bc -q fac.bc )
39988

real	0m52.071s
user	0m51.061s
sys	0m0.897s
pi@pi400-64:~ $