Problem:  For any prime $p$, the number $N(p,q)$ is defined by $N(p,q)=\sum_{n=0}^q T_np^n$ with $T_n$ generated by the following random number generator:

$S_0=290797,\quad S_{n+1}=S_n^2\text{ mod }50515093,\quad T_n=S_n\text{ 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)\text{ mod }3^{20}=624955285$. Find $NF(61,10^7)\text{ mod }61^{10}$.

### Approach

I want to make sure I understand the problem correctly. The first thing I am going to do is try to reproduce the $NF(3,10000)\text{ mod }3^{20}$ result. We need a few tools in order to do that. This may not be the fastest approach, but the sequence $S_i$ can be stored in a list and used for the generation of each $T_i$. That can be done with the following code. (We’ll first get something that we know works, and then try to make it faster.)

""" Compute and return a list of elements in the sequence S_i """ def SeqGen(p,q): S = [290797] for i in range(q): S.append(S[i]**2 % 50515093) return S

Now, how long does it take to form the sequences involved? The check solution is where $p=3$ and $q=10000$.

import time start = time.time() S = SeqGen(3,10000) elapsed = time.time() - start print "%s seconds" % elapsed
0.0155670642853 seconds

The problem we’re trying to solve is where $p=61$ and $q=10^7$. That should take a bit longer to calculate.

start = time.time() S = SeqGen(61,10**7) elapsed = time.time() - start print "%s seconds" % elapsed
4.84225988388 seconds

So, we have the sequence generated, and there’s really no way to get around that. (At least we’re storing the sequence for repeated use instead of computing each $S_i$ and $T_i$ from scratch when they are needed.) But, the problem here is that we’re considering a sum of $T_i$ that will end up being absolutely huge, and then we’re supposed to consider the factorial of that! Um, how big is that number that we’re going to be taking the factorial of?

start = time.time() p = 3 q = 10000 S = SeqGen(p,q) summand = 0 for i in range(q + 1): summand += (S[i] % p) * p**i elapsed = time.time() - start print "%s found in %s seconds" % (summand,elapsed)

That gives us the following HUGE number.

730559684346872487356308424468789381739382472243104512256910727234612146\ 700534188007434920975662394367813322448421745790130109894995629026500374\ 895011594186028875546881525383692190901857070383833601818290785177074579\ 911743377542081623320139528346723997438895367511366325656859170069501892\ 782888954053372302132596260888870616595592246062506668938053947188269010\ 400553898377000621315460037720267470765385723925294615221799419878348634\ 941697100134717918504949375867315531887209110886823403867771125638452655\ 419963863463676818960699876392271245411287465034804904299149434829574363\ 770468053954258743419442793921037036133984801241622614666871400079764805\ 900359305411895238129224243327849475948701119129980806157988796433970429\ 767434356953911596347357745083143620729117746678220096159021424980589161\ 244984215592267871406818552160832420391632566752817475768883436736038842\ 522000945336450020910807258558940388286502953874414804758436489924366403\ 947302687002094708689116314943448323020960512070613984468171097011413793\ 815513010319096054358875584978751325763116736267690167927829753762880555\ 556492702423860742161506131781311489223910177639224622324264799978209792\ 137516266242000584912257651924323790209254150744514288660959513060668035\ 973225327239076026979768488390491444263297199193699839519881450073826859\ 308373986510545986393307440612730124553757195296560232991330312047030365\ 439140122350528090824215785078244449718422816025901223160225749799081881\ 988430664995829601897673857754766692861499127643968276276623422007471468\ 172908403505783118651898682362325563123083437839559958668622113095448097\ 838918074208742911144089970606803829614386725740786000635368447203667950\ 919203899690655183541386287564659803875910034486817250700675054831856742\ 035247475944246669523244511489033077185989751177302056006622400102322784\ 989845436319501972904106278352402251267862987843212251575486159493488621\ 419004070699786286867072388160751058892527781090980214514736036820101984\ 021575154483845837502770313982784775558983857516011444428943761280375624\ 001952338205797561885398526808997109816279089774072149087781516188880702\ 608233247734329861677612419964265181378869947517853199977037079229291120\ 464625858155739642990433998513357423180305979877365154314842676895851795\ 760576022818502562105123566673123108580786247852033516496746978362884952\ 942857915868486448373175276670655594702334217465802036651469475179854123\ 576334647440171805787579774605496977007910876530142161323040167528309301\ 393422160543176919470220970372232485279534778999286622468261538067226561\ 555097826633524553445830259964898613993557192047585752617190210902952965\ 023859119496947566516318311636273280475133934551257178698148121016914217\ 601635282811274292951764677139952566037328156035799548655023949268547284\ 809803753032119590215375306603046215706177234216996024131857639663170520\ 040443209001844694801103623963784228403949731061696509229838383141346412\ 396840736097533127298186296546021026067585925153039919837286307851902378\ 608106584940719518001399972148818734755824896042240388003248509954172230\ 457875696005298711270591168072531988900269290984533265472277713834821552\ 688627849442548426730677565761857748700949921996084459144768478980887361\ 308570015798008752191153537433637864164115257325737451302878323350513806\ 294637073800220289831038499134290753821824937995156685004102694409025050\ 677506077860480650134181164350989085407457828530788590469734236728651105\ 676262006578011380829166976653790493746638376475989405224819273383529385\ 532863549938170077470214978095447081535900000723694048581893421673406131\ 109533480459713318059149162064159952801277745096099653953042257676645734\ 033409098349440753210352634374498292831878619296112463252781883559685421\ 067134390006626690231502559111513270668170841237110979180067720007739317\ 585970923377914844103905254029879477702295672447225795344081751213856674\ 211906754788178078616402870316493923729160455051745371549667665511104695\ 586653367386191427354784072444333174771530684870185759160433248337062693\ 241539671265073012920758822756618788271313509549860869792768997274766995\ 313142628609896166587869086036134423061561270679939737235411692366776742\ 931830299523846661763825091147013824676353740657480723066545922030438851\ 336813999896498288327780291320769916220990607522889614427866105922093737\ 605505823866292382841719392400184241564923478503374222476232576395511449\ 778688491299558506777311849692101013641615409477217600277465271536417044\ 126580096565238246605862418216134296755489469606953685229425485095505145\ 354340285520945045486011095376520333083833266407481003815982570630207580\ 757117499879596617711378303148935948900974043533904686493402149617885320\ 110564180214417176241073001553538422708853000908450781412877954930208499\ 833275276736919712641870420034477342646070330401801478161590958515312176\ 8998208978920939146 found in 0.629118919373 seconds

Clearly, we’re going to need to do things differently, since we can’t and don’t want to take the factorial of a number that big. So, let’s consider the mathematics behind this.

### Organizing a Solution

First, the main point of the problem is to count the number of factors of some prime $p$ in $N!$ where $N$ is some big number. We don’t actually need to compute $N!$ in order to do that. Instead, we can find the factors of $p$ in $N!$ by just looking at the factors of $p$ in $n$ for $1\le n\le N$. Moreover, we’re considering this total number of factors of $p$ modulo some integer, and so we can trim down our calculations a bit. In fact, we don’t even need to consider the sum of the $T_i$, since we know that $T_i = S_i\text{ mod }p$ must be somewhere between 0 and $p-1$. We can put this all together in a relatively efficient Python script.

### Python Solution

Here’s a first attempt at achieving the correct solution for the test case.

import time   """ Compute and return a list of elements in the sequence T_i """ def T(p,q): S = [290797] for i in range(q): S.append(S[i]**2 % 50515093) T = [(i % p) for i in S] return T   """ How many factors of p does N! have modulo mod? """ def num_p_in_factorial(p,N,mod=None): quo = N//p num = quo while quo >= p: quo //= p; num += quo if mod is not None and num > mod: num = num % mod return num   def Nfac(p,q,n): # form the sequence S Tseq = T(p,q) iter = q num_factors = 0 while iter >= 0: if Tseq[iter] == 0: iter -= 1 continue elif iter >= n: num_factors += num_p_in_factorial(p,Tseq[iter]*p**n) else: num_factors += num_p_in_factorial(p,Tseq[iter]*p**iter) iter -= 1 if num_factors >= p**n: num_factors = num_factors % p**n return num_factors   start = time.time() N = Nfac(3,10000,20)   elapsed = (time.time() - start) print "found %s in %s seconds" % (N,elapsed)

When executed, we get the following result.

found 624955285 in 0.059849023819 seconds

This would work to solve the bigger problem, so let’s see how fast an initial solution is.

start = time.time() N = Nfac(61,10**7,10)   elapsed = (time.time() - start) print "found %s in %s seconds" % (N,elapsed)

And we find the solution to the problem in around 72 seconds.

found 605857431263981935 in 72.1750190258 seconds

I’m not happy with that timing. Since we’re really only calculating the num_p_in_factorial routine for exponents of $p$ that are below $n$ and we know that each $T_i$ is below $p$, we should be able to cache some of these results to limit redundant calculations. Let’s do that.

import time   """ Compute and return a list of elements in the sequence T_i """ def T(p,q): S = [290797] for i in range(q): S.append(S[i]**2 % 50515093) T = [(i % p) for i in S] return T   """ How many factors of p does N! have modulo mod? """ def num_p_in_factorial(p,N,mod=None): quo = N//p num = quo while quo >= p: quo //= p; num += quo if mod is not None and num > mod: num = num % mod return num   def Nfac(p,q,n): # form the sequence S Tseq = T(p,q) iter = q num_factors = 0 num_factors_cached = [None] * p while iter >= 0: t = Tseq[iter] if t == 0: iter -= 1 continue elif iter >= n: if num_factors_cached[t] is None: num_factors_cached[t] = num_p_in_factorial(p,t*p**n) num_factors += num_factors_cached[t] else: num_factors += num_p_in_factorial(p,Tseq[iter]*p**iter) iter -= 1 if num_factors >= p**n: num_factors = num_factors % p**n return num_factors   start = time.time() #N = Nfac(3,10000,20) N = Nfac(61,10**7,10)   elapsed = (time.time() - start) print "found %s in %s seconds" % (N,elapsed)

Now, we get something that is a bit more reasonable.

found 605857431263981935 in 10.3537530899 seconds