Detailed mathematical approach

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