Problem Summary

For each positive integer \(n\), define the Gauss factorial

$$\operatorname{GF}(n)=\prod_{\substack{1\le k\lt n\\ \gcd(k,n)=1}} k,$$

with \(\operatorname{GF}(1)=1\). The problem asks for

$$P(N)=\prod_{n=1}^{N}\operatorname{GF}(n)\pmod{10^9+7}$$

at \(N=10^8\). A direct computation would enumerate far too many coprime pairs, so the solution rewrites the whole product into a divisor-based formula that can be evaluated in linear time.

Mathematical Approach

The key idea is to convert each Gauss factorial into a Möbius-inverted divisor product, then combine all \(n\le N\) into blocks where \(\left\lfloor N/d\right\rfloor\) is constant.

Step 1: Group the Factors of \((n-1)!\) by Their GCD with \(n\)

Fix \(n\). Every integer \(k\) with \(1\le k\lt n\) has some gcd \(d=\gcd(k,n)\), so it can be written as

$$k=d\,m,\qquad \gcd\!\left(m,\frac{n}{d}\right)=1,\qquad 1\le m\lt\frac{n}{d}.$$

For a fixed divisor \(d\mid n\), the product of all such terms is

$$d^{\phi(n/d)}\operatorname{GF}\!\left(\frac{n}{d}\right).$$

Multiplying over all possible gcd values yields

$$ (n-1)! = \prod_{\substack{d \mid n\\ d\lt n}} d^{\phi(n/d)}\operatorname{GF}\!\left(\frac{n}{d}\right). $$

After the substitution \(m=n/d\), this becomes

$$ (n-1)! = \prod_{\substack{m \mid n\\ m\gt 1}} \left(\frac{n}{m}\right)^{\phi(m)}\operatorname{GF}(m). $$

Step 2: Normalize and Apply Möbius Inversion

Because

$$\sum_{\substack{m \mid n\\ m\gt 1}}\phi(m)=n-1,$$

we can divide both sides by \(n^{n-1}\) and rewrite the identity as

$$ \frac{(n-1)!}{n^{n-1}}=\prod_{m \mid n}\frac{\operatorname{GF}(m)}{m^{\phi(m)}}. $$

Define

$$ A(n)=\frac{(n-1)!}{n^{n-1}},\qquad H(n)=\frac{\operatorname{GF}(n)}{n^{\phi(n)}}. $$

Then \(A(n)=\prod_{d\mid n}H(d)\), so the multiplicative form of Möbius inversion gives

$$ H(n)=\prod_{d \mid n} A(d)^{\mu(n/d)} =\prod_{d \mid n}\left(\frac{(d-1)!}{d^{d-1}}\right)^{\mu(n/d)}. $$

Multiplying back by \(n^{\phi(n)}\) and simplifying the powers of \(n\) gives a more implementation-friendly identity:

$$ \operatorname{GF}(n)=\prod_{d \mid n}\left(\left(\frac{n}{d}-1\right)! \, d^{\,n/d-1}\right)^{\mu(d)}. $$

This formula is valid for \(n=1\) as well, because the only divisor is \(d=1\) and the factor is \(0!\cdot 1^0=1\).

Step 3: Multiply over All \(n\le N\) and Swap the Order

Now multiply the previous identity for every \(1\le n\le N\):

$$ P(N)=\prod_{n=1}^{N}\prod_{d \mid n}\left(\left(\frac{n}{d}-1\right)! \, d^{\,n/d-1}\right)^{\mu(d)}. $$

Write \(n=d\,m\). Then the outer product becomes

$$ P(N)=\prod_{d=1}^{N}\prod_{m=1}^{\lfloor N/d\rfloor}\left((m-1)! \, d^{\,m-1}\right)^{\mu(d)}. $$

For a fixed \(d\), let \(q=\left\lfloor N/d\right\rfloor\). The exponent of \(d\) is

$$ \sum_{m=1}^{q}(m-1)=\frac{q(q-1)}{2}=\binom{q}{2}, $$

and the factorial part is

$$ \prod_{m=1}^{q}(m-1)! = \prod_{t=1}^{q-1} t!. $$

Step 4: Introduce the Superfactorial

Define the superfactorial

$$ \operatorname{SF}(k)=\prod_{t=1}^{k} t!,\qquad \operatorname{SF}(0)=1. $$

Then the whole product collapses to

$$ \boxed{ P(N)\equiv \prod_{d=1}^{N}\left(d^{\binom{\lfloor N/d\rfloor}{2}}\operatorname{SF}\!\left(\lfloor N/d\rfloor-1\right)\right)^{\mu(d)} \pmod{10^9+7}. } $$

This is exactly the mathematical formula implemented by the programs.

Step 5: Group Equal Floor Quotients

The quantity

$$q=\left\lfloor \frac{N}{d}\right\rfloor$$

is constant on an interval

$$ d\in [l,r],\qquad r=\left\lfloor\frac{N}{q}\right\rfloor. $$

So instead of handling each \(d\) in isolation, we process one maximal interval at a time. Over such an interval the superfactorial term and the exponent \(\binom{q}{2}\) are shared. Only three aggregated quantities are needed:

$$ P_{+}=\prod_{\substack{d=l\\ \mu(d)=1}}^{r} d,\qquad P_{-}=\prod_{\substack{d=l\\ \mu(d)=-1}}^{r} d,\qquad S_{\mu}=\sum_{d=l}^{r}\mu(d). $$

The interval contribution is then

$$ P_{+}^{\binom{q}{2}} \cdot P_{-}^{-\binom{q}{2}} \cdot \operatorname{SF}(q-1)^{S_{\mu}}. $$

Because the modulus \(10^9+7\) is prime and \(N<10^9+7\), every base appearing here is invertible modulo the modulus. Therefore negative Möbius exponents can be handled with modular inverses, or equivalently with exponents reduced modulo \(10^9+6\).

Worked Example: \(N=10\)

For \(N=10\), the quotient plateaus are

$$ \left\lfloor\frac{10}{d}\right\rfloor= \begin{cases} 10,& d=1,\\ 5,& d=2,\\ 3,& d=3,\\ 2,& d=4,5,\\ 1,& d=6,7,8,9,10. \end{cases} $$

The Möbius values are

$$ \mu(1)=1,\ \mu(2)=-1,\ \mu(3)=-1,\ \mu(4)=0,\ \mu(5)=-1,\ \mu(6)=1,\ \mu(7)=-1,\ \mu(8)=0,\ \mu(9)=0,\ \mu(10)=1. $$

All terms with \(q=1\) contribute \(d^0\operatorname{SF}(0)=1\), and every term with \(\mu(d)=0\) disappears. So

$$ P(10)=\frac{\operatorname{SF}(9)}{2^{10}\operatorname{SF}(4)\cdot 3^3\operatorname{SF}(2)\cdot 5}. $$

Numerically,

$$ \operatorname{SF}(9)=1!\,2!\,3!\cdots 9!=1834933472251084800000, $$

and the denominator equals

$$ 2^{10}\operatorname{SF}(4)\cdot 3^3\operatorname{SF}(2)\cdot 5 = 79626240. $$

Hence

$$ P(10)=23044331520000, $$

and modulo \(10^9+7\) this becomes

$$ P(10)\equiv 331358692. $$

This matches the checkpoint used by the implementation.

How the Code Works

The C++, Python, and Java implementations follow the same structure. First they compute the Möbius function \(\mu(1),\dots,\mu(N)\) with a linear sieve, so each integer is classified as squarefree with sign \(\pm1\) or discarded when \(\mu(d)=0\).

Next they enumerate the maximal intervals on which \(q=\lfloor N/d\rfloor\) is constant, and collect the corresponding keys \(q-1\). They then build factorials and superfactorials incrementally from \(1\) upward, storing the superfactorial values only at the keys that will actually be queried later.

For each quotient interval, the implementation multiplies together the \(d\)-values with positive Möbius sign, multiplies together the \(d\)-values with negative Möbius sign, and also sums the Möbius values over that interval. With those three aggregated numbers, one interval update reproduces the entire contribution

$$ P_{+}^{\binom{q}{2}} \cdot P_{-}^{-\binom{q}{2}} \cdot \operatorname{SF}(q-1)^{S_{\mu}} $$

modulo \(10^9+7\). Fast binary exponentiation is used for every modular power, and exponents are reduced modulo \(10^9+6\) using Fermat's little theorem.

Complexity Analysis

The linear sieve for \(\mu\) costs \(O(N)\) time and \(O(N)\) memory. Building factorial and superfactorial values up to \(N-1\) is another \(O(N)\) pass. The quotient-interval traversal visits each \(d\) exactly once overall, with only \(O(\sqrt{N})\) distinct plateaus and logarithmic-time modular exponentiations per plateau. Therefore the total running time is \(O(N)\), and the memory usage is \(O(N)\).

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=754
  2. Möbius function: Wikipedia — Möbius function
  3. Möbius inversion formula: Wikipedia — Möbius inversion formula
  4. Euler's totient function: Wikipedia — Euler's totient function
  5. Superfactorial: Wikipedia — Superfactorial
  6. Linear sieve: cp-algorithms — Linear Sieve

Problemzusammenfassung

Für jede positive ganze Zahl \(n\) sei das Gauß-Fakultätsprodukt definiert durch

$$\operatorname{GF}(n)=\prod_{\substack{1\le k\lt n\\ \gcd(k,n)=1}} k,$$

mit \(\operatorname{GF}(1)=1\). Gesucht ist

$$P(N)=\prod_{n=1}^{N}\operatorname{GF}(n)\pmod{10^9+7}$$

für \(N=10^8\). Ein direktes Durchlaufen aller teilerfremden Paare wäre viel zu langsam, daher formt die Lösung das gesamte Produkt in eine divisorbasierte Formel um, die sich linear auswerten lässt.

Mathematischer Ansatz

Die Grundidee besteht darin, jedes Gauß-Fakultätsprodukt per Möbius-Inversion in ein Teilerprodukt zu verwandeln und danach alle \(n\le N\) in Blöcken mit konstantem \(\left\lfloor N/d\right\rfloor\) zusammenzufassen.

Schritt 1: Die Faktoren von \((n-1)!\) nach ihrem ggT mit \(n\) gruppieren

Fixiere \(n\). Jede Zahl \(k\) mit \(1\le k\lt n\) besitzt einen ggT \(d=\gcd(k,n)\) und kann geschrieben werden als

$$k=d\,m,\qquad \gcd\!\left(m,\frac{n}{d}\right)=1,\qquad 1\le m\lt\frac{n}{d}.$$

Für einen festen Teiler \(d\mid n\) ist das Produkt aller solchen Zahlen gleich

$$d^{\phi(n/d)}\operatorname{GF}\!\left(\frac{n}{d}\right).$$

Über alle möglichen ggT-Werte multipliziert ergibt das

$$ (n-1)! = \prod_{\substack{d \mid n\\ d\lt n}} d^{\phi(n/d)}\operatorname{GF}\!\left(\frac{n}{d}\right). $$

Mit der Substitution \(m=n/d\) folgt

$$ (n-1)! = \prod_{\substack{m \mid n\\ m\gt 1}} \left(\frac{n}{m}\right)^{\phi(m)}\operatorname{GF}(m). $$

Schritt 2: Normalisieren und Möbius-Inversion anwenden

Da

$$\sum_{\substack{m \mid n\\ m\gt 1}}\phi(m)=n-1,$$

können wir durch \(n^{n-1}\) teilen und erhalten

$$ \frac{(n-1)!}{n^{n-1}}=\prod_{m \mid n}\frac{\operatorname{GF}(m)}{m^{\phi(m)}}. $$

Setze

$$ A(n)=\frac{(n-1)!}{n^{n-1}},\qquad H(n)=\frac{\operatorname{GF}(n)}{n^{\phi(n)}}. $$

Dann gilt \(A(n)=\prod_{d\mid n}H(d)\). Die multiplikative Möbius-Inversion liefert daher

$$ H(n)=\prod_{d \mid n} A(d)^{\mu(n/d)} =\prod_{d \mid n}\left(\frac{(d-1)!}{d^{d-1}}\right)^{\mu(n/d)}. $$

Nach dem Zurückmultiplizieren von \(n^{\phi(n)}\) und dem Vereinfachen der Potenzen von \(n\) entsteht die für die Implementierung bequeme Form

$$ \operatorname{GF}(n)=\prod_{d \mid n}\left(\left(\frac{n}{d}-1\right)! \, d^{\,n/d-1}\right)^{\mu(d)}. $$

Auch für \(n=1\) stimmt die Formel, denn dann ist der einzige Teiler \(d=1\) und der Faktor ist \(0!\cdot 1^0=1\).

Schritt 3: Über alle \(n\le N\) multiplizieren und die Reihenfolge tauschen

Nun multiplizieren wir die vorige Identität für alle \(1\le n\le N\):

$$ P(N)=\prod_{n=1}^{N}\prod_{d \mid n}\left(\left(\frac{n}{d}-1\right)! \, d^{\,n/d-1}\right)^{\mu(d)}. $$

Mit \(n=d\,m\) wird daraus

$$ P(N)=\prod_{d=1}^{N}\prod_{m=1}^{\lfloor N/d\rfloor}\left((m-1)! \, d^{\,m-1}\right)^{\mu(d)}. $$

Für festes \(d\) sei \(q=\left\lfloor N/d\right\rfloor\). Dann ist der Exponent von \(d\)

$$ \sum_{m=1}^{q}(m-1)=\frac{q(q-1)}{2}=\binom{q}{2}, $$

und der Fakultätsteil lautet

$$ \prod_{m=1}^{q}(m-1)! = \prod_{t=1}^{q-1} t!. $$

Schritt 4: Die Superfakultät einführen

Definiere die Superfakultät durch

$$ \operatorname{SF}(k)=\prod_{t=1}^{k} t!,\qquad \operatorname{SF}(0)=1. $$

Dann verdichtet sich das gesamte Produkt zu

$$ \boxed{ P(N)\equiv \prod_{d=1}^{N}\left(d^{\binom{\lfloor N/d\rfloor}{2}}\operatorname{SF}\!\left(\lfloor N/d\rfloor-1\right)\right)^{\mu(d)} \pmod{10^9+7}. } $$

Genau diese Formel setzen die Programme um.

Schritt 5: Gleiche Quotienten gemeinsam behandeln

Die Größe

$$q=\left\lfloor \frac{N}{d}\right\rfloor$$

ist auf einem Intervall

$$ d\in [l,r],\qquad r=\left\lfloor\frac{N}{q}\right\rfloor, $$

konstant. Daher verarbeitet man jeweils ein maximales Intervall statt jedes \(d\) einzeln. Auf einem solchen Block sind sowohl \(\binom{q}{2}\) als auch \(\operatorname{SF}(q-1)\) gemeinsam. Benötigt werden nur drei aggregierte Größen:

$$ P_{+}=\prod_{\substack{d=l\\ \mu(d)=1}}^{r} d,\qquad P_{-}=\prod_{\substack{d=l\\ \mu(d)=-1}}^{r} d,\qquad S_{\mu}=\sum_{d=l}^{r}\mu(d). $$

Der Beitrag des Intervalls ist dann

$$ P_{+}^{\binom{q}{2}} \cdot P_{-}^{-\binom{q}{2}} \cdot \operatorname{SF}(q-1)^{S_{\mu}}. $$

Da der Modul \(10^9+7\) prim ist und \(N<10^9+7\) gilt, ist jede auftretende Basis modulo dem Modul invertierbar. Negative Möbius-Exponenten werden also über modulare Inversen beziehungsweise über Exponenten modulo \(10^9+6\) behandelt.

Durchgerechnetes Beispiel: \(N=10\)

Für \(N=10\) sind die Quotientenplateaus

$$ \left\lfloor\frac{10}{d}\right\rfloor= \begin{cases} 10,& d=1,\\ 5,& d=2,\\ 3,& d=3,\\ 2,& d=4,5,\\ 1,& d=6,7,8,9,10. \end{cases} $$

Die Möbius-Werte lauten

$$ \mu(1)=1,\ \mu(2)=-1,\ \mu(3)=-1,\ \mu(4)=0,\ \mu(5)=-1,\ \mu(6)=1,\ \mu(7)=-1,\ \mu(8)=0,\ \mu(9)=0,\ \mu(10)=1. $$

Alle Terme mit \(q=1\) liefern \(d^0\operatorname{SF}(0)=1\), und jede Stelle mit \(\mu(d)=0\) fällt weg. Somit gilt

$$ P(10)=\frac{\operatorname{SF}(9)}{2^{10}\operatorname{SF}(4)\cdot 3^3\operatorname{SF}(2)\cdot 5}. $$

Numerisch ist

$$ \operatorname{SF}(9)=1!\,2!\,3!\cdots 9!=1834933472251084800000, $$

und der Nenner ist

$$ 2^{10}\operatorname{SF}(4)\cdot 3^3\operatorname{SF}(2)\cdot 5 = 79626240. $$

Daher folgt

$$ P(10)=23044331520000, $$

und modulo \(10^9+7\) erhält man

$$ P(10)\equiv 331358692. $$

Genau dieser Kontrollwert wird von der Implementierung bestätigt.

Wie der Code arbeitet

Die Implementierungen in C++, Python und Java folgen derselben Struktur. Zuerst berechnen sie mit einem linearen Sieb die Möbius-Funktion \(\mu(1),\dots,\mu(N)\), sodass jede Zahl als quadratfrei mit Vorzeichen \(\pm1\) klassifiziert oder bei \(\mu(d)=0\) verworfen wird.

Danach werden die maximalen Intervalle bestimmt, auf denen \(q=\lfloor N/d\rfloor\) konstant ist, und die dazugehörigen Schlüssel \(q-1\) gesammelt. Anschließend bauen die Programme Fakultäten und Superfakultäten inkrementell ab \(1\) auf und speichern die Superfakultätswerte nur an den später tatsächlich benötigten Stellen.

Für jedes Quotientenintervall werden die \(d\)-Werte mit positivem Möbius-Vorzeichen multipliziert, die \(d\)-Werte mit negativem Vorzeichen ebenfalls multipliziert, und zugleich die Möbius-Werte über das Intervall summiert. Mit diesen drei aggregierten Größen reproduziert ein einziges Update den gesamten Beitrag

$$ P_{+}^{\binom{q}{2}} \cdot P_{-}^{-\binom{q}{2}} \cdot \operatorname{SF}(q-1)^{S_{\mu}} $$

modulo \(10^9+7\). Für alle modularen Potenzen wird schnelle binäre Exponentiation verwendet; die Exponenten werden per kleinem Satz von Fermat modulo \(10^9+6\) reduziert.

Komplexitätsanalyse

Das lineare Sieb für \(\mu\) kostet \(O(N)\) Zeit und \(O(N)\) Speicher. Der Aufbau von Fakultäten und Superfakultäten bis \(N-1\) ist ein weiterer \(O(N)\)-Durchlauf. Die Quotientenblock-Verarbeitung besucht jedes \(d\) insgesamt genau einmal; dazu kommen nur \(O(\sqrt{N})\) verschiedene Plateaus und logarithmische Kosten für modulare Potenzen pro Plateau. Insgesamt ergibt sich also Laufzeit \(O(N)\) und Speicherbedarf \(O(N)\).

Fußnoten und Referenzen

  1. Problemseite: https://projecteuler.net/problem=754
  2. Möbius-Funktion: Wikipedia — Möbius function
  3. Möbius-Inversionsformel: Wikipedia — Möbius inversion formula
  4. Eulersche Phi-Funktion: Wikipedia — Euler's totient function
  5. Superfakultät: Wikipedia — Superfactorial
  6. Lineares Sieb: cp-algorithms — Linear Sieve

Problem Özeti

Her pozitif \(n\) için Gauss faktöriyeli

$$\operatorname{GF}(n)=\prod_{\substack{1\le k\lt n\\ \gcd(k,n)=1}} k$$

olarak tanımlansın; ayrıca \(\operatorname{GF}(1)=1\). İstenen büyüklük

$$P(N)=\prod_{n=1}^{N}\operatorname{GF}(n)\pmod{10^9+7}$$

ve problemde \(N=10^8\) alınır. Doğrudan yaklaşım, tüm aralarında asal çiftleri tek tek dolaşmak zorunda kalacağı için pratik değildir. Çözüm bunun yerine tüm çarpımı bölenler üzerinden yeniden yazar ve hesabı lineer zamana indirir.

Matematiksel Yaklaşım

Ana fikir, her bir Gauss faktöriyelini Möbius terslemesiyle bölen bazlı bir çarpıma dönüştürmek ve sonra tüm \(n\le N\) terimlerini \(\left\lfloor N/d\right\rfloor\) sabit kalan bloklarda toplamaktır.

Adım 1: \((n-1)!\) içindeki terimleri \(n\) ile EBOB'larına göre grupla

Sabit bir \(n\) seçelim. \(1\le k\lt n\) aralığındaki her sayı için \(d=\gcd(k,n)\) yazabiliriz ve bu durumda

$$k=d\,m,\qquad \gcd\!\left(m,\frac{n}{d}\right)=1,\qquad 1\le m\lt\frac{n}{d}$$

olur. Belirli bir \(d\mid n\) böleni için bu sınıftaki tüm sayıların çarpımı

$$d^{\phi(n/d)}\operatorname{GF}\!\left(\frac{n}{d}\right)$$

şeklindedir. Bütün olası EBOB sınıflarını çarptığımızda

$$ (n-1)! = \prod_{\substack{d \mid n\\ d\lt n}} d^{\phi(n/d)}\operatorname{GF}\!\left(\frac{n}{d}\right) $$

elde edilir. \(m=n/d\) değişkeniyle bu ifade

$$ (n-1)! = \prod_{\substack{m \mid n\\ m\gt 1}} \left(\frac{n}{m}\right)^{\phi(m)}\operatorname{GF}(m) $$

biçimini alır.

Adım 2: Normalize et ve Möbius terslemesini uygula

Şu toplam bilinir:

$$\sum_{\substack{m \mid n\\ m\gt 1}}\phi(m)=n-1.$$

Bu nedenle her iki tarafı \(n^{n-1}\)'e bölersek

$$ \frac{(n-1)!}{n^{n-1}}=\prod_{m \mid n}\frac{\operatorname{GF}(m)}{m^{\phi(m)}} $$

olur. Şimdi

$$ A(n)=\frac{(n-1)!}{n^{n-1}},\qquad H(n)=\frac{\operatorname{GF}(n)}{n^{\phi(n)}} $$

tanımlarını yapalım. O zaman \(A(n)=\prod_{d\mid n}H(d)\) ve çarpımsal Möbius terslemesi

$$ H(n)=\prod_{d \mid n} A(d)^{\mu(n/d)} =\prod_{d \mid n}\left(\frac{(d-1)!}{d^{d-1}}\right)^{\mu(n/d)} $$

sonucunu verir. \(n^{\phi(n)}\) tekrar eklenip \(n\)'nin üsleri sadeleştirildiğinde kod için daha kullanışlı olan

$$ \operatorname{GF}(n)=\prod_{d \mid n}\left(\left(\frac{n}{d}-1\right)! \, d^{\,n/d-1}\right)^{\mu(d)} $$

kimliğini elde ederiz. \(n=1\) için de tek bölen \(d=1\) olduğundan faktör \(0!\cdot 1^0=1\) çıkar.

Adım 3: Tüm \(n\le N\) için çarp ve sırayı değiştir

Önceki formülü \(1\le n\le N\) için çarparsak

$$ P(N)=\prod_{n=1}^{N}\prod_{d \mid n}\left(\left(\frac{n}{d}-1\right)! \, d^{\,n/d-1}\right)^{\mu(d)} $$

olur. \(n=d\,m\) yazarsak dış çarpım

$$ P(N)=\prod_{d=1}^{N}\prod_{m=1}^{\lfloor N/d\rfloor}\left((m-1)! \, d^{\,m-1}\right)^{\mu(d)} $$

şeklinde düzenlenir. Sabit bir \(d\) için \(q=\left\lfloor N/d\right\rfloor\) olsun. O zaman \(d\)'nin toplam üssü

$$ \sum_{m=1}^{q}(m-1)=\frac{q(q-1)}{2}=\binom{q}{2} $$

ve faktöriyel kısmı

$$ \prod_{m=1}^{q}(m-1)! = \prod_{t=1}^{q-1} t! $$

olur.

Adım 4: Süperfaktöriyeli tanımla

Süperfaktöriyeli

$$ \operatorname{SF}(k)=\prod_{t=1}^{k} t!,\qquad \operatorname{SF}(0)=1 $$

olarak tanımlayalım. Böylece tüm ifade

$$ \boxed{ P(N)\equiv \prod_{d=1}^{N}\left(d^{\binom{\lfloor N/d\rfloor}{2}}\operatorname{SF}\!\left(\lfloor N/d\rfloor-1\right)\right)^{\mu(d)} \pmod{10^9+7} } $$

haline gelir. Uygulamanın hesapladığı ana formül tam olarak budur.

Adım 5: Aynı bölüm değerlerini blok halinde işle

\(q=\left\lfloor N/d\right\rfloor\) değeri

$$ d\in [l,r],\qquad r=\left\lfloor\frac{N}{q}\right\rfloor $$

aralığında sabittir. Bu yüzden her \(d\) için ayrı ayrı uğraşmak yerine maksimum bloklar dolaşılır. Böyle bir blokta hem \(\binom{q}{2}\) hem de \(\operatorname{SF}(q-1)\) ortaktır. Gerekli üç toplulaştırma şunlardır:

$$ P_{+}=\prod_{\substack{d=l\\ \mu(d)=1}}^{r} d,\qquad P_{-}=\prod_{\substack{d=l\\ \mu(d)=-1}}^{r} d,\qquad S_{\mu}=\sum_{d=l}^{r}\mu(d). $$

Blok katkısı

$$ P_{+}^{\binom{q}{2}} \cdot P_{-}^{-\binom{q}{2}} \cdot \operatorname{SF}(q-1)^{S_{\mu}} $$

şeklindedir. Modül \(10^9+7\) asaldır ve \(N<10^9+7\) olduğundan ortaya çıkan tüm tabanlar modül altında terslenebilir. Böylece negatif Möbius üsleri ya modüler terslerle ya da üsleri \(10^9+6\) modunda indirerek güvenle uygulanır.

Çalışılmış Örnek: \(N=10\)

\(N=10\) için bölüm platoları

$$ \left\lfloor\frac{10}{d}\right\rfloor= \begin{cases} 10,& d=1,\\ 5,& d=2,\\ 3,& d=3,\\ 2,& d=4,5,\\ 1,& d=6,7,8,9,10 \end{cases} $$

şeklindedir. Möbius değerleri ise

$$ \mu(1)=1,\ \mu(2)=-1,\ \mu(3)=-1,\ \mu(4)=0,\ \mu(5)=-1,\ \mu(6)=1,\ \mu(7)=-1,\ \mu(8)=0,\ \mu(9)=0,\ \mu(10)=1 $$

olur. \(q=1\) olan tüm terimler \(d^0\operatorname{SF}(0)=1\) verdiği için etkisizdir; \(\mu(d)=0\) olanlar da zaten düşer. Dolayısıyla

$$ P(10)=\frac{\operatorname{SF}(9)}{2^{10}\operatorname{SF}(4)\cdot 3^3\operatorname{SF}(2)\cdot 5} $$

elde edilir. Sayısal olarak

$$ \operatorname{SF}(9)=1!\,2!\,3!\cdots 9!=1834933472251084800000 $$

ve payda

$$ 2^{10}\operatorname{SF}(4)\cdot 3^3\operatorname{SF}(2)\cdot 5 = 79626240 $$

olduğundan

$$ P(10)=23044331520000 $$

çıkar. Bunun \(10^9+7\) modundaki değeri

$$ P(10)\equiv 331358692 $$

olup uygulamadaki kontrol sonucu ile aynıdır.

Kod Nasıl Çalışır

C++, Python ve Java uygulamaları aynı iskeleti izler. Önce lineer elek ile \(\mu(1),\dots,\mu(N)\) hesaplanır; böylece her sayı karefree olup olmadığına ve işaretine göre sınıflanır, \(\mu(d)=0\) olanlar baştan elenir.

Daha sonra \(q=\lfloor N/d\rfloor\) sabit kalan maksimum aralıklar çıkarılır ve ilgili \(q-1\) anahtarları toplanır. Ardından faktöriyel ve süperfaktöriyel değerleri \(1\)'den başlayarak artımlı biçimde üretilir; süperfaktöriyel yalnızca ileride gerçekten sorgulanacak noktalarda saklanır.

Her bölüm bloğu için, Möbius işareti pozitif olan \(d\)'ler çarpılır, negatif olan \(d\)'ler ayrı çarpılır ve aynı bloktaki Möbius değerleri toplanır. Bu üç toplu bilgi ile blok katkısı

$$ P_{+}^{\binom{q}{2}} \cdot P_{-}^{-\binom{q}{2}} \cdot \operatorname{SF}(q-1)^{S_{\mu}} $$

tek güncellemede uygulanır. Tüm modüler üs alma işlemleri hızlı ikili üs alma ile yapılır; üsler de Fermat sayesinde \(10^9+6\) modunda indirgenir.

Karmaşıklık Analizi

\(\mu\) için lineer elek \(O(N)\) zaman ve \(O(N)\) bellek kullanır. Faktöriyel ve süperfaktöriyellerin \(N-1\)'e kadar hazırlanması da ayrı bir \(O(N)\) geçiştir. Bölüm bloklarının işlenmesinde her \(d\) toplamda tam bir kez ziyaret edilir; ayrıca yalnızca \(O(\sqrt{N})\) farklı plato vardır ve her plato için logaritmik maliyetli modüler üs alma yapılır. Bu nedenle toplam çalışma süresi \(O(N)\), bellek kullanımı da \(O(N)\) olur.

Dipnotlar ve Referanslar

  1. Problem sayfası: https://projecteuler.net/problem=754
  2. Möbius fonksiyonu: Wikipedia — Möbius function
  3. Möbius tersleme formülü: Wikipedia — Möbius inversion formula
  4. Euler totient fonksiyonu: Wikipedia — Euler's totient function
  5. Süperfaktöriyel: Wikipedia — Superfactorial
  6. Lineer elek: cp-algorithms — Linear Sieve

Resumen del Problema

Para cada entero positivo \(n\), definimos el factorial de Gauss por

$$\operatorname{GF}(n)=\prod_{\substack{1\le k\lt n\\ \gcd(k,n)=1}} k,$$

con \(\operatorname{GF}(1)=1\). El objetivo es calcular

$$P(N)=\prod_{n=1}^{N}\operatorname{GF}(n)\pmod{10^9+7}$$

cuando \(N=10^8\). Una enumeración directa de todos los pares coprimos sería demasiado costosa, así que la solución reescribe el producto completo como una expresión sobre divisores y luego la evalúa en tiempo lineal.

Enfoque Matemático

La idea central es convertir cada factorial de Gauss en un producto sobre divisores mediante inversión de Möbius, y después reagrupar todos los términos según los bloques donde \(\left\lfloor N/d\right\rfloor\) permanece constante.

Paso 1: Agrupar los factores de \((n-1)!\) por su mcd con \(n\)

Fijemos \(n\). Cada entero \(k\) con \(1\le k\lt n\) tiene algún divisor común \(d=\gcd(k,n)\), de modo que puede escribirse como

$$k=d\,m,\qquad \gcd\!\left(m,\frac{n}{d}\right)=1,\qquad 1\le m\lt\frac{n}{d}.$$

Para un divisor fijo \(d\mid n\), el producto de todos los términos de esa clase es

$$d^{\phi(n/d)}\operatorname{GF}\!\left(\frac{n}{d}\right).$$

Multiplicando sobre todos los posibles valores del mcd obtenemos

$$ (n-1)! = \prod_{\substack{d \mid n\\ d\lt n}} d^{\phi(n/d)}\operatorname{GF}\!\left(\frac{n}{d}\right). $$

Si hacemos el cambio \(m=n/d\), esto se convierte en

$$ (n-1)! = \prod_{\substack{m \mid n\\ m\gt 1}} \left(\frac{n}{m}\right)^{\phi(m)}\operatorname{GF}(m). $$

Paso 2: Normalizar y aplicar inversión de Möbius

Como

$$\sum_{\substack{m \mid n\\ m\gt 1}}\phi(m)=n-1,$$

podemos dividir por \(n^{n-1}\) y reescribir la identidad como

$$ \frac{(n-1)!}{n^{n-1}}=\prod_{m \mid n}\frac{\operatorname{GF}(m)}{m^{\phi(m)}}. $$

Definimos

$$ A(n)=\frac{(n-1)!}{n^{n-1}},\qquad H(n)=\frac{\operatorname{GF}(n)}{n^{\phi(n)}}. $$

Entonces \(A(n)=\prod_{d\mid n}H(d)\), y la forma multiplicativa de la inversión de Möbius da

$$ H(n)=\prod_{d \mid n} A(d)^{\mu(n/d)} =\prod_{d \mid n}\left(\frac{(d-1)!}{d^{d-1}}\right)^{\mu(n/d)}. $$

Al volver a multiplicar por \(n^{\phi(n)}\) y simplificar las potencias de \(n\), llegamos a la forma que usan las implementaciones:

$$ \operatorname{GF}(n)=\prod_{d \mid n}\left(\left(\frac{n}{d}-1\right)! \, d^{\,n/d-1}\right)^{\mu(d)}. $$

También funciona para \(n=1\), porque el único divisor es \(d=1\) y el factor resultante es \(0!\cdot 1^0=1\).

Paso 3: Multiplicar sobre todos los \(n\le N\) e intercambiar el orden

Multiplicamos la identidad anterior para cada \(1\le n\le N\):

$$ P(N)=\prod_{n=1}^{N}\prod_{d \mid n}\left(\left(\frac{n}{d}-1\right)! \, d^{\,n/d-1}\right)^{\mu(d)}. $$

Escribiendo \(n=d\,m\), el producto total pasa a ser

$$ P(N)=\prod_{d=1}^{N}\prod_{m=1}^{\lfloor N/d\rfloor}\left((m-1)! \, d^{\,m-1}\right)^{\mu(d)}. $$

Para un \(d\) fijo, sea \(q=\left\lfloor N/d\right\rfloor\). Entonces el exponente total de \(d\) es

$$ \sum_{m=1}^{q}(m-1)=\frac{q(q-1)}{2}=\binom{q}{2}, $$

y la parte factorial es

$$ \prod_{m=1}^{q}(m-1)! = \prod_{t=1}^{q-1} t!. $$

Paso 4: Introducir el superfactorial

Definimos

$$ \operatorname{SF}(k)=\prod_{t=1}^{k} t!,\qquad \operatorname{SF}(0)=1. $$

Con esa notación, todo el producto se resume en

$$ \boxed{ P(N)\equiv \prod_{d=1}^{N}\left(d^{\binom{\lfloor N/d\rfloor}{2}}\operatorname{SF}\!\left(\lfloor N/d\rfloor-1\right)\right)^{\mu(d)} \pmod{10^9+7}. } $$

Ésta es exactamente la identidad que evalúa el algoritmo.

Paso 5: Agrupar los cocientes iguales

La cantidad

$$q=\left\lfloor \frac{N}{d}\right\rfloor$$

permanece constante en un intervalo

$$ d\in [l,r],\qquad r=\left\lfloor\frac{N}{q}\right\rfloor. $$

Por eso conviene procesar cada intervalo máximo de una sola vez. En un bloque así, tanto \(\binom{q}{2}\) como \(\operatorname{SF}(q-1)\) son comunes. Sólo hacen falta tres agregados:

$$ P_{+}=\prod_{\substack{d=l\\ \mu(d)=1}}^{r} d,\qquad P_{-}=\prod_{\substack{d=l\\ \mu(d)=-1}}^{r} d,\qquad S_{\mu}=\sum_{d=l}^{r}\mu(d). $$

La contribución del bloque es entonces

$$ P_{+}^{\binom{q}{2}} \cdot P_{-}^{-\binom{q}{2}} \cdot \operatorname{SF}(q-1)^{S_{\mu}}. $$

Como el módulo \(10^9+7\) es primo y \(N<10^9+7\), todas las bases que aparecen son invertibles módulo el primo. Por eso los exponentes negativos de Möbius pueden manejarse con inversos modulares, o de forma equivalente reduciendo exponentes módulo \(10^9+6\).

Ejemplo Desarrollado: \(N=10\)

Para \(N=10\), los bloques de cociente son

$$ \left\lfloor\frac{10}{d}\right\rfloor= \begin{cases} 10,& d=1,\\ 5,& d=2,\\ 3,& d=3,\\ 2,& d=4,5,\\ 1,& d=6,7,8,9,10. \end{cases} $$

Los valores de Möbius son

$$ \mu(1)=1,\ \mu(2)=-1,\ \mu(3)=-1,\ \mu(4)=0,\ \mu(5)=-1,\ \mu(6)=1,\ \mu(7)=-1,\ \mu(8)=0,\ \mu(9)=0,\ \mu(10)=1. $$

Todos los términos con \(q=1\) aportan \(d^0\operatorname{SF}(0)=1\), y los términos con \(\mu(d)=0\) desaparecen. Por tanto,

$$ P(10)=\frac{\operatorname{SF}(9)}{2^{10}\operatorname{SF}(4)\cdot 3^3\operatorname{SF}(2)\cdot 5}. $$

Numéricamente,

$$ \operatorname{SF}(9)=1!\,2!\,3!\cdots 9!=1834933472251084800000, $$

mientras que el denominador es

$$ 2^{10}\operatorname{SF}(4)\cdot 3^3\operatorname{SF}(2)\cdot 5 = 79626240. $$

Así,

$$ P(10)=23044331520000, $$

y reduciendo módulo \(10^9+7\) obtenemos

$$ P(10)\equiv 331358692. $$

Ese valor coincide con el punto de verificación usado por la implementación.

Cómo Funciona el Código

Las implementaciones en C++, Python y Java siguen la misma organización. Primero calculan la función de Möbius \(\mu(1),\dots,\mu(N)\) con una criba lineal, de modo que cada entero queda clasificado como libre de cuadrados con signo \(\pm1\), o bien descartado cuando \(\mu(d)=0\).

Luego enumeran los intervalos máximos en los que \(q=\lfloor N/d\rfloor\) es constante y recopilan las claves \(q-1\) asociadas. Después construyen factoriales y superfactoriales de forma incremental desde \(1\), guardando los valores de superfactorial sólo en los puntos que realmente se usarán más adelante.

En cada bloque de cociente, la implementación multiplica los \(d\) con signo de Möbius positivo, multiplica por separado los \(d\) con signo negativo y además suma los valores de Möbius del bloque. Con esos tres acumulados realiza una sola actualización que reproduce

$$ P_{+}^{\binom{q}{2}} \cdot P_{-}^{-\binom{q}{2}} \cdot \operatorname{SF}(q-1)^{S_{\mu}} $$

módulo \(10^9+7\). Todas las potencias modulares se calculan con exponenciación binaria rápida, y los exponentes se reducen módulo \(10^9+6\) usando el pequeño teorema de Fermat.

Análisis de Complejidad

La criba lineal para \(\mu\) cuesta \(O(N)\) tiempo y \(O(N)\) memoria. Construir factoriales y superfactoriales hasta \(N-1\) requiere otro recorrido \(O(N)\). El procesamiento por bloques de cociente visita cada \(d\) exactamente una vez en total; además sólo hay \(O(\sqrt{N})\) mesetas distintas y un número logarítmico de multiplicaciones por potencia modular en cada una. Por ello, el tiempo total es \(O(N)\) y el uso de memoria es \(O(N)\).

Notas y Referencias

  1. Página del problema: https://projecteuler.net/problem=754
  2. Función de Möbius: Wikipedia — Möbius function
  3. Fórmula de inversión de Möbius: Wikipedia — Möbius inversion formula
  4. Función phi de Euler: Wikipedia — Euler's totient function
  5. Superfactorial: Wikipedia — Superfactorial
  6. Criba lineal: cp-algorithms — Linear Sieve

问题概述

对每个正整数 \(n\),定义 Gauss 阶乘为

$$\operatorname{GF}(n)=\prod_{\substack{1\le k\lt n\\ \gcd(k,n)=1}} k,$$

并约定 \(\operatorname{GF}(1)=1\)。题目要求计算

$$P(N)=\prod_{n=1}^{N}\operatorname{GF}(n)\pmod{10^9+7}$$

其中 \(N=10^8\)。如果直接枚举所有满足 \(\gcd(k,n)=1\) 的配对,规模会大得不可行。真正可行的做法是先把单个 Gauss 阶乘改写成与除数有关的乘积,再把整个总乘积压缩成可以线性处理的形式。

数学方法

核心思路分成两步:先对单个 \(\operatorname{GF}(n)\) 做 Möbius 反演,再把所有 \(n\le N\) 的贡献按照 \(\left\lfloor N/d\right\rfloor\) 相同的区间合并。

步骤 1:按 \(\gcd(k,n)\) 对 \((n-1)!\) 中的项分组

固定 \(n\)。任意满足 \(1\le k\lt n\) 的整数,都有某个公因子

$$d=\gcd(k,n).$$

因此它可以写成

$$k=d\,m,\qquad \gcd\!\left(m,\frac{n}{d}\right)=1,\qquad 1\le m\lt\frac{n}{d}.$$

对固定的除数 \(d\mid n\),这一类所有数的乘积就是

$$d^{\phi(n/d)}\operatorname{GF}\!\left(\frac{n}{d}\right).$$

把所有可能的 \(\gcd\) 类别乘起来,就得到

$$ (n-1)! = \prod_{\substack{d \mid n\\ d\lt n}} d^{\phi(n/d)}\operatorname{GF}\!\left(\frac{n}{d}\right). $$

把变量改成 \(m=n/d\) 以后,上式变成

$$ (n-1)! = \prod_{\substack{m \mid n\\ m\gt 1}} \left(\frac{n}{m}\right)^{\phi(m)}\operatorname{GF}(m). $$

步骤 2:归一化并应用 Möbius 反演

由于

$$\sum_{\substack{m \mid n\\ m\gt 1}}\phi(m)=n-1,$$

所以两边同时除以 \(n^{n-1}\) 可以写成

$$ \frac{(n-1)!}{n^{n-1}}=\prod_{m \mid n}\frac{\operatorname{GF}(m)}{m^{\phi(m)}}. $$

定义

$$ A(n)=\frac{(n-1)!}{n^{n-1}},\qquad H(n)=\frac{\operatorname{GF}(n)}{n^{\phi(n)}}. $$

于是有 \(A(n)=\prod_{d\mid n}H(d)\)。对这个乘法型关系使用 Möbius 反演,可得

$$ H(n)=\prod_{d \mid n} A(d)^{\mu(n/d)} =\prod_{d \mid n}\left(\frac{(d-1)!}{d^{d-1}}\right)^{\mu(n/d)}. $$

再乘回 \(n^{\phi(n)}\),并把 \(n\) 的幂整理掉,就得到更适合实现的形式:

$$ \operatorname{GF}(n)=\prod_{d \mid n}\left(\left(\frac{n}{d}-1\right)! \, d^{\,n/d-1}\right)^{\mu(d)}. $$

当 \(n=1\) 时也完全成立,因为唯一的除数是 \(d=1\),对应因子就是 \(0!\cdot 1^0=1\)。

步骤 3:对所有 \(n\le N\) 求积并交换乘积顺序

把上面的公式对每个 \(1\le n\le N\) 都乘起来:

$$ P(N)=\prod_{n=1}^{N}\prod_{d \mid n}\left(\left(\frac{n}{d}-1\right)! \, d^{\,n/d-1}\right)^{\mu(d)}. $$

令 \(n=d\,m\),就能把总乘积改写为

$$ P(N)=\prod_{d=1}^{N}\prod_{m=1}^{\lfloor N/d\rfloor}\left((m-1)! \, d^{\,m-1}\right)^{\mu(d)}. $$

对固定的 \(d\),记

$$q=\left\lfloor \frac{N}{d}\right\rfloor.$$

那么 \(d\) 的总指数是

$$ \sum_{m=1}^{q}(m-1)=\frac{q(q-1)}{2}=\binom{q}{2}, $$

而阶乘部分则是

$$ \prod_{m=1}^{q}(m-1)! = \prod_{t=1}^{q-1} t!. $$

步骤 4:引入超阶乘

定义超阶乘

$$ \operatorname{SF}(k)=\prod_{t=1}^{k} t!,\qquad \operatorname{SF}(0)=1. $$

这样总公式就压缩成

$$ \boxed{ P(N)\equiv \prod_{d=1}^{N}\left(d^{\binom{\lfloor N/d\rfloor}{2}}\operatorname{SF}\!\left(\lfloor N/d\rfloor-1\right)\right)^{\mu(d)} \pmod{10^9+7}. } $$

这就是实现真正计算的数学恒等式。

步骤 5:按相同的整除商分块

$$q=\left\lfloor \frac{N}{d}\right\rfloor$$

会在一个区间内保持不变。若当前商为 \(q\),那么这个最大区间是

$$ d\in [l,r],\qquad r=\left\lfloor\frac{N}{q}\right\rfloor. $$

因此没有必要逐个 \(d\) 重新计算相同的 \(\binom{q}{2}\) 和 \(\operatorname{SF}(q-1)\)。在同一个区间里,只需维护三个聚合量:

$$ P_{+}=\prod_{\substack{d=l\\ \mu(d)=1}}^{r} d,\qquad P_{-}=\prod_{\substack{d=l\\ \mu(d)=-1}}^{r} d,\qquad S_{\mu}=\sum_{d=l}^{r}\mu(d). $$

这样整个区间的贡献就变成

$$ P_{+}^{\binom{q}{2}} \cdot P_{-}^{-\binom{q}{2}} \cdot \operatorname{SF}(q-1)^{S_{\mu}}. $$

由于模数 \(10^9+7\) 是素数,且 \(N<10^9+7\),这里出现的每个底数在模意义下都可逆。所以 Möbius 产生的负指数可以通过模逆处理,也可以等价地把指数按 \(10^9+6\) 取模后再做快速幂。

计算示例:\(N=10\)

当 \(N=10\) 时,\(\left\lfloor 10/d\right\rfloor\) 的分块为

$$ \left\lfloor\frac{10}{d}\right\rfloor= \begin{cases} 10,& d=1,\\ 5,& d=2,\\ 3,& d=3,\\ 2,& d=4,5,\\ 1,& d=6,7,8,9,10. \end{cases} $$

相应的 Möbius 值是

$$ \mu(1)=1,\ \mu(2)=-1,\ \mu(3)=-1,\ \mu(4)=0,\ \mu(5)=-1,\ \mu(6)=1,\ \mu(7)=-1,\ \mu(8)=0,\ \mu(9)=0,\ \mu(10)=1. $$

所有 \(q=1\) 的项都只给出 \(d^0\operatorname{SF}(0)=1\),而 \(\mu(d)=0\) 的项直接消失,所以

$$ P(10)=\frac{\operatorname{SF}(9)}{2^{10}\operatorname{SF}(4)\cdot 3^3\operatorname{SF}(2)\cdot 5}. $$

数值上,

$$ \operatorname{SF}(9)=1!\,2!\,3!\cdots 9!=1834933472251084800000, $$

分母为

$$ 2^{10}\operatorname{SF}(4)\cdot 3^3\operatorname{SF}(2)\cdot 5 = 79626240. $$

因此

$$ P(10)=23044331520000, $$

再取模 \(10^9+7\) 后得到

$$ P(10)\equiv 331358692. $$

这与实现中的校验值完全一致。

代码如何工作

C++、Python 和 Java 实现遵循同一套结构。首先用线性筛计算 \(\mu(1),\dots,\mu(N)\),把每个整数分类为平方因子自由且符号为 \(\pm1\),或者在 \(\mu(d)=0\) 时直接忽略。

接着枚举所有使 \(q=\lfloor N/d\rfloor\) 保持常数的最大区间,并收集会用到的 \(q-1\) 键值。然后从 \(1\) 开始递推阶乘与超阶乘,只在后面真正要查询的点上保存超阶乘结果。

对每个商区间,实现会分别把 Möbius 值为 \(+1\) 的 \(d\) 相乘、把 Möbius 值为 \(-1\) 的 \(d\) 相乘,同时累计该区间的 Möbius 和。依靠这三个聚合量,就能一次性完成

$$ P_{+}^{\binom{q}{2}} \cdot P_{-}^{-\binom{q}{2}} \cdot \operatorname{SF}(q-1)^{S_{\mu}} $$

这一整块更新。所有模幂都通过快速二进制幂完成,指数则利用费马小定理按 \(10^9+6\) 约化。

复杂度分析

计算 \(\mu\) 的线性筛需要 \(O(N)\) 时间和 \(O(N)\) 内存。把阶乘和超阶乘顺次建到 \(N-1\) 还需要一次 \(O(N)\) 扫描。整除商分块阶段总体上对每个 \(d\) 只访问一次;不同的商区间只有 \(O(\sqrt{N})\) 个,每个区间只额外做对数复杂度的模幂运算。所以总时间复杂度为 \(O(N)\),空间复杂度也为 \(O(N)\)。

脚注与参考资料

  1. 题目页面: https://projecteuler.net/problem=754
  2. Möbius 函数: Wikipedia — Möbius function
  3. Möbius 反演公式: Wikipedia — Möbius inversion formula
  4. Euler φ 函数: Wikipedia — Euler's totient function
  5. 超阶乘: Wikipedia — Superfactorial
  6. 线性筛: cp-algorithms — Linear Sieve

Краткое описание задачи

Для каждого положительного целого \(n\) определим факториал Гаусса

$$\operatorname{GF}(n)=\prod_{\substack{1\le k\lt n\\ \gcd(k,n)=1}} k,$$

причем \(\operatorname{GF}(1)=1\). Требуется вычислить

$$P(N)=\prod_{n=1}^{N}\operatorname{GF}(n)\pmod{10^9+7}$$

при \(N=10^8\). Прямой перебор всех взаимно простых пар слишком дорог, поэтому решение сначала переписывает отдельный факториал Гаусса в виде произведения по делителям, а затем сжимает общий ответ до линейного по времени алгоритма.

Математический подход

Основная идея такова: сначала применить инверсию Мёбиуса к формуле для одного \(\operatorname{GF}(n)\), а затем сгруппировать все вклады по интервалам, на которых \(\left\lfloor N/d\right\rfloor\) не меняется.

Шаг 1: Разбить множители \((n-1)!\) по их НОД с \(n\)

Зафиксируем \(n\). Любое число \(k\) при \(1\le k\lt n\) имеет некоторый НОД

$$d=\gcd(k,n),$$

поэтому его можно представить в виде

$$k=d\,m,\qquad \gcd\!\left(m,\frac{n}{d}\right)=1,\qquad 1\le m\lt\frac{n}{d}.$$

Для фиксированного делителя \(d\mid n\) произведение всех чисел из этого класса равно

$$d^{\phi(n/d)}\operatorname{GF}\!\left(\frac{n}{d}\right).$$

Перемножая по всем возможным значениям НОД, получаем

$$ (n-1)! = \prod_{\substack{d \mid n\\ d\lt n}} d^{\phi(n/d)}\operatorname{GF}\!\left(\frac{n}{d}\right). $$

После замены переменной \(m=n/d\) это становится

$$ (n-1)! = \prod_{\substack{m \mid n\\ m\gt 1}} \left(\frac{n}{m}\right)^{\phi(m)}\operatorname{GF}(m). $$

Шаг 2: Нормализация и инверсия Мёбиуса

Так как

$$\sum_{\substack{m \mid n\\ m\gt 1}}\phi(m)=n-1,$$

можно разделить обе части на \(n^{n-1}\) и переписать равенство так:

$$ \frac{(n-1)!}{n^{n-1}}=\prod_{m \mid n}\frac{\operatorname{GF}(m)}{m^{\phi(m)}}. $$

Введем обозначения

$$ A(n)=\frac{(n-1)!}{n^{n-1}},\qquad H(n)=\frac{\operatorname{GF}(n)}{n^{\phi(n)}}. $$

Тогда \(A(n)=\prod_{d\mid n}H(d)\), а мультипликативная форма инверсии Мёбиуса дает

$$ H(n)=\prod_{d \mid n} A(d)^{\mu(n/d)} =\prod_{d \mid n}\left(\frac{(d-1)!}{d^{d-1}}\right)^{\mu(n/d)}. $$

Если затем вернуть множитель \(n^{\phi(n)}\) и сократить степени \(n\), получится удобная для реализации формула:

$$ \operatorname{GF}(n)=\prod_{d \mid n}\left(\left(\frac{n}{d}-1\right)! \, d^{\,n/d-1}\right)^{\mu(d)}. $$

Она верна и для \(n=1\), потому что единственный делитель \(d=1\) дает фактор \(0!\cdot 1^0=1\).

Шаг 3: Перемножить по всем \(n\le N\) и поменять порядок произведений

Теперь перемножим эту формулу по всем \(1\le n\le N\):

$$ P(N)=\prod_{n=1}^{N}\prod_{d \mid n}\left(\left(\frac{n}{d}-1\right)! \, d^{\,n/d-1}\right)^{\mu(d)}. $$

Положим \(n=d\,m\). Тогда общее произведение примет вид

$$ P(N)=\prod_{d=1}^{N}\prod_{m=1}^{\lfloor N/d\rfloor}\left((m-1)! \, d^{\,m-1}\right)^{\mu(d)}. $$

Для фиксированного \(d\) обозначим

$$q=\left\lfloor \frac{N}{d}\right\rfloor.$$

Тогда показатель степени при \(d\) равен

$$ \sum_{m=1}^{q}(m-1)=\frac{q(q-1)}{2}=\binom{q}{2}, $$

а факториальная часть равна

$$ \prod_{m=1}^{q}(m-1)! = \prod_{t=1}^{q-1} t!. $$

Шаг 4: Ввести суперфакториал

Определим суперфакториал:

$$ \operatorname{SF}(k)=\prod_{t=1}^{k} t!,\qquad \operatorname{SF}(0)=1. $$

Тогда вся формула сжимается до

$$ \boxed{ P(N)\equiv \prod_{d=1}^{N}\left(d^{\binom{\lfloor N/d\rfloor}{2}}\operatorname{SF}\!\left(\lfloor N/d\rfloor-1\right)\right)^{\mu(d)} \pmod{10^9+7}. } $$

Именно это произведение и вычисляет программа.

Шаг 5: Объединить одинаковые значения \(\left\lfloor N/d\right\rfloor\)

Величина

$$q=\left\lfloor \frac{N}{d}\right\rfloor$$

постоянна на интервале

$$ d\in [l,r],\qquad r=\left\lfloor\frac{N}{q}\right\rfloor. $$

Поэтому выгодно обрабатывать один максимальный интервал целиком. На таком блоке одинаковы и \(\binom{q}{2}\), и \(\operatorname{SF}(q-1)\). Нужны только три агрегата:

$$ P_{+}=\prod_{\substack{d=l\\ \mu(d)=1}}^{r} d,\qquad P_{-}=\prod_{\substack{d=l\\ \mu(d)=-1}}^{r} d,\qquad S_{\mu}=\sum_{d=l}^{r}\mu(d). $$

Вклад интервала равен

$$ P_{+}^{\binom{q}{2}} \cdot P_{-}^{-\binom{q}{2}} \cdot \operatorname{SF}(q-1)^{S_{\mu}}. $$

Поскольку модуль \(10^9+7\) прост, а \(N<10^9+7\), все встречающиеся основания обратимы по модулю. Значит, отрицательные показатели от функции Мёбиуса можно обрабатывать через модульные обратные, либо эквивалентно сводить показатели по модулю \(10^9+6\).

Разобранный пример: \(N=10\)

Для \(N=10\) блоки по значению \(\left\lfloor 10/d\right\rfloor\) таковы:

$$ \left\lfloor\frac{10}{d}\right\rfloor= \begin{cases} 10,& d=1,\\ 5,& d=2,\\ 3,& d=3,\\ 2,& d=4,5,\\ 1,& d=6,7,8,9,10. \end{cases} $$

Значения функции Мёбиуса:

$$ \mu(1)=1,\ \mu(2)=-1,\ \mu(3)=-1,\ \mu(4)=0,\ \mu(5)=-1,\ \mu(6)=1,\ \mu(7)=-1,\ \mu(8)=0,\ \mu(9)=0,\ \mu(10)=1. $$

Все слагаемые с \(q=1\) дают \(d^0\operatorname{SF}(0)=1\), а члены с \(\mu(d)=0\) просто исчезают. Поэтому

$$ P(10)=\frac{\operatorname{SF}(9)}{2^{10}\operatorname{SF}(4)\cdot 3^3\operatorname{SF}(2)\cdot 5}. $$

Численно

$$ \operatorname{SF}(9)=1!\,2!\,3!\cdots 9!=1834933472251084800000, $$

а знаменатель равен

$$ 2^{10}\operatorname{SF}(4)\cdot 3^3\operatorname{SF}(2)\cdot 5 = 79626240. $$

Следовательно,

$$ P(10)=23044331520000, $$

и по модулю \(10^9+7\) получаем

$$ P(10)\equiv 331358692. $$

Это совпадает с контрольным значением, используемым в реализации.

Как работает код

Реализации на C++, Python и Java устроены одинаково. Сначала они вычисляют значения \(\mu(1),\dots,\mu(N)\) с помощью линейного решета, то есть помечают каждое число как свободное от квадратов со знаком \(\pm1\) либо отбрасывают его, если \(\mu(d)=0\).

Затем перечисляются максимальные интервалы, на которых \(q=\lfloor N/d\rfloor\) постоянно, и собираются все ключи вида \(q-1\). После этого факториалы и суперфакториалы строятся последовательно от \(1\) вверх, причем значения суперфакториала сохраняются только в тех точках, которые действительно понадобятся позже.

Для каждого интервала программа перемножает все \(d\) с положительным значением \(\mu(d)\), отдельно перемножает все \(d\) с отрицательным значением и одновременно суммирует \(\mu(d)\) на интервале. На основе этих трех агрегатов выполняется одно обновление, эквивалентное полному вкладу

$$ P_{+}^{\binom{q}{2}} \cdot P_{-}^{-\binom{q}{2}} \cdot \operatorname{SF}(q-1)^{S_{\mu}} $$

по модулю \(10^9+7\). Все модульные степени считаются быстрым двоичным возведением в степень, а показатели сокращаются по модулю \(10^9+6\) по малой теореме Ферма.

Анализ сложности

Линейное решето для \(\mu\) требует \(O(N)\) времени и \(O(N)\) памяти. Построение факториалов и суперфакториалов до \(N-1\) дает еще один проход за \(O(N)\). Обработка блоков по частному в сумме посещает каждый \(d\) ровно один раз; различных плато только \(O(\sqrt{N})\), а на каждое плато приходится лишь логарифмическое число операций быстрого возведения в степень. Поэтому суммарная сложность равна \(O(N)\), а память также составляет \(O(N)\).

Сноски и ссылки

  1. Страница задачи: https://projecteuler.net/problem=754
  2. Функция Мёбиуса: Wikipedia — Möbius function
  3. Формула инверсии Мёбиуса: Wikipedia — Möbius inversion formula
  4. Функция Эйлера \(\phi\): Wikipedia — Euler's totient function
  5. Суперфакториал: Wikipedia — Superfactorial
  6. Линейное решето: cp-algorithms — Linear Sieve

ملخص المسألة

لكل عدد صحيح موجب \(n\)، نعرّف عاملي غاوس بالصيغة

$$\operatorname{GF}(n)=\prod_{\substack{1\le k\lt n\\ \gcd(k,n)=1}} k,$$

مع الاتفاق على أن \(\operatorname{GF}(1)=1\). المطلوب هو حساب

$$P(N)=\prod_{n=1}^{N}\operatorname{GF}(n)\pmod{10^9+7}$$

عندما يكون \(N=10^8\). التنفيذ المباشر غير عملي لأنه يمر على عدد هائل من الأزواج المتباينة نسبيًا، لذلك يعيد الحل كتابة الناتج كله بصيغة تعتمد على القواسم، ثم يستفيد من هذه الصيغة للوصول إلى خوارزمية خطية.

المنهج الرياضي

الفكرة الأساسية هي تحويل كل \(\operatorname{GF}(n)\) إلى حاصل ضرب على القواسم باستخدام قلب موبيوس، ثم تجميع الحدود بحسب المجالات التي تبقى فيها القيمة \(\left\lfloor N/d\right\rfloor\) ثابتة.

الخطوة 1: تجميع حدود \((n-1)!\) بحسب القاسم المشترك الأكبر مع \(n\)

ثبّت \(n\). أي عدد \(k\) يحقق \(1\le k\lt n\) له قيمة

$$d=\gcd(k,n),$$

ومن ثم يمكن كتابته على الصورة

$$k=d\,m,\qquad \gcd\!\left(m,\frac{n}{d}\right)=1,\qquad 1\le m\lt\frac{n}{d}.$$

إذا ثبّتنا قاسمًا \(d\mid n\)، فإن حاصل ضرب كل الأعداد في هذه الفئة يساوي

$$d^{\phi(n/d)}\operatorname{GF}\!\left(\frac{n}{d}\right).$$

وبضرب كل الفئات الممكنة نحصل على

$$ (n-1)! = \prod_{\substack{d \mid n\\ d\lt n}} d^{\phi(n/d)}\operatorname{GF}\!\left(\frac{n}{d}\right). $$

وبالتبديل \(m=n/d\) تصبح الصيغة

$$ (n-1)! = \prod_{\substack{m \mid n\\ m\gt 1}} \left(\frac{n}{m}\right)^{\phi(m)}\operatorname{GF}(m). $$

الخطوة 2: التطبيع وتطبيق قلب موبيوس

بما أن

$$\sum_{\substack{m \mid n\\ m\gt 1}}\phi(m)=n-1,$$

فإن قسمة الطرفين على \(n^{n-1}\) تعطي

$$ \frac{(n-1)!}{n^{n-1}}=\prod_{m \mid n}\frac{\operatorname{GF}(m)}{m^{\phi(m)}}. $$

نعرّف الآن

$$ A(n)=\frac{(n-1)!}{n^{n-1}},\qquad H(n)=\frac{\operatorname{GF}(n)}{n^{\phi(n)}}. $$

عندئذٍ يكون \(A(n)=\prod_{d\mid n}H(d)\)، وبصيغة القلب الضربي لموبيوس نحصل على

$$ H(n)=\prod_{d \mid n} A(d)^{\mu(n/d)} =\prod_{d \mid n}\left(\frac{(d-1)!}{d^{d-1}}\right)^{\mu(n/d)}. $$

بعد إعادة ضرب العامل \(n^{\phi(n)}\) وتبسيط قوى \(n\)، نصل إلى الصيغة الأنسب للتنفيذ:

$$ \operatorname{GF}(n)=\prod_{d \mid n}\left(\left(\frac{n}{d}-1\right)! \, d^{\,n/d-1}\right)^{\mu(d)}. $$

وهذه الصيغة صحيحة أيضًا عندما \(n=1\)، لأن القاسم الوحيد هو \(d=1\) والعامل الناتج يساوي \(0!\cdot 1^0=1\).

الخطوة 3: الضرب على جميع \(n\le N\) ثم تبديل ترتيب الحواصل

إذا ضربنا الصيغة السابقة لكل \(1\le n\le N\)، نحصل على

$$ P(N)=\prod_{n=1}^{N}\prod_{d \mid n}\left(\left(\frac{n}{d}-1\right)! \, d^{\,n/d-1}\right)^{\mu(d)}. $$

نكتب \(n=d\,m\)، وعندها يصبح الناتج الكلي

$$ P(N)=\prod_{d=1}^{N}\prod_{m=1}^{\lfloor N/d\rfloor}\left((m-1)! \, d^{\,m-1}\right)^{\mu(d)}. $$

ولقيمة ثابتة من \(d\)، إذا وضعنا

$$q=\left\lfloor \frac{N}{d}\right\rfloor,$$

فإن أس \(d\) الكلي يساوي

$$ \sum_{m=1}^{q}(m-1)=\frac{q(q-1)}{2}=\binom{q}{2}, $$

والجزء الخاص بالعوامل يساوي

$$ \prod_{m=1}^{q}(m-1)! = \prod_{t=1}^{q-1} t!. $$

الخطوة 4: إدخال دالة السوبرفاكتوريال

نعرّف السوبرفاكتوريال بالصيغة

$$ \operatorname{SF}(k)=\prod_{t=1}^{k} t!,\qquad \operatorname{SF}(0)=1. $$

وبذلك يختصر الناتج الكامل إلى

$$ \boxed{ P(N)\equiv \prod_{d=1}^{N}\left(d^{\binom{\lfloor N/d\rfloor}{2}}\operatorname{SF}\!\left(\lfloor N/d\rfloor-1\right)\right)^{\mu(d)} \pmod{10^9+7}. } $$

وهذه هي الصيغة الرياضية التي تعتمد عليها التطبيقات مباشرة.

الخطوة 5: تجميع القيم المتساوية لـ \(\left\lfloor N/d\right\rfloor\)

القيمة

$$q=\left\lfloor \frac{N}{d}\right\rfloor$$

تبقى ثابتة على مجال من الشكل

$$ d\in [l,r],\qquad r=\left\lfloor\frac{N}{q}\right\rfloor. $$

لذلك لا حاجة لإعادة الحساب لكل \(d\) على حدة. في كل مجال أعظمي تكون كل من \(\binom{q}{2}\) و\(\operatorname{SF}(q-1)\) ثابتتين، ونحتاج فقط إلى ثلاث كميات مجمعة:

$$ P_{+}=\prod_{\substack{d=l\\ \mu(d)=1}}^{r} d,\qquad P_{-}=\prod_{\substack{d=l\\ \mu(d)=-1}}^{r} d,\qquad S_{\mu}=\sum_{d=l}^{r}\mu(d). $$

وعندها تصبح مساهمة المجال

$$ P_{+}^{\binom{q}{2}} \cdot P_{-}^{-\binom{q}{2}} \cdot \operatorname{SF}(q-1)^{S_{\mu}}. $$

ولأن \(10^9+7\) عدد أولي، ولأن \(N<10^9+7\)، فإن كل أساس يظهر هنا قابل للعكس بترديد هذا العدد. لذلك يمكن التعامل مع الأسس السالبة الناتجة عن \(\mu(d)=-1\) باستخدام المعكوسات الضربية، أو بشكل مكافئ باختزال الأسس modulo \(10^9+6\).

مثال محلول: \(N=10\)

عندما \(N=10\)، تكون مجالات القيم \(\left\lfloor 10/d\right\rfloor\) كما يلي:

$$ \left\lfloor\frac{10}{d}\right\rfloor= \begin{cases} 10,& d=1,\\ 5,& d=2,\\ 3,& d=3,\\ 2,& d=4,5,\\ 1,& d=6,7,8,9,10. \end{cases} $$

وقيم دالة موبيوس هي

$$ \mu(1)=1,\ \mu(2)=-1,\ \mu(3)=-1,\ \mu(4)=0,\ \mu(5)=-1,\ \mu(6)=1,\ \mu(7)=-1,\ \mu(8)=0,\ \mu(9)=0,\ \mu(10)=1. $$

كل الحدود التي فيها \(q=1\) تعطي \(d^0\operatorname{SF}(0)=1\)، وكل حد له \(\mu(d)=0\) يختفي. لذا

$$ P(10)=\frac{\operatorname{SF}(9)}{2^{10}\operatorname{SF}(4)\cdot 3^3\operatorname{SF}(2)\cdot 5}. $$

عدديًا،

$$ \operatorname{SF}(9)=1!\,2!\,3!\cdots 9!=1834933472251084800000, $$

أما المقام فهو

$$ 2^{10}\operatorname{SF}(4)\cdot 3^3\operatorname{SF}(2)\cdot 5 = 79626240. $$

إذن

$$ P(10)=23044331520000, $$

وبتخفيضه modulo \(10^9+7\) نحصل على

$$ P(10)\equiv 331358692. $$

وهذه القيمة تطابق نقطة التحقق المستخدمة في التنفيذ.

كيف يعمل الكود

تتبع تطبيقات C++ وPython وJava البنية نفسها. في البداية تُحسب قيم \(\mu(1),\dots,\mu(N)\) باستخدام غربال خطي، بحيث يُصنَّف كل عدد إما كعدد خالٍ من مربعات بعلامة \(\pm1\)، أو يُهمَل إذا كانت \(\mu(d)=0\).

بعد ذلك تُستخرج المجالات العظمى التي تكون فيها \(q=\lfloor N/d\rfloor\) ثابتة، وتُجمع القيم المطلوبة من النوع \(q-1\). ثم تُبنى قيم المضروب والسوبرفاكتوريال تدريجيًا ابتداءً من \(1\)، مع تخزين قيم السوبرفاكتوريال فقط عند النقاط التي سيحتاجها التنفيذ لاحقًا.

في كل مجال، تضرب الخوارزمية قيم \(d\) التي لها \(\mu(d)=1\) في حاصل واحد، وتضرب قيم \(d\) التي لها \(\mu(d)=-1\) في حاصل ثانٍ، كما تجمع قيم \(\mu(d)\) على المجال كله. بهذه الكميات الثلاث يمكن تنفيذ التحديث المكافئ للحد الكامل

$$ P_{+}^{\binom{q}{2}} \cdot P_{-}^{-\binom{q}{2}} \cdot \operatorname{SF}(q-1)^{S_{\mu}} $$

بترديد \(10^9+7\). وتُستخدم الأسس الثنائية السريعة لكل عمليات القوى، مع اختزال الأسس modulo \(10^9+6\) اعتمادًا على مبرهنة فيرما الصغرى.

تحليل التعقيد

الغربال الخطي لحساب \(\mu\) يستهلك زمنًا \(O(N)\) وذاكرة \(O(N)\). وبناء قيم المضروب والسوبرفاكتوريال حتى \(N-1\) هو مرور آخر بتكلفة \(O(N)\). أما مرحلة تجميع المجالات حسب خارج القسمة فتمر على كل \(d\) مرة واحدة إجمالًا، ولا يوجد سوى \(O(\sqrt{N})\) مجالات مميزة، مع عدد لوغاريتمي من عمليات الأس المعياري في كل مجال. لذلك يكون التعقيد الزمني الكلي \(O(N)\)، واستهلاك الذاكرة \(O(N)\).

حواشٍ ومراجع

  1. صفحة المسألة: https://projecteuler.net/problem=754
  2. دالة موبيوس: Wikipedia — Möbius function
  3. صيغة قلب موبيوس: Wikipedia — Möbius inversion formula
  4. دالة أويلر \(\phi\): Wikipedia — Euler's totient function
  5. السوبرفاكتوريال: Wikipedia — Superfactorial
  6. الغربال الخطي: cp-algorithms — Linear Sieve