An enormous factorial

Author: Shlomi Fish

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);
}