# An enormous factorial

Author: Shlomi Fish

http://projecteuler.net/problem=288

For any prime `p` the number `N(p,q)` is defined by `N(p,q) = ∑_n=0 to q` `T_n*p_n` with `T_n` generated by the following random number generator:

```S0 = 290797
Sn+1 = Sn2 mod 50515093
Tn = Sn mod p
```

Let `Nfac(p,q)` be the factorial of `N(p,q)`. Let `NF(p,q)` be the number of factors `p` in `Nfac(p,q)`.

You are given that `NF(3,10000) mod 3^20=624955285`.

Find `NF(61,10^7) mod 61^10`.

Source code: prob288-shlomif.p6

```use v6;

sub factorial_factor_exp(\$n , \$f)
{
if (\$n < \$f)
{
return 0;
}
else
{
my \$div = \$n / \$f;
return \$div + factorial_factor_exp(\$div, \$f);
}
}

sub MAIN(:\$verbose = False) {
my int @t_n;

my \$N_LIM = 10;
my int \$BASE = 61;
my \$LIM = 10_000_000;

my \$S_0 = 290797;
my int \$s = \$S_0;

for (0 .. \$N_LIM-1) -> \$n {
@t_n.push(\$s % \$BASE);
\$s = ((\$s * \$s) % 50515093);
}

my int \$sum = 0;
for (\$N_LIM .. \$LIM) -> int \$n {
if \$n % 10_000 == 0 {
say "Reached \$n" if \$verbose;
}
\$sum += (\$s % \$BASE);
\$s = ((\$s * \$s) % 50515093);
}

sub f(\$n)
{
return factorial_factor_exp(\$n, (\$BASE)) % ((\$BASE) ** \$N_LIM);
}

say "Solution == ", +([+] flat(
(map { f((\$BASE) ** \$_) * @t_n[\$_] }, 1 .. @t_n-1),
\$sum * f((\$BASE) ** \$N_LIM)
)) % ((\$BASE) ** \$N_LIM);
}
```