Testar skriva lösning till ett symboltungt matte/programmerings-problem. Det är från Project Euler. Jag tänker inte skriva vilket, för att PE inte tycker om att lösningar finns på internet. Jag ändrar lite på konstanterna för att göra det mindre sökbart. Det går fortfarande att hitta vilket med lite ansträngning. Jag har fastnat på att skriva klart om Peter Lundgren och klimatet och gör något som känns mindre träligt istället.
Nuförtiden brukar jag lösa några enstaka problem på somrarna, tidigare har jag försökt hänga med och lösa de senaste när de kommer:
Definition
Vi definierar \(B(n)\) som \[B(n) = \prod_{k=0}^n \binom{n}{k}\] och \(D(n)\) som summan av alla delare till \(B(n)\). Hitta summan av alla \(D(n)\) för \(n\) mellan \(1\) och 10000. Den är lite för stor för att mata in i ett textfält, så räkna ut dess rest efter division med ett stort primtal. Säg \(P=10000000207\).
Vi ‘förenklar’ en smula (det ser inte mycket enklare ut efteråt heller): \[B(n) = \prod_{k=0}^n\frac{n!}{k!(n-k)!} = \frac{n!^{n+1}}{\prod_{k=0}^n(k!)^2}\]
Att bara använda definitionen rakt av går inte, för att talet blir för stort.
Hur stort blir B(n)?
Jag uppskattar hur stort. Det behövs inte egentligen, men det blir mycket formler ☺. Med Stirling-approximation blir \(\log B(n)\) i storleksordningen \[\log B(n) = (n+1)\cdot(n \log n - n + O(\log n))\] \[-2\sum_{k=0}^n(k \cdot \log k - k) + O\left(\sum_{k=1}^n\log k\right)\] Summan \(2\sum_{k=0}^n(k\cdot\log k)\) uppskattar jag med en integral. Den är (nästan helt säker; den är inte monotont växande mellan 0 och 1, men det är säkert sant ändå): \[2\sum_{k=0}^n(k \cdot\log k) < 2n\log n + 2 \int_0^{n} x\log x dx = 2n\log n + \left[x^2\log x - \frac{x^2}{2}\right]_{0}^{n}\] \[= 2n\log n + n^2\log n - n^2/2\] Vi vet att \(O\left(\sum_{0}^n \log k\right) = O(n \log n)\), så hela uttrycket blir då åtminstone \[(n+1)\cdot(n \log n - n + O(\log n)) - 2n\log n - n^2\log n + n^2/2 - n\cdot(n+1) + O(n \log n)\] \[ = n^2/2 + O(n\log n)\] \(B(n)\) växer alltså som \(e^{n^2/2}\). Vilket är snabbt, och inte något som går att faktorisera naivt upp till \(n=10000\).
Tillbaka till problemet
Problemet är egentligen enklare än den där slarviga Stirling-beräkningen ovan. Jag hade bara med den för att jag tänkte på det, och för att det blev mycket symboler.
Formeln för \(B(n)\) har en rekursiv struktur: vi kan skriva \(B(n+1)\) i termer av \(B(n)\): \[B(n+1) = \frac{{(n+1)!}^{n+2}}{\prod_{k=0}^{n+1}(k!)^2} = \frac{ n!^{n+2} \cdot (n+1)^{n+2} }{ \prod_{k=0}^n(k!)^2 \cdot (n+1)!^2 } \] \[=\frac{n!^{n+1}} {\prod_{k=0}^n(k!)^2} \cdot \frac{n!(n+1)^{n+2}}{n!(n+1) \cdot n!(n+1)} = B(n) \cdot \frac{(n+1)^n}{n!} \]
Den rekursiva formeln gör att man enkelt kan hålla koll på talets faktorisering. Säg att vi har en map-struktur som tar primtal på deras multipliciter i faktoriseringen, fct[p] = # antalet ggr p delar B(n). En sån kan enkelt uppdateras med den rekursiva formeln. Så här blir koden:
# Uppdaterar dicts B_fct, factorial_fct till att ha
# faktorisering av B(n+1), (n+1)!
def update_fcts(B_fct, factorial_fct, n):
n_plus_1_fct = sympy.factorint(n+1)
# Vi har B(n+1) = B(n) * (n+1)^n / n!
for p in set(factorial_fct) | set(n_plus_1_fct):
B_fct[p] = B_fct.get(p, 0) + n * n_plus_1_fct.get(p, 0) - factorial_fct.get(p, 0)
for p in set(factorial_fct) | set(n_plus_1_fct):
factorial_fct[p] = factorial_fct.get(p, 0) + n_plus_1_fct.get(p, 0)Om man kör den med
b, f = {}, {}
for i in range(1, 6):
print("B({}) = {}, {}! = {}".format(i, b, i, f))
update_fcts(b, f, i)blir resultatet
B(1) = {}, 1! = {}
B(2) = {2: 1}, 2! = {2: 1}
B(3) = {2: 0, 3: 2}, 3! = {2: 1, 3: 1}
B(4) = {2: 5, 3: 1}, 4! = {2: 3, 3: 1}
B(5) = {2: 2, 3: 0, 5: 4}, 5! = {2: 3, 3: 1, 5: 1}vilket är rätt.
Summan D(n) av alla delare.
Det är en talteoretiskt multiplikativ funktion, som är ganska enkel att räkna ut från faktoriseringen: \(D(p^a q^b r^c) = \frac{p^{a+1}-1}{p-1} \frac{q^{b+1}-1}{q-1} \frac{r^{c+1}-1}{r-1}\). Vi skulle räkna ut summan av alla \(D\) efter division med \(P=10000000207\). Python har \(pow(base, exp, MOD)\), och division med \(p-1\) kan man göra om till multiplikation med \(pow(p-1, P-2, P)\) genom Fermats lilla sats.
Såhär gjorde jag:
inv = {}
for p in sympy.primerange(1, LIM):
inv[p-1] = pow(p-1, PRIME-2, PRIME)
def compute_D(fct):
res = 1
for p in fct:
res *= (pow(p, fct[p]+1, PRIME) - 1) # p^(fct[p]+1) - 1
res *= inv[p-1] # MOD-invers av p-1 MOD PRIME
res %= PRIME
return resDet tar 15.5 sekunder att lägga ihop alla \(D(n)\) upp till n=10000. När jag tog bort skapandet av set i looparna, tog det 14.5 sekunder. Jag tror att det framför allt för att faktoriseringen med Sympy är långsam.