Problem 559 asks for the value of \(Q(50000)\) modulo the prime
$$M=1000000123.$$
The matrix-counting statement is encoded through auxiliary quantities \(P(k,r,n)\). The C++, Python, and Java implementations never try to enumerate the matrices directly. Instead, they transform the count into a coefficient-extraction problem built from inverse factorial powers, alternating signs, and one outer sum over \(k\).
Fix \(n\), \(r\), and \(k\). The implementations first compute a normalized core value and multiply by \((n!)^r\) only at the end. That normalization turns the combinatorial problem into a clean reciprocal-series recurrence.
Because \(M\) is prime and the computation only needs factorials up to \(n=50000\lt M\), every \(m!\) is invertible modulo \(M\). Define
$$u_m=(m!)^{-r}\pmod{M}\qquad (0\le m\le n).$$
These numbers are the basic weights of the method. All divisions by factorials are replaced by modular inverses, and the missing factor \((n!)^r\) is restored only after the normalized coefficient has been computed.
Write
$$n=bk+s,\qquad b=\left\lfloor\frac{n}{k}\right\rfloor,\qquad 0\le s\lt k.$$
The recurrence advances in jumps of size \(k\). The quotient \(b\) counts how many full \(k\)-steps fit into \(n\), while the remainder \(s\) records the final incomplete part. This is why the formulas separate naturally into a full-block part and a remainder correction.
Introduce the formal power series
$$F_k(z)=\sum_{j\ge 0}(-1)^j u_{jk}z^j=1-u_k z+u_{2k}z^2-u_{3k}z^3+\cdots.$$
Now define its reciprocal
$$A_k(z)=\frac{1}{F_k(z)}=\sum_{i\ge 0} a_i z^i.$$
Comparing coefficients in \(A_k(z)F_k(z)=1\) gives
$$a_0=1,$$
$$a_i=\sum_{j=1}^{i}(-1)^{j+1}a_{i-j}u_{jk}\qquad (i\ge 1).$$
This alternating convolution is the heart of the algorithm. Once the sequence \(a_0,a_1,\dots,a_b\) is known, the rest is just a short finishing sum.
If \(s=0\), the normalized core contribution is simply
$$\gamma(n,k,r)=a_b.$$
If \(s\gt 0\), one shifted series is needed:
$$G_{k,s}(z)=\sum_{j\ge 0}(-1)^j u_{s+jk}z^j.$$
The desired core is then the coefficient of \(z^b\) in the product \(A_k(z)G_{k,s}(z)\):
$$\gamma(n,k,r)=[z^b](A_k(z)G_{k,s}(z))=\sum_{j=0}^{b}(-1)^j a_{b-j}u_{s+jk}.$$
So the same recurrence handles all complete \(k\)-blocks, and the shifted tail series accounts for the leftover \(s\) entries.
After the normalized core has been computed, the required count is
$$P(k,r,n)=(n!)^r\gamma(n,k,r)\pmod{M}.$$
For the Project Euler target, the implementations set \(r=n\) and sum over every \(k\):
$$Q(n)=(n!)^n\sum_{k=1}^{n}\gamma(n,k,n)\pmod{M}.$$
The full problem is therefore reduced to evaluating one reciprocal-series coefficient for each \(k\), summing those normalized cores, and multiplying once at the end.
Take \(k=1\), \(r=2\), and \(n=3\). Then \(b=3\) and \(s=0\). The relevant weights are
$$u_1=\frac{1}{(1!)^2}=1,\qquad u_2=\frac{1}{(2!)^2}=\frac14,\qquad u_3=\frac{1}{(3!)^2}=\frac1{36}.$$
Now apply the recurrence:
$$a_0=1,$$
$$a_1=a_0u_1=1,$$
$$a_2=a_1u_1-a_0u_2=1-\frac14=\frac34,$$
$$a_3=a_2u_1-a_1u_2+a_0u_3=\frac34-\frac14+\frac1{36}=\frac{19}{36}.$$
Because \(s=0\), we have \(\gamma(3,1,2)=a_3\). Therefore
$$P(1,2,3)=(3!)^2\cdot \frac{19}{36}=36\cdot \frac{19}{36}=19,$$
which matches the built-in checkpoint used by the implementations.
The C++, Python, and Java implementations follow the same mathematics. They first precompute factorials and inverse factorials, then raise each inverse factorial to the exponent \(r\) to obtain the table \(u_m\). For a fixed \(k\), they fill the coefficient sequence \(a_0,a_1,\dots,a_b\) with the alternating convolution above and evaluate the shifted finishing sum when \(s>0\).
To compute \(P(k,r,n)\), the implementation multiplies the normalized core by \((n!)^r\). To compute \(Q(n)\), it uses \(r=n\), repeats the same per-\(k\) core computation for every \(1\le k\le n\), adds the cores, and only then multiplies by \((n!)^n\). The C++ version parallelizes the outer loop over \(k\), the Java version runs the same logic serially, and the Python version acts as a thin bridge that invokes the compiled C++ solver and returns its numeric result.
Precomputing factorials and inverse factorials costs \(O(n)\). Building the table \(u_m=(m!)^{-r}\) with repeated fast exponentiation costs \(O(n\log r)\). For one fixed \(k\), the recurrence length is \(b=\lfloor n/k\rfloor\), and the nested alternating sum costs \(O(b^2)\) time. Therefore
$$\sum_{k=1}^{n} O\!\left(\left\lfloor\frac{n}{k}\right\rfloor^2\right)=O(n^2),$$
so \(Q(n)\) is evaluated in overall \(O(n^2)\) time. The main tables use \(O(n)\) memory; the multithreaded C++ version adds one reusable work array per worker thread.
Problem 559 verlangt den Wert von \(Q(50000)\) modulo der Primzahl
$$M=1000000123.$$
Die eigentliche Matrixzählung wird über Hilfsgrößen \(P(k,r,n)\) formuliert. Die C++, Python- und Java-Implementierungen enumerieren die Matrizen nicht direkt, sondern wandeln die Aufgabe in ein Koeffizientenproblem mit inversen Fakultätspotenzen, alternierenden Vorzeichen und einer äußeren Summe über \(k\) um.
Fixiere \(n\), \(r\) und \(k\). Zuerst wird ein normierter Kernwert berechnet; der Faktor \((n!)^r\) kommt erst am Ende zurück. Dadurch wird die kombinatorische Struktur zu einer sauberen Rekursion für eine reziproke Reihe.
Da \(M\) prim ist und nur Fakultäten bis \(n=50000\lt M\) benötigt werden, ist jedes \(m!\) modulo \(M\) invertierbar. Definiere
$$u_m=(m!)^{-r}\pmod{M}\qquad (0\le m\le n).$$
Diese Werte sind die Grundgewichte der Methode. Alle Divisionen durch Fakultäten werden durch modulare Inverse ersetzt, und der fehlende Faktor \((n!)^r\) wird erst nach der normierten Koeffizientenberechnung wieder eingesetzt.
Schreibe
$$n=bk+s,\qquad b=\left\lfloor\frac{n}{k}\right\rfloor,\qquad 0\le s\lt k.$$
Die Rekursion springt immer in Schritten der Größe \(k\). Der Quotient \(b\) zählt also die vollständigen \(k\)-Schritte, und der Rest \(s\) beschreibt den letzten unvollständigen Teil. Deshalb zerfallen die Formeln natürlich in einen Vollblock-Teil und eine Restkorrektur.
Führe die formale Potenzreihe
$$F_k(z)=\sum_{j\ge 0}(-1)^j u_{jk}z^j=1-u_k z+u_{2k}z^2-u_{3k}z^3+\cdots$$
ein. Definiere dann ihre Kehrreihe
$$A_k(z)=\frac{1}{F_k(z)}=\sum_{i\ge 0} a_i z^i.$$
Aus dem Koeffizientenvergleich in \(A_k(z)F_k(z)=1\) folgt
$$a_0=1,$$
$$a_i=\sum_{j=1}^{i}(-1)^{j+1}a_{i-j}u_{jk}\qquad (i\ge 1).$$
Diese alternierende Faltung ist der Kern des Algorithmus. Sobald \(a_0,a_1,\dots,a_b\) bekannt sind, bleibt nur noch eine kurze Endsumme.
Falls \(s=0\), ist der normierte Kernbeitrag einfach
$$\gamma(n,k,r)=a_b.$$
Falls \(s>0\), braucht man zusätzlich die verschobene Reihe
$$G_{k,s}(z)=\sum_{j\ge 0}(-1)^j u_{s+jk}z^j.$$
Der gesuchte Kern ist dann der Koeffizient von \(z^b\) im Produkt \(A_k(z)G_{k,s}(z)\):
$$\gamma(n,k,r)=[z^b](A_k(z)G_{k,s}(z))=\sum_{j=0}^{b}(-1)^j a_{b-j}u_{s+jk}.$$
Damit verarbeitet dieselbe Rekursion alle vollständigen \(k\)-Blöcke, und die verschobene Restreihe erfasst genau die letzten \(s\) Elemente.
Nachdem der normierte Kern berechnet ist, lautet die gesuchte Anzahl
$$P(k,r,n)=(n!)^r\gamma(n,k,r)\pmod{M}.$$
Für das Project-Euler-Ziel setzen die Implementierungen \(r=n\) und summieren über alle \(k\):
$$Q(n)=(n!)^n\sum_{k=1}^{n}\gamma(n,k,n)\pmod{M}.$$
Das gesamte Problem reduziert sich also auf einen reziproken Reihenkoeffizienten pro \(k\), danach auf eine einzige modulare Gesamtsumme.
Nimm \(k=1\), \(r=2\) und \(n=3\). Dann gilt \(b=3\) und \(s=0\). Die relevanten Gewichte sind
$$u_1=\frac{1}{(1!)^2}=1,\qquad u_2=\frac{1}{(2!)^2}=\frac14,\qquad u_3=\frac{1}{(3!)^2}=\frac1{36}.$$
Nun wendet man die Rekursion an:
$$a_0=1,$$
$$a_1=a_0u_1=1,$$
$$a_2=a_1u_1-a_0u_2=1-\frac14=\frac34,$$
$$a_3=a_2u_1-a_1u_2+a_0u_3=\frac34-\frac14+\frac1{36}=\frac{19}{36}.$$
Wegen \(s=0\) ist \(\gamma(3,1,2)=a_3\). Daher
$$P(1,2,3)=(3!)^2\cdot \frac{19}{36}=36\cdot \frac{19}{36}=19,$$
genau der eingebaute Kontrollwert der Implementierungen.
Die C++, Python- und Java-Implementierungen folgen derselben Mathematik. Zuerst werden Fakultäten und inverse Fakultäten vorab berechnet. Danach wird jede inverse Fakultät auf den Exponenten \(r\) potenziert, um die Tabelle \(u_m\) zu erhalten. Für festes \(k\) füllen die Programme die Koeffizientenfolge \(a_0,a_1,\dots,a_b\) mit der alternierenden Faltung und werten bei \(s>0\) zusätzlich die verschobene Endsumme aus.
Für \(P(k,r,n)\) multipliziert die Implementierung den normierten Kern mit \((n!)^r\). Für \(Q(n)\) wird \(r=n\) verwendet, derselbe Kern für jedes \(1\le k\le n\) berechnet, alles aufsummiert und erst dann mit \((n!)^n\) multipliziert. Die C++-Version parallelisiert die äußere Schleife über \(k\), die Java-Version führt dieselbe Logik seriell aus, und die Python-Version fungiert als schlanke Brücke zum kompilierten C++-Löser.
Das Vorberechnen der Fakultäten und inversen Fakultäten kostet \(O(n)\). Das Aufbauen der Tabelle \(u_m=(m!)^{-r}\) per schneller modularer Potenzierung kostet \(O(n\log r)\). Für festes \(k\) hat die Rekursion die Länge \(b=\lfloor n/k\rfloor\), und die verschachtelte alternierende Summe benötigt \(O(b^2)\) Zeit. Deshalb gilt
$$\sum_{k=1}^{n} O\!\left(\left\lfloor\frac{n}{k}\right\rfloor^2\right)=O(n^2),$$
sodass \(Q(n)\) insgesamt in \(O(n^2)\) Zeit berechnet wird. Die Haupttabellen brauchen \(O(n)\) Speicher; die mehrfädige C++-Version ergänzt pro Worker-Thread ein wiederverwendbares Arbeitsfeld.
Problem 559,
$$M=1000000123$$
asal modülü altında \(Q(50000)\) değerini ister. Asıl matris sayımı, yardımcı \(P(k,r,n)\) nicelikleri üzerinden kodlanır. C++, Python ve Java uygulamaları matrisleri doğrudan üretmez; bunun yerine sayımı, ters faktöriyel kuvvetlerinden oluşan ağırlıklar, dönüşümlü işaretler ve \(k\) üzerinde tek bir dış toplam ile kurulan bir katsayı çıkarma problemine dönüştürür.
\(n\), \(r\) ve \(k\) sabit olsun. Uygulamalar önce normalize edilmiş bir çekirdek değer hesaplar, \((n!)^r\) çarpanı ise en sonda geri eklenir. Bu normalizasyon, kombinatorik yapıyı açık bir ters-seri bağıntısına indirger.
\(M\) asal olduğu ve hesaplama yalnızca \(n=50000\lt M\) için gereken faktöriyellere kadar çıktığı için her \(m!\) değeri modulo \(M\) içinde terslenebilir. Şunu tanımlayalım:
$$u_m=(m!)^{-r}\pmod{M}\qquad (0\le m\le n).$$
Bu sayılar yöntemin temel ağırlıklarıdır. Faktöriyelle bölme işlemleri modüler terslerle temsil edilir; eksik kalan \((n!)^r\) çarpanı ise normalize katsayı bulunduktan sonra geri konur.
Şöyle yazalım:
$$n=bk+s,\qquad b=\left\lfloor\frac{n}{k}\right\rfloor,\qquad 0\le s\lt k.$$
Bağıntı yalnızca \(k\) büyüklüğünde sıçramalarla ilerler. Bu yüzden \(b\), kaç tam \(k\)-adımı olduğunu; \(s\) ise en sondaki eksik parçayı gösterir. Formüllerin bir tam-blok kısmı ve bir artık düzeltmesi olarak ayrılması tam olarak bundan kaynaklanır.
Önce şu formal kuvvet serisini tanımlayalım:
$$F_k(z)=\sum_{j\ge 0}(-1)^j u_{jk}z^j=1-u_k z+u_{2k}z^2-u_{3k}z^3+\cdots.$$
Şimdi bunun tersini alalım:
$$A_k(z)=\frac{1}{F_k(z)}=\sum_{i\ge 0} a_i z^i.$$
\(A_k(z)F_k(z)=1\) eşitliğinde katsayıları karşılaştırınca
$$a_0=1,$$
$$a_i=\sum_{j=1}^{i}(-1)^{j+1}a_{i-j}u_{jk}\qquad (i\ge 1)$$
elde edilir. Algoritmanın kalbi bu dönüşümlü evrişimdir. \(a_0,a_1,\dots,a_b\) dizisi hesaplanınca geriye yalnızca kısa bir bitiş toplamı kalır.
Eğer \(s=0\) ise normalize çekirdek doğrudan
$$\gamma(n,k,r)=a_b$$
olur. Eğer \(s>0\) ise bir de kaydırılmış seri gerekir:
$$G_{k,s}(z)=\sum_{j\ge 0}(-1)^j u_{s+jk}z^j.$$
İstenen çekirdek, \(A_k(z)G_{k,s}(z)\) çarpımındaki \(z^b\) katsayısıdır:
$$\gamma(n,k,r)=[z^b](A_k(z)G_{k,s}(z))=\sum_{j=0}^{b}(-1)^j a_{b-j}u_{s+jk}.$$
Böylece aynı bağıntı tüm tam \(k\)-bloklarını işler; kaydırılmış son seri de geride kalan \(s\) elemanlık parçayı hesaba katar.
Normalize çekirdek hesaplandıktan sonra istenen sayı
$$P(k,r,n)=(n!)^r\gamma(n,k,r)\pmod{M}$$
şeklindedir. Project Euler hedefinde uygulamalar \(r=n\) alır ve tüm \(k\) değerleri üzerinde toplar:
$$Q(n)=(n!)^n\sum_{k=1}^{n}\gamma(n,k,n)\pmod{M}.$$
Dolayısıyla tüm problem, her \(k\) için bir ters-seri katsayısı hesaplamaya, bu normalize çekirdekleri toplamaya ve sonda tek bir çarpım yapmaya indirgenir.
\(k=1\), \(r=2\), \(n=3\) alalım. Bu durumda \(b=3\) ve \(s=0\). Gerekli ağırlıklar
$$u_1=\frac{1}{(1!)^2}=1,\qquad u_2=\frac{1}{(2!)^2}=\frac14,\qquad u_3=\frac{1}{(3!)^2}=\frac1{36}$$
olur. Şimdi bağıntıyı uygularız:
$$a_0=1,$$
$$a_1=a_0u_1=1,$$
$$a_2=a_1u_1-a_0u_2=1-\frac14=\frac34,$$
$$a_3=a_2u_1-a_1u_2+a_0u_3=\frac34-\frac14+\frac1{36}=\frac{19}{36}.$$
\(s=0\) olduğu için \(\gamma(3,1,2)=a_3\) olur. O halde
$$P(1,2,3)=(3!)^2\cdot \frac{19}{36}=36\cdot \frac{19}{36}=19,$$
ve bu, uygulamalardaki yerleşik kontrol değeriyle aynıdır.
C++, Python ve Java uygulamaları aynı matematiği izler. Önce faktöriyeller ve modüler ters faktöriyeller önceden hazırlanır. Ardından her ters faktöriyel, gerekli \(r\) üssüne yükseltilerek \(u_m\) tablosu üretilir. Sabit bir \(k\) için programlar yukarıdaki dönüşümlü evrişimle \(a_0,a_1,\dots,a_b\) dizisini doldurur ve \(s>0\) ise kaydırılmış bitiş toplamını hesaplar.
\(P(k,r,n)\) için uygulama normalize çekirdeği \((n!)^r\) ile çarpar. \(Q(n)\) için \(r=n\) alınır, aynı çekirdek hesabı her \(1\le k\le n\) için tekrarlanır, çekirdekler toplanır ve en sonda \((n!)^n\) ile çarpılır. C++ sürümü dıştaki \(k\) döngüsünü paralelleştirir, Java sürümü aynı mantığı tek iş parçacığında çalıştırır, Python sürümü ise derlenmiş C++ çözücüsünü çağıran ince bir köprü görevi görür.
Faktöriyel ve ters faktöriyel tablolarını hazırlamak \(O(n)\) zaman alır. \(u_m=(m!)^{-r}\) tablosunu tekrarlı hızlı üs alma ile kurmak \(O(n\log r)\) maliyetlidir. Sabit bir \(k\) için bağıntı uzunluğu \(b=\lfloor n/k\rfloor\) olur ve iç içe dönüşümlü toplam \(O(b^2)\) zaman ister. Bu nedenle
$$\sum_{k=1}^{n} O\!\left(\left\lfloor\frac{n}{k}\right\rfloor^2\right)=O(n^2)$$
elde edilir; yani \(Q(n)\) toplamda \(O(n^2)\) zamanda hesaplanır. Ana tablolar \(O(n)\) bellek kullanır; çok iş parçacıklı C++ sürümü buna ek olarak her işçi için yeniden kullanılan bir çalışma dizisi tutar.
El problema 559 pide el valor de \(Q(50000)\) módulo el primo
$$M=1000000123.$$
La formulación de conteo de matrices se codifica mediante las cantidades auxiliares \(P(k,r,n)\). Las implementaciones en C++, Python y Java no intentan enumerar matrices de forma directa; transforman el conteo en un problema de extracción de coeficientes con potencias de factoriales inversos, signos alternantes y una suma exterior sobre \(k\).
Fijemos \(n\), \(r\) y \(k\). Primero se calcula un núcleo normalizado y solo al final se vuelve a multiplicar por \((n!)^r\). Esa normalización convierte el problema combinatorio en una recurrencia limpia de series recíprocas.
Como \(M\) es primo y la computación solo necesita factoriales hasta \(n=50000\lt M\), cada \(m!\) es invertible módulo \(M\). Definimos
$$u_m=(m!)^{-r}\pmod{M}\qquad (0\le m\le n).$$
Estos valores son los pesos básicos del método. Toda división por factoriales se reemplaza por inversos modulares, y el factor faltante \((n!)^r\) se restaura únicamente después de obtener el coeficiente normalizado.
Escribimos
$$n=bk+s,\qquad b=\left\lfloor\frac{n}{k}\right\rfloor,\qquad 0\le s\lt k.$$
La recurrencia avanza en saltos de tamaño \(k\). El cociente \(b\) cuenta cuántos pasos completos caben en \(n\), mientras que el resto \(s\) registra la parte final incompleta. Por eso las fórmulas se dividen de manera natural en una parte de bloques completos y una corrección final.
Introducimos la serie formal
$$F_k(z)=\sum_{j\ge 0}(-1)^j u_{jk}z^j=1-u_k z+u_{2k}z^2-u_{3k}z^3+\cdots.$$
Ahora definimos su recíproca
$$A_k(z)=\frac{1}{F_k(z)}=\sum_{i\ge 0} a_i z^i.$$
Al comparar coeficientes en \(A_k(z)F_k(z)=1\), obtenemos
$$a_0=1,$$
$$a_i=\sum_{j=1}^{i}(-1)^{j+1}a_{i-j}u_{jk}\qquad (i\ge 1).$$
Esta convolución alternante es el núcleo del algoritmo. Una vez conocida la sucesión \(a_0,a_1,\dots,a_b\), solo queda una suma final corta.
Si \(s=0\), el núcleo normalizado es simplemente
$$\gamma(n,k,r)=a_b.$$
Si \(s>0\), hace falta una serie desplazada adicional:
$$G_{k,s}(z)=\sum_{j\ge 0}(-1)^j u_{s+jk}z^j.$$
Entonces el núcleo buscado es el coeficiente de \(z^b\) en el producto \(A_k(z)G_{k,s}(z)\):
$$\gamma(n,k,r)=[z^b](A_k(z)G_{k,s}(z))=\sum_{j=0}^{b}(-1)^j a_{b-j}u_{s+jk}.$$
Así, la misma recurrencia maneja todos los bloques completos de tamaño \(k\), y la serie desplazada corrige exactamente las \(s\) posiciones sobrantes.
Una vez calculado el núcleo normalizado, la cantidad pedida es
$$P(k,r,n)=(n!)^r\gamma(n,k,r)\pmod{M}.$$
Para el objetivo de Project Euler, las implementaciones toman \(r=n\) y suman sobre todos los \(k\):
$$Q(n)=(n!)^n\sum_{k=1}^{n}\gamma(n,k,n)\pmod{M}.$$
En consecuencia, todo el problema se reduce a evaluar un coeficiente de serie recíproca para cada \(k\), sumar esos núcleos normalizados y multiplicar una sola vez al final.
Tomemos \(k=1\), \(r=2\) y \(n=3\). Entonces \(b=3\) y \(s=0\). Los pesos relevantes son
$$u_1=\frac{1}{(1!)^2}=1,\qquad u_2=\frac{1}{(2!)^2}=\frac14,\qquad u_3=\frac{1}{(3!)^2}=\frac1{36}.$$
Aplicamos ahora la recurrencia:
$$a_0=1,$$
$$a_1=a_0u_1=1,$$
$$a_2=a_1u_1-a_0u_2=1-\frac14=\frac34,$$
$$a_3=a_2u_1-a_1u_2+a_0u_3=\frac34-\frac14+\frac1{36}=\frac{19}{36}.$$
Como \(s=0\), tenemos \(\gamma(3,1,2)=a_3\). Por lo tanto
$$P(1,2,3)=(3!)^2\cdot \frac{19}{36}=36\cdot \frac{19}{36}=19,$$
lo que coincide con el valor de comprobación incorporado en las implementaciones.
Las implementaciones en C++, Python y Java siguen la misma matemática. Primero precalculan factoriales e inversos factoriales modulares. Después elevan cada inverso factorial al exponente \(r\) para construir la tabla \(u_m\). Para un \(k\) fijo, generan la sucesión \(a_0,a_1,\dots,a_b\) mediante la convolución alternante anterior y, si \(s>0\), evalúan la suma final desplazada.
Para calcular \(P(k,r,n)\), la implementación multiplica el núcleo normalizado por \((n!)^r\). Para calcular \(Q(n)\), usa \(r=n\), repite el mismo cálculo del núcleo para cada \(1\le k\le n\), suma todos los núcleos y solo entonces multiplica por \((n!)^n\). La versión en C++ paraleliza el bucle exterior sobre \(k\), la versión en Java ejecuta la misma lógica en serie y la versión en Python actúa como un puente ligero que invoca el solucionador compilado en C++ y devuelve el resultado numérico.
Precalcular factoriales e inversos factoriales cuesta \(O(n)\). Construir la tabla \(u_m=(m!)^{-r}\) con exponenciación rápida repetida cuesta \(O(n\log r)\). Para un \(k\) fijo, la longitud de la recurrencia es \(b=\lfloor n/k\rfloor\), y la suma alternante anidada requiere \(O(b^2)\) tiempo. Por tanto,
$$\sum_{k=1}^{n} O\!\left(\left\lfloor\frac{n}{k}\right\rfloor^2\right)=O(n^2),$$
de modo que \(Q(n)\) se obtiene en \(O(n^2)\) tiempo total. Las tablas principales usan \(O(n)\) memoria; la versión multihilo en C++ añade un arreglo de trabajo reutilizable por cada hilo.
第 559 题要求计算
$$Q(50000)\pmod{M},\qquad M=1000000123,$$
其中 \(M\) 是素数。题目的矩阵计数被编码为辅助量 \(P(k,r,n)\)。C++、Python 和 Java 实现都不直接枚举矩阵，而是把问题改写成一个系数提取问题：核心对象是逆阶乘幂、交替符号的卷积，以及对所有 \(k\) 的一次外层求和。
固定 \(n\)、\(r\) 和 \(k\)。实现先求一个“归一化核心值”，最后再乘回 \((n!)^r\)。这样做以后，中间式子会变得非常整齐，本质上就是一个倒数形式幂级数的系数递推。
因为 \(M\) 是素数，而且实际计算只需要到 \(n=50000\lt M\) 为止的阶乘，所以每个 \(m!\) 在模 \(M\) 下都有乘法逆元。定义
$$u_m=(m!)^{-r}\pmod{M}\qquad (0\le m\le n).$$
这些 \(u_m\) 就是整个算法的基本权重。所有“除以阶乘”的操作都被替换成模逆，等到归一化系数算完以后，再统一乘上 \((n!)^r\) 还原成题目需要的计数。
写成
$$n=bk+s,\qquad b=\left\lfloor\frac{n}{k}\right\rfloor,\qquad 0\le s\lt k.$$
递推每次只按大小为 \(k\) 的步长前进，因此 \(b\) 表示能放下多少个完整的 \(k\) 块，而 \(s\) 记录最后剩下的不完整部分。这也是为什么公式会自然拆成“完整块部分”和“尾部修正部分”。
先定义形式幂级数
$$F_k(z)=\sum_{j\ge 0}(-1)^j u_{jk}z^j=1-u_k z+u_{2k}z^2-u_{3k}z^3+\cdots.$$
再取它的倒数：
$$A_k(z)=\frac{1}{F_k(z)}=\sum_{i\ge 0} a_i z^i.$$
把 \(A_k(z)F_k(z)=1\) 的同次项系数逐项比较，可得
$$a_0=1,$$
$$a_i=\sum_{j=1}^{i}(-1)^{j+1}a_{i-j}u_{jk}\qquad (i\ge 1).$$
这条交替卷积递推正是实现真正执行的核心步骤。只要把 \(a_0,a_1,\dots,a_b\) 依次求出，后面的工作就只是一个较短的收尾求和。
如果 \(s=0\)，那么归一化核心值直接就是
$$\gamma(n,k,r)=a_b.$$
如果 \(s>0\)，则还需要一个“平移后”的级数：
$$G_{k,s}(z)=\sum_{j\ge 0}(-1)^j u_{s+jk}z^j.$$
这时目标核心值就是乘积 \(A_k(z)G_{k,s}(z)\) 中 \(z^b\) 的系数：
$$\gamma(n,k,r)=[z^b](A_k(z)G_{k,s}(z))=\sum_{j=0}^{b}(-1)^j a_{b-j}u_{s+jk}.$$
也就是说，完整的 \(k\) 块仍然由同一条递推处理，而最后剩下的 \(s\) 个位置通过这个平移级数被精确修正进去。
当归一化核心值求出以后，真正需要的计数是
$$P(k,r,n)=(n!)^r\gamma(n,k,r)\pmod{M}.$$
对于本题目标，三种实现都取 \(r=n\)，并对所有 \(k\) 求和：
$$Q(n)=(n!)^n\sum_{k=1}^{n}\gamma(n,k,n)\pmod{M}.$$
因此整个问题可以概括成三步：对每个 \(k\) 计算一个倒数幂级数系数，累加这些归一化核心值，最后乘上统一的 \((n!)^n\) 因子。
取 \(k=1\)、\(r=2\)、\(n=3\)。此时 \(b=3\)、\(s=0\)。相关权重为
$$u_1=\frac{1}{(1!)^2}=1,\qquad u_2=\frac{1}{(2!)^2}=\frac14,\qquad u_3=\frac{1}{(3!)^2}=\frac1{36}.$$
按照递推逐步计算：
$$a_0=1,$$
$$a_1=a_0u_1=1,$$
$$a_2=a_1u_1-a_0u_2=1-\frac14=\frac34,$$
$$a_3=a_2u_1-a_1u_2+a_0u_3=\frac34-\frac14+\frac1{36}=\frac{19}{36}.$$
由于 \(s=0\)，所以 \(\gamma(3,1,2)=a_3\)。于是
$$P(1,2,3)=(3!)^2\cdot \frac{19}{36}=36\cdot \frac{19}{36}=19,$$
这与实现内置的校验值完全一致。这个例子也说明了为什么先做归一化会更方便：中间可以用非常短的递推式，最后再把 \((n!)^r\) 一次性乘回去。
C++、Python 和 Java 实现遵循完全相同的数学结构。第一步是预计算阶乘与模逆阶乘。第二步是把每个逆阶乘提升到所需指数 \(r\)，得到权重表 \(u_m\)。接下来对固定的 \(k\)，用上面的交替卷积递推填出 \(a_0,a_1,\dots,a_b\)，如果 \(s>0\)，再额外计算一次平移后的收尾求和。
在求 \(P(k,r,n)\) 时，实现会把归一化核心值乘上 \((n!)^r\)。在求 \(Q(n)\) 时，它们把 \(r\) 设成 \(n\)，对每个 \(1\le k\le n\) 重复同样的核心计算，把所有核心值先加起来，最后才乘上 \((n!)^n\)。其中 C++ 版本把外层的 \(k\) 循环并行化；Java 版本以单线程方式执行同样的流程；Python 版本则不重复实现整套递推，而是负责调用编译后的 C++ 求解器并返回最终数值。
预计算阶乘和逆阶乘需要 \(O(n)\) 时间。利用快速幂构造 \(u_m=(m!)^{-r}\) 表需要 \(O(n\log r)\) 时间。对固定的 \(k\)，递推长度是 \(b=\lfloor n/k\rfloor\)，嵌套的交替求和代价是 \(O(b^2)\)。因此总时间满足
$$\sum_{k=1}^{n} O\!\left(\left\lfloor\frac{n}{k}\right\rfloor^2\right)=O(n^2).$$
也就是说，计算 \(Q(n)\) 的总体时间复杂度为 \(O(n^2)\)。主要表格占用 \(O(n)\) 内存；多线程 C++ 版本还会为每个工作线程保留一个可复用的工作数组。
Задача 559 требует вычислить \(Q(50000)\) по модулю простого числа
$$M=1000000123.$$
Подсчёт матриц переписывается через вспомогательные величины \(P(k,r,n)\). Реализации на C++, Python и Java не перебирают матрицы напрямую. Вместо этого они сводят задачу к извлечению коэффициента из ряда, построенного из степеней обратных факториалов, чередующихся знаков и внешней суммы по всем \(k\).
Зафиксируем \(n\), \(r\) и \(k\). Сначала вычисляется нормированное ядро, а множитель \((n!)^r\) возвращается только в самом конце. Благодаря этому комбинаторная структура превращается в прозрачную рекурсию для обратного степенного ряда.
Так как \(M\) простое и в вычислении нужны только факториалы до \(n=50000\lt M\), каждый \(m!\) обратим по модулю \(M\). Введём
$$u_m=(m!)^{-r}\pmod{M}\qquad (0\le m\le n).$$
Именно эти значения являются базовыми весами метода. Деление на факториалы заменяется модульными обратными, а общий множитель \((n!)^r\) возвращается уже после вычисления нормированного коэффициента.
Запишем
$$n=bk+s,\qquad b=\left\lfloor\frac{n}{k}\right\rfloor,\qquad 0\le s\lt k.$$
Рекурсия движется шагами длины \(k\). Число \(b\) показывает, сколько полных \(k\)-блоков помещается в \(n\), а остаток \(s\) отвечает за финальную неполную часть. Поэтому формулы естественно распадаются на полноблочную часть и отдельную поправку для хвоста.
Введём формальный степенной ряд
$$F_k(z)=\sum_{j\ge 0}(-1)^j u_{jk}z^j=1-u_k z+u_{2k}z^2-u_{3k}z^3+\cdots.$$
Теперь рассмотрим его обратный ряд
$$A_k(z)=\frac{1}{F_k(z)}=\sum_{i\ge 0} a_i z^i.$$
Сравнение коэффициентов в равенстве \(A_k(z)F_k(z)=1\) даёт
$$a_0=1,$$
$$a_i=\sum_{j=1}^{i}(-1)^{j+1}a_{i-j}u_{jk}\qquad (i\ge 1).$$
Эта знакопеременная свёртка и есть центральная часть алгоритма. После того как найдены \(a_0,a_1,\dots,a_b\), остаётся только короткая завершающая сумма.
Если \(s=0\), нормированное ядро равно
$$\gamma(n,k,r)=a_b.$$
Если \(s>0\), требуется ещё один сдвинутый ряд:
$$G_{k,s}(z)=\sum_{j\ge 0}(-1)^j u_{s+jk}z^j.$$
Тогда нужное ядро есть коэффициент при \(z^b\) в произведении \(A_k(z)G_{k,s}(z)\):
$$\gamma(n,k,r)=[z^b](A_k(z)G_{k,s}(z))=\sum_{j=0}^{b}(-1)^j a_{b-j}u_{s+jk}.$$
Тем самым одна и та же рекурсия обрабатывает все полные \(k\)-блоки, а сдвинутый хвост точно учитывает оставшиеся \(s\) позиций.
После вычисления нормированного ядра искомое количество равно
$$P(k,r,n)=(n!)^r\gamma(n,k,r)\pmod{M}.$$
Для целевой задачи реализации берут \(r=n\) и суммируют по всем \(k\):
$$Q(n)=(n!)^n\sum_{k=1}^{n}\gamma(n,k,n)\pmod{M}.$$
Итак, вся задача сводится к вычислению одного коэффициента обратного ряда для каждого \(k\), суммированию нормированных ядер и одному финальному умножению.
Возьмём \(k=1\), \(r=2\), \(n=3\). Тогда \(b=3\) и \(s=0\). Нужные веса:
$$u_1=\frac{1}{(1!)^2}=1,\qquad u_2=\frac{1}{(2!)^2}=\frac14,\qquad u_3=\frac{1}{(3!)^2}=\frac1{36}.$$
Применяем рекурсию:
$$a_0=1,$$
$$a_1=a_0u_1=1,$$
$$a_2=a_1u_1-a_0u_2=1-\frac14=\frac34,$$
$$a_3=a_2u_1-a_1u_2+a_0u_3=\frac34-\frac14+\frac1{36}=\frac{19}{36}.$$
Так как \(s=0\), получаем \(\gamma(3,1,2)=a_3\). Следовательно,
$$P(1,2,3)=(3!)^2\cdot \frac{19}{36}=36\cdot \frac{19}{36}=19,$$
что совпадает со встроенной проверкой в реализациях.
Реализации на C++, Python и Java используют одну и ту же математическую схему. Сначала они заранее вычисляют факториалы и обратные факториалы по модулю. Затем каждый обратный факториал возводится в степень \(r\), чтобы получить таблицу \(u_m\). Для фиксированного \(k\) программы строят последовательность \(a_0,a_1,\dots,a_b\) по знакопеременной свёртке, а при \(s>0\) дополнительно вычисляют сдвинутую завершающую сумму.
Чтобы найти \(P(k,r,n)\), реализация умножает нормированное ядро на \((n!)^r\). Чтобы найти \(Q(n)\), берётся \(r=n\), один и тот же расчёт ядра повторяется для каждого \(1\le k\le n\), все ядра суммируются, и только потом результат умножается на \((n!)^n\). Версия на C++ распараллеливает внешний цикл по \(k\), версия на Java выполняет ту же логику последовательно, а версия на Python служит тонким мостом: она запускает скомпилированный решатель на C++ и возвращает итоговое число.
Предварительное вычисление факториалов и обратных факториалов занимает \(O(n)\). Построение таблицы \(u_m=(m!)^{-r}\) с помощью повторной быстрой экспоненты занимает \(O(n\log r)\). Для фиксированного \(k\) длина рекурсии равна \(b=\lfloor n/k\rfloor\), а вложенная знакопеременная сумма стоит \(O(b^2)\) времени. Поэтому
$$\sum_{k=1}^{n} O\!\left(\left\lfloor\frac{n}{k}\right\rfloor^2\right)=O(n^2),$$
и общая трудоёмкость вычисления \(Q(n)\) составляет \(O(n^2)\). Основные таблицы требуют \(O(n)\) памяти; многопоточная версия на C++ добавляет по одному переиспользуемому рабочему массиву на поток.
المطلوب في المسألة 559 هو حساب \(Q(50000)\) بترديد العدد الأولي
$$M=1000000123.$$
صياغة عدّ المصفوفات تُكتب عبر الكميات المساعدة \(P(k,r,n)\). لا تحاول تطبيقات ++C وPython وJava تعداد المصفوفات مباشرة، بل تحوّل العد إلى مسألة استخراج معامل من متسلسلة مبنية على قوى مقلوبات المضروب، مع إشارات متناوبة ومجموع خارجي على جميع قيم \(k\).
لنثبّت \(n\) و\(r\) و\(k\). تحسب التطبيقات أولاً نواة مُطبَّعة، ثم تعيد ضرب النتيجة في \((n!)^r\) في النهاية فقط. هذا التطبيع يجعل البنية التوافقية تظهر على شكل علاقة تكرارية واضحة لسلسلة مقلوبة.
بما أن \(M\) أولي، ولأن الحساب لا يحتاج إلا إلى مضاريب حتى \(n=50000\lt M\)، فإن كل \(m!\) قابل للعكس بترديد \(M\). نعرّف
$$u_m=(m!)^{-r}\pmod{M}\qquad (0\le m\le n).$$
هذه القيم هي الأوزان الأساسية في الطريقة. كل قسمة على مضروب تُستبدل بمعكوس ضربي بترديد \(M\)، ثم يُعاد العامل \((n!)^r\) بعد انتهاء حساب المعامل المُطبَّع.
نكتب
$$n=bk+s,\qquad b=\left\lfloor\frac{n}{k}\right\rfloor,\qquad 0\le s\lt k.$$
العلاقة التكرارية تتحرك بخطوات طولها \(k\). لذلك فإن \(b\) يساوي عدد الكتل الكاملة، بينما يسجل \(s\) الجزء الأخير غير المكتمل. ولهذا السبب تنقسم الصيغ طبيعياً إلى جزء خاص بالكتل الكاملة وجزء آخر لتصحيح الباقي.
لنعرّف السلسلة الشكلية
$$F_k(z)=\sum_{j\ge 0}(-1)^j u_{jk}z^j=1-u_k z+u_{2k}z^2-u_{3k}z^3+\cdots.$$
ثم نأخذ مقلوبها:
$$A_k(z)=\frac{1}{F_k(z)}=\sum_{i\ge 0} a_i z^i.$$
بمقارنة المعاملات في العلاقة \(A_k(z)F_k(z)=1\) نحصل على
$$a_0=1,$$
$$a_i=\sum_{j=1}^{i}(-1)^{j+1}a_{i-j}u_{jk}\qquad (i\ge 1).$$
هذا الالتفاف ذو الإشارات المتناوبة هو قلب الخوارزمية. وما إن تُحسَب القيم \(a_0,a_1,\dots,a_b\)، حتى لا يبقى إلا مجموع ختامي قصير.
إذا كان \(s=0\)، فإن النواة المُطبَّعة تساوي مباشرة
$$\gamma(n,k,r)=a_b.$$
أما إذا كان \(s>0\)، فنحتاج إلى سلسلة مزاحة إضافية:
$$G_{k,s}(z)=\sum_{j\ge 0}(-1)^j u_{s+jk}z^j.$$
وعندئذ تكون النواة المطلوبة هي معامل \(z^b\) في حاصل الضرب \(A_k(z)G_{k,s}(z)\):
$$\gamma(n,k,r)=[z^b](A_k(z)G_{k,s}(z))=\sum_{j=0}^{b}(-1)^j a_{b-j}u_{s+jk}.$$
وهكذا تعالج العلاقة التكرارية جميع الكتل الكاملة، بينما تُدخِل السلسلة المزاحة أثر العناصر \(s\) المتبقية في النهاية.
بعد حساب النواة المُطبَّعة تصبح الكمية المطلوبة
$$P(k,r,n)=(n!)^r\gamma(n,k,r)\pmod{M}.$$
وفي الهدف الخاص بمسألة Project Euler تأخذ التطبيقات \(r=n\) ثم تجمع على جميع قيم \(k\):
$$Q(n)=(n!)^n\sum_{k=1}^{n}\gamma(n,k,n)\pmod{M}.$$
إذن المسألة كلها تختزل إلى حساب معامل واحد من سلسلة مقلوبة لكل \(k\)، ثم جمع هذه الأنوية المُطبَّعة، ثم إجراء ضرب نهائي واحد.
خذ \(k=1\) و\(r=2\) و\(n=3\). عندئذٍ \(b=3\) و\(s=0\). الأوزان اللازمة هي
$$u_1=\frac{1}{(1!)^2}=1,\qquad u_2=\frac{1}{(2!)^2}=\frac14,\qquad u_3=\frac{1}{(3!)^2}=\frac1{36}.$$
نطبّق الآن العلاقة التكرارية:
$$a_0=1,$$
$$a_1=a_0u_1=1,$$
$$a_2=a_1u_1-a_0u_2=1-\frac14=\frac34,$$
$$a_3=a_2u_1-a_1u_2+a_0u_3=\frac34-\frac14+\frac1{36}=\frac{19}{36}.$$
وبما أن \(s=0\)، فإن \(\gamma(3,1,2)=a_3\). لذلك
$$P(1,2,3)=(3!)^2\cdot \frac{19}{36}=36\cdot \frac{19}{36}=19,$$
وهو تماماً نفس مقدار التحقق الداخلي المستخدم في التطبيقات.
تتبع تطبيقات ++C وPython وJava البنية الرياضية نفسها. تبدأ بحساب المضاريب ومقلوبات المضاريب بترديد \(M\)، ثم ترفع كل مقلوب مضروب إلى الأس \(r\) لبناء الجدول \(u_m\). بعد ذلك، ومن أجل \(k\) ثابت، تُولَّد المتتالية \(a_0,a_1,\dots,a_b\) بواسطة الالتفاف المتناوب السابق، وإذا كان \(s>0\) تُحسب أيضاً النهاية المزاحة.
لحساب \(P(k,r,n)\)، تضرب الشيفرة النواة المُطبَّعة في \((n!)^r\). ولحساب \(Q(n)\)، تضع \(r=n\)، وتكرر الحساب نفسه لكل \(1\le k\le n\)، ثم تجمع جميع الأنوية قبل أن تضرب مرة واحدة في \((n!)^n\). نسخة ++C توازي الحلقة الخارجية على \(k\)، ونسخة Java تنفذ المنطق نفسه على نحو تسلسلي، أما نسخة Python فتعمل كجسر خفيف يستدعي محلل ++C المترجم ويعيد الناتج العددي.
يساوي بناء جداول المضروب ومقلوبه \(O(n)\). أما إنشاء الجدول \(u_m=(m!)^{-r}\) باستخدام الأس السريع المتكرر فيكلف \(O(n\log r)\). ولـ \(k\) ثابت، يكون طول العلاقة التكرارية \(b=\lfloor n/k\rfloor\)، بينما تكلف الحلقة المتداخلة ذات الإشارات المتناوبة \(O(b^2)\) زمناً. ومن ثم
$$\sum_{k=1}^{n} O\!\left(\left\lfloor\frac{n}{k}\right\rfloor^2\right)=O(n^2),$$
أي إن حساب \(Q(n)\) كله يتم في \(O(n^2)\) زمنياً. الجداول الرئيسية تستهلك \(O(n)\) من الذاكرة، وتضيف نسخة ++C متعددة الخيوط مصفوفة عمل قابلة لإعادة الاستخدام لكل خيط.