Detailed mathematical approach

Problem Summary

Let \(M_q = 2^q - 1\), where \(q\) is the prime from the problem. The goal is to determine the canonical square root of \(q\) modulo the Mersenne number \(M_q\), and then report that enormous integer modulo \(10^9+7\).

A direct modular-square-root computation is hopeless because \(M_q\) itself has \(q\) binary digits. For the target value \(q = 74{,}207{,}281\), even writing down \(M_q\) is already infeasible. The implementations therefore avoid generic number-theoretic machinery and instead use a character sum whose binary digits are determined by quadratic residues modulo \(q\).

Mathematical Approach

The central idea is that powers of 2 modulo a Mersenne number behave like roots of unity, so a classical quadratic Gauss sum turns into an explicit binary expansion of the desired square root.

A \(q\)-th Root of Unity Hidden in the Mersenne Modulus

Work modulo

\[ M_q = 2^q - 1. \]

Inside this ring we have \(2^q \equiv 1 \pmod{M_q}\), so exponents of 2 can always be reduced modulo \(q\). In other words, the element 2 plays the role of a nontrivial \(q\)-th root of unity.

Now define the quadratic character on the nonzero residue classes modulo \(q\):

\[ \chi(a)= \begin{cases} 1, & a \text{ is a nonzero quadratic residue mod } q, \\ -1, & a \text{ is a quadratic non-residue mod } q. \end{cases} \]

The sum used by the implementations is

\[ G=\sum_{a=1}^{q-1}\chi(a)\,2^a. \]

This is not an arbitrary weighted sum. It is the exact object that encodes the square root.

Why the Gauss Sum Squares to \(q\)

Expand the square of \(G\) modulo \(M_q\):

\[ G^2 = \sum_{a=1}^{q-1}\sum_{b=1}^{q-1}\chi(a)\chi(b)\,2^{a+b}. \]

For each fixed nonzero \(a\), write \(b \equiv at \pmod q\). Since \(\chi(ab)=\chi(a)\chi(b)=\chi(t)\) and exponents may be reduced modulo \(q\), this becomes

\[ G^2 \equiv \sum_{a=1}^{q-1}\sum_{t=1}^{q-1}\chi(t)\,2^{a(1+t)} \pmod{M_q}. \]

Now split into two cases.

If \(t \not\equiv -1 \pmod q\), then multiplication by \(1+t\) permutes the nonzero residue classes modulo \(q\), so

\[ \sum_{a=1}^{q-1}2^{a(1+t)} \equiv \sum_{u=1}^{q-1}2^u = 2^q-2 = M_q-1 \equiv -1 \pmod{M_q}. \]

If \(t \equiv -1 \pmod q\), then every exponent is \(0\) modulo \(q\), so the inner sum is simply \(q-1\).

Therefore

\[ G^2 \equiv \chi(-1)(q-1) - \sum_{t \ne -1}\chi(t) \pmod{M_q}. \]

Because the quadratic character has equally many \(+1\) and \(-1\) values,

\[ \sum_{t=1}^{q-1}\chi(t)=0, \]

so \(\sum_{t \ne -1}\chi(t) = -\chi(-1)\). Substituting gives the standard Gauss-sum identity

\[ G^2 \equiv \chi(-1)\,q = (-1)^{(q-1)/2}q \pmod{M_q}. \]

For the prime used here, \(q \equiv 1 \pmod 4\), hence \(\chi(-1)=1\) and

\[ G^2 \equiv q \pmod{M_q}. \]

So \(G\) is one of the two square roots of \(q\) modulo \(M_q\).

Recovering the Smaller of the Two Roots

Let

\[ R=\sum_{\chi(a)=1}2^a, \qquad N=\sum_{\chi(a)=-1}2^a. \]

Then \(G = R - N\), while

\[ R+N=\sum_{a=1}^{q-1}2^a = M_q - 1. \]

Hence

\[ G = M_q - 1 - 2N, \qquad M_q - G = 1 + 2N. \]

So once the non-residue sum \(N\) is known, the opposite square root is available for free.

The remaining question is: which of \(G\) and \(M_q-G\) is the smaller positive root? The answer depends on \(q \bmod 8\).

The top bit \(2^{q-1}\) always appears in \(G\) with positive sign because

\[ \chi(q-1)=\chi(-1)=1. \]

The next bit \(2^{q-2}\) has sign \(\chi(q-2)=\chi(-2)=\chi(2)\), and the supplementary law gives

\[ \chi(2)=(-1)^{(q^2-1)/8} = \begin{cases} 1, & q \equiv 1 \pmod 8, \\ -1, & q \equiv 5 \pmod 8. \end{cases} \]

If \(q \equiv 1 \pmod 8\), then even in the most pessimistic case for the lower bits,

\[ G \ge 2^{q-1}+2^{q-2}-(2^{q-2}-2)=2^{q-1}+2, \]

so \(G > M_q/2\). Therefore the smaller root is

\[ M_q-G = 1+2N. \]

If \(q \equiv 5 \pmod 8\), then even giving every lower bit the favorable sign,

\[ G \le 2^{q-1}-2^{q-2}+(2^{q-2}-2)=2^{q-1}-2, \]

so \(G < M_q/2\). In that case \(G\) itself is the smaller root.

Worked Examples

For \(q=5\), the nonzero quadratic residues are \(\{1,4\}\) and the non-residues are \(\{2,3\}\). Therefore

\[ G = 2^1 - 2^2 - 2^3 + 2^4 = 2 - 4 - 8 + 16 = 6. \]

Since \(M_5 = 31\), we get

\[ 6^2 = 36 \equiv 5 \pmod{31}. \]

Also \(5 \equiv 5 \pmod 8\), so the smaller root is exactly \(G=6\).

For \(q=17\), the same construction yields \(G=83{,}502\). The opposite root is

\[ M_{17}-G = 2^{17}-1-83{,}502 = 47{,}569. \]

Because \(17 \equiv 1 \pmod 8\), the smaller root is \(47{,}569\), and indeed

\[ 47{,}569^2 \equiv 17 \pmod{131{,}071}. \]

How the Code Works

Marking Quadratic Residues in One Linear Pass

The implementation first builds a boolean table for the nonzero quadratic residues modulo \(q\). It only needs to iterate \(x=1,2,\dots,(q-1)/2\), because \(x^2 \equiv (-x)^2 \pmod q\), so these values already cover every nonzero residue exactly once.

Instead of recomputing each square from scratch, it maintains the invariant

\[ (x+1)^2 \equiv x^2 + (2x+1) \pmod q. \]

Since \(2x+1 \le q\) throughout that loop, each update needs at most one subtraction by \(q\). That is why the residue-generation phase is a tight \(O(q)\) scan with very small constant factors.

Accumulating the Gauss Sum and the Non-Residue Sum

Next, the implementation sweeps through \(n=1,2,\dots,q-1\). At each step it updates the current power \(2^n \bmod 10^9+7\) by a single doubling.

If \(n\) is marked as a quadratic residue, that power is added to the running value of \(G\). If \(n\) is a non-residue, the same power is subtracted from \(G\) and added to the running value of \(N\). So the whole mathematical formula is evaluated in one pass, with only modular additions and subtractions.

Crucially, the code never constructs \(M_q\) and never performs a square-root routine modulo a gigantic number. It only evaluates the explicit binary expansion of the root and reduces the result modulo \(10^9+7\) as it goes.

Applying the \(q \bmod 8\) Branch

After the sweep finishes, the implementation has exactly the two quantities needed for the final decision: \(G\) and \(N\), both already reduced modulo \(10^9+7\).

If \(q \equiv 1 \pmod 8\), the smaller root is \(1+2N\), so that value is returned. Otherwise the smaller root is \(G\), so the implementation returns the signed Gauss sum directly. The C++, Python, and Java versions all follow this same logic.

Complexity Analysis

The algorithm performs two linear scans of length \(q-1\): one to mark quadratic residues, and one to accumulate the weighted powers of 2. Therefore the running time is \(O(q)\).

The memory usage is \(O(q)\) for the residue table. Aside from that table, the algorithm stores only a few modular accumulators and loop counters. This is exactly the right trade-off for the target prime: a byte-level table of length \(q\) is practical, while any attempt to represent \(M_q\) itself would be astronomically out of reach.

Footnotes and References

1. Project Euler problem page: Problem 942
2. Mersenne numbers: Wikipedia - Mersenne number
3. Quadratic residues: Wikipedia - Quadratic residue
4. Legendre symbol: Wikipedia - Legendre symbol
5. Gauss sums: Wikipedia - Gauss sum
6. Supplementary law for \(\left(\frac{2}{q}\right)\): Wikipedia - Quadratic reciprocity