!YS-v0

# Use *' (vs. *) to allow arbitrary length arithmetic:
mul =: value("*'")

# Number of random times to test each n:
TIMES =: 1000000

defn main(max=20):
  say: "\tAvg\tExp\tDiff"

  doseq n (1 .. max):
    anal =: analytical(n):F
    avg  =: avg-cycle-length(n TIMES):F
    diff =: abs(100 * (1 - (avg / anal)))
    say:
      format "%3d\t%.4f\t%.4f\t%.2f%%": n avg anal diff

defn factorial(n):
  :: n!
  range(1 n.++): .mul(*)

defn analytical(n):
  :: Analytical computation
  1 .. n:
    .map(\(factorial(n) / pow(n _) / factorial(n - _)))
    .sum()

defn single-test-cycle-length(n):
  :: Single random test of cycle length
  loop count 0, bits 0, x 1:
    if bit-and(x bits) == 0:
      recur: count.++ bit-or(bits x) bit-shift-left(1 rand-int(n))
      else: count

defn avg-cycle-length(n times):
  :: Average results of single tests of cycle lengths
  div _ times:
    sum:
      for i (range times):
        single-test-cycle-length: n