Detailed mathematical approach

Problem Summary

For a prime modulus \(M\), an \((N,M)\)-functional inverse of \(x^2\) is a lowest-degree monic polynomial \(Q(x)\) whose coefficients are represented in \(\{0,\ldots,M-1\}\) and which satisfies

\[Q(n^2)\equiv n\pmod M\qquad(0\le n\le N).\]

Here \(N=10^7\), \(M=p=10^9+7\), and only the coefficient of \(x^{10}\) is required. Constructing the degree-\(N\) interpolation polynomial, or even storing all its coefficients, is unnecessary: the calculation can be truncated rigorously at degree \(9\).

Conventions and the meaning of the inverse

The notation \([x^r]F(x)\) means the coefficient of \(x^r\) in \(F\). All equalities below are in \(\mathbb F_p\), unless explicitly presented as rational identities. A fraction \(a/b\) means \(a\,b^{-1}\pmod p\), not integer division or floating-point division. Since \(p\) is prime, every nonzero denominator has a unique inverse. Choosing its representative between \(0\) and \(p-1\) gives the coefficient format required by the problem.

The task does not ask for an inverse of squaring on the whole field: \(n\) and \(-n\) have the same square. It prescribes one value only at each of the selected square nodes. The bound \(2N<p\) ensures that two selected positive indices never supply opposite square roots of the same node. This is why interpolation, rather than a power-series expansion of the real square root, is the appropriate starting point.

Mathematical Approach

Interpolation at the square nodes

Work in the field \(\mathbb F_p\). The nodes \(0^2,1^2,\ldots,N^2\) are distinct: if \(i^2\equiv j^2\pmod p\), then \(p\mid(i-j)(i+j)\), but \(0\le i,j\le N\) and \(i+j\le2N<p\). Hence there is a unique polynomial \(P_N\) of degree at most \(N\) such that \(P_N(n^2)=n\) for every required node.

Write that polynomial in Newton form:

\[P_N(x)=\sum_{k=1}^{N}d_k\prod_{j=0}^{k-1}(x-j^2).\]

The term \(k=0\) vanishes because the prescribed value at \(0\) is \(0\).

Here is an explicit existence and uniqueness argument. For a fixed \(k\), define the Lagrange basis below. At \(x=j^2\), its \(i\)-th member is \(1\) when \(i=j\) and \(0\) otherwise. Thus their weighted sum has the prescribed values. If two degree-at-most-\(k\) polynomials had those values, their difference would have \(k+1\) distinct roots; a nonzero polynomial of degree at most \(k\) cannot do so.

\[L_{i,k}(x)=\prod_{\substack{0\le j\le k\\j\ne i}}\frac{x-j^2}{i^2-j^2},\qquad P_k(x)=\sum_{i=0}^{k}iL_{i,k}(x).\]

Newton's form adds the nodes one at a time. The difference \(P_k-P_{k-1}\) vanishes at the first \(k\) nodes, so it is a constant multiple of their monic product \(B_k\). Evaluating at the new node determines that constant. Since \(P_{k-1}\) has degree at most \(k-1\), \(d_k\) is also the coefficient of \(x^k\) in the Lagrange expression for \(P_k\). Reading off that leading coefficient gives the divided-difference sum used next.

\[\begin{aligned}B_k(x)&=\prod_{j=0}^{k-1}(x-j^2),\\P_k(x)&=P_{k-1}(x)+d_kB_k(x),\\d_k&=\frac{k-P_{k-1}(k^2)}{B_k(k^2)}=[x^k]P_k(x).\end{aligned}\]

Closed form for the Newton coefficients

The divided difference on \(0^2,\ldots,k^2\) is

\[d_k=\sum_{i=1}^{k}\frac{i} {\displaystyle\prod_{\substack{0\le j\le k\\j\ne i}}(i^2-j^2)}.\]

To evaluate the denominator, keep the two linear products separate. In the first product there are exactly \(k-i\) negative factors. In the second, the consecutive product from \(i\) to \(i+k\) must have its missing factor \(2i\) removed. For \(1\le i\le k\), this gives

\[\begin{aligned}\prod_{\substack{0\le j\le k\\j\ne i}}(i-j)&=i!\,(-1)^{k-i}(k-i)!,\\\prod_{\substack{0\le j\le k\\j\ne i}}(i+j)&=\frac{(i+k)!}{(i-1)!\,2i}.\end{aligned}\]

Multiplying these expressions cancels \(i!/(i-1)!=i\) against \(2i\). All factorials involved are invertible modulo \(p\), since their arguments are at most \(2k\le2N<p\). Thus

\[\prod_{\substack{0\le j\le k\\j\ne i}}(i^2-j^2)=\frac{(-1)^{k-i}}{2}(k-i)!(k+i)!.\]

\[d_k=\sum_{i=1}^{k}\frac{2i(-1)^{k-i}}{(k-i)!(k+i)!}.\]

The remaining sum can be evaluated without leaving an unexplained cancellation. Multiply by \((2k)!\), put \(r=k-i\), and name the result \(S_k\). The factorial denominator becomes a binomial coefficient:

\[S_k=(2k)!d_k=\sum_{r=0}^{k-1}(-1)^r\,2(k-r)\binom{2k}{r}.\]

\[2(k-r)\binom{2k}{r}=2k\left(\binom{2k-1}{r}-\binom{2k-1}{r-1}\right).\]

To handle the alternating sums, use Pascal's identity in each summand. After shifting the second sum by one index, every interior term cancels and only the last term remains. We use the convention that a binomial coefficient with lower index \(-1\) is zero and that an empty sum is zero.

\[\begin{aligned}\sum_{r=0}^{m}(-1)^r\binom{n}{r}&=\sum_{r=0}^{m}(-1)^r\binom{n-1}{r}-\sum_{r=0}^{m-1}(-1)^r\binom{n-1}{r}\\&=(-1)^m\binom{n-1}{m}.\end{aligned}\]

Substitute the binomial difference into \(S_k\). The minus sign in that difference becomes a plus after the index shift, because the sign also changes. Apply the partial-sum identity twice, and then use the ratio of the two adjacent binomial coefficients. The calculation is valid for \(k=1\) as well under the empty-sum convention:

\[\begin{aligned}\frac{S_k}{2k}&=\sum_{r=0}^{k-1}(-1)^r\binom{2k-1}{r}+\sum_{r=0}^{k-2}(-1)^r\binom{2k-1}{r}\\&=(-1)^{k-1}\left(\binom{2k-2}{k-1}-\binom{2k-2}{k-2}\right)\\&=\frac{(-1)^{k-1}}{k}\binom{2k-2}{k-1}.\end{aligned}\]

\[d_k=\frac{2(-1)^{k-1}}{(2k)!}\binom{2k-2}{k-1}=\boxed{\frac{(-1)^{k-1}}{(2k-1)k!(k-1)!}}.\]

Normalize the products and retain only degree \(9\)

Define

\[A_k(x)=\prod_{j=1}^{k}\left(1-\frac{x}{j^2}\right), \qquad A_0(x)=1.\]

Since

\[\prod_{j=0}^{k-1}(x-j^2) =x(-1)^{k-1}((k-1)!)^2A_{k-1}(x),\]

the signs and factorials cancel against \(d_k\), leaving the particularly simple identity

\[\boxed{\frac{P_N(x)}x= \sum_{k=1}^{N}\frac{A_{k-1}(x)}{k(2k-1)}}.\]

Consequently \([x^{10}]P_N=[x^9](P_N/x)\). Put \(a_d^{(k)}=[x^d]A_k\) and let \(s_d^{(k)}\) denote the accumulated coefficient after terms \(1,\ldots,k\). For \(1\le d\le9\), each step is

\[s_d^{(k)}=s_d^{(k-1)} +\frac{a_d^{(k-1)}}{k(2k-1)},\]

\[a_d^{(k)}=a_d^{(k-1)} -\frac{a_{d-1}^{(k-1)}}{k^2}.\]

The \(a_d\) array is updated from degree \(9\) down to \(1\), so every right-hand side still belongs to \(A_{k-1}\). Dropping terms above degree \(9\) is exact, not an approximation: later multiplications are by polynomials with nonnegative degree and can never move a discarded coefficient down to degree \(9\).

More explicitly, the scalar multiplying \(A_{k-1}\) is \(d_k(-1)^{k-1}((k-1)!)^2=1/(k(2k-1))\). Initialize the coefficient state and the sum as follows; the constant coefficient of \(A_k\) remains \(1\) throughout, so it never needs updating.

\[a_0^{(0)}=1,\qquad a_d^{(0)}=0\ (1\le d\le9),\qquad s_d^{(0)}=0\ (0\le d\le9).\]

There are two separate ordering requirements. First accumulate the contribution from \(A_{k-1}\), then multiply by the new factor to obtain \(A_k\). Inside that multiplication, update degrees downwards. An upward update would use the already changed \(a_{d-1}\). Starting from \(1\), for example, it would incorrectly create an \(x^2/k^4\) term after multiplying by just \(1-x/k^2\), even though that product has degree one.

The formal reason truncation is safe is that we compute in the quotient ring \(\mathbb F_p[x]/(x^{10})\). Reduction modulo \(x^{10}\) commutes with addition and multiplication. Equivalently, the coefficient of degree \(d\le9\) in a product uses only input coefficients of degrees at most \(d\). No discarded high-degree term can affect any stored coefficient:

\[[x^d](F(x)G(x))=\sum_{r=0}^{d}([x^r]F(x))([x^{d-r}]G(x)),\qquad 0\le d\le9.\]

A serial run therefore stores ten product coefficients but only the scalar \(s_9\). The other sum coefficients never feed into \(s_9\). All ten product coefficients are still essential: repeated updates can move information from degree zero through the intermediate degrees to degree nine. Also, \(a_9^{(k-1)}=0\) for \(k<10\), which explains why the first nine summands contribute nothing to the interpolation part of the requested coefficient.

A small worked interpolation example

Take \(N=3\) and keep rational fractions temporarily to make the algebra readable. These identities can then be reduced modulo any prime \(p>6\). The first normalized products are

\[A_0=1,\qquad A_1=1-x,\qquad A_2=(1-x)(1-x/4)=1-\frac54x+\frac14x^2.\]

The weights for \(k=1,2,3\) are \(1,1/6,1/15\). Each weight uses the product before its own factor has been inserted. Substituting and collecting coefficients yields

\[\begin{aligned}\frac{P_3(x)}x&=1+\frac{1-x}{6}+\frac{1-5x/4+x^2/4}{15}\\&=\frac{37}{30}-\frac{x}{4}+\frac{x^2}{60},\\P_3(x)&=\frac{37}{30}x-\frac14x^2+\frac1{60}x^3.\end{aligned}\]

\[P_3(0)=0,\qquad P_3(1)=1,\qquad P_3(4)=2,\qquad P_3(9)=3.\]

This checks both the prescribed values and the coefficient formula \(d_3=1/60\). It also exposes the distinction between interpolation and the monic condition: for the target modulus, this small interpolant has leading coefficient \(1/60\ne1\), so interpolation alone is not the complete construction. The monic correction is derived below.

Batch inversion removes millions of exponentiations

Let \(u_k=k(2k-1)\). All \(u_k\) are nonzero in \(\mathbb F_p\), because \(1\le k\le N<p\) and \(2k-1<p\). In a batch \(u_1,\ldots,u_m\), form prefix products \(U_0=1\), \(U_i=U_{i-1}u_i\). Compute only \(U_m^{-1}=U_m^{p-2}\) and walk backwards:

\[u_i^{-1}=U_{i-1}U_i^{-1}, \qquad U_{i-1}^{-1}=U_i^{-1}u_i.\]

Thus one Fermat exponentiation supplies every reciprocal in a batch. Moreover \(k^{-1}=(2k-1)u_k^{-1}\), so squaring it gives \(k^{-2}\) for the product update. The implementations use batches of \(32768\) entries.

The subscripts in the batch formula are local to that batch; its first global index need not be \(1\). During the backward pass, maintain the invariant that the running inverse equals \(U_i^{-1}\). Multiplying by the still-unmodified prefix \(U_{i-1}\) yields \(u_i^{-1}\), and multiplying the running inverse by \(u_i\) prepares the next iteration. The prefix slot can then be overwritten by its recovered inverse, because subsequent iterations only need earlier slots.

As a concrete check of this general trick, invert \(2,3,5\) modulo \(11\). The prefix products are \(1,2,6,8\), and \(8^{-1}=7\). The backward pass first gives \(5^{-1}=6\cdot7=9\), then \(3^{-1}=2\cdot2=4\), then \(2^{-1}=1\cdot6=6\), all modulo \(11\). Only one exponentiation was required. In the actual algorithm the inputs are the nonzero values \(k(2k-1)\), and every multiplication is reduced modulo \(p\).

Making the interpolant monic

Because each Newton basis polynomial is monic, the leading coefficient of \(P_N\) is

\[c_N=d_N= \frac{(-1)^{N-1}N}{(2N-1)(N!)^2}.\]

The polynomial that vanishes at every interpolation node is

\[V_N(x)=\prod_{j=0}^{N}(x-j^2) =(-1)^N(N!)^2xA_N(x).\]

The formula for \(d_N\) is nonzero: its numerator is a sign and none of its denominator factors vanishes. Therefore \(P_N\) has degree exactly \(N\), not merely at most \(N\). To describe every possible solution, let \(F\) be any polynomial with the prescribed values. Each linear factor \(x-j^2\) divides \(F-P_N\); because the roots are distinct, their product divides it. Hence

\[F(x)=P_N(x)+V_N(x)R(x),\qquad R(x)\in\mathbb F_p[x].\]

If \(R\ne0\), then \(\deg(V_NR)=N+1+\deg R>\deg P_N\), so the highest term cannot cancel. Thus the only solution of degree at most \(N\) is \(P_N\). If it is not monic, a monic solution of the next degree requires \(R\) to be a constant, and its leading coefficient forces that constant to be exactly \(1\). This proves both minimality and uniqueness. Simply dividing \(P_N\) by \(c_N\) would instead change the required value at \(1\) from \(1\) to \(c_N^{-1}\), so scaling is not an admissible repair.

If \(c_N=1\), the interpolant itself is already the unique lowest-degree monic solution. Otherwise no polynomial of degree at most \(N\) can be another solution, and the lowest possible monic degree is \(N+1\); adding the monic polynomial \(V_N\) gives

\[Q(x)=P_N(x)+V_N(x).\]

For the target instance the code verifies \(c_N\ne1\). The desired coefficient is therefore obtained without constructing either large polynomial:

\[\boxed{[x^{10}]Q=s_9^{(N)} +(-1)^N(N!)^2a_9^{(N)}}.\]

Note on the published small example. Read literally, the “lowest-degree” wording has a special edge case at \((N,M)=(2,7)\): there \(c_2=1\), so \(P_2(x)=x^2\) is already monic and satisfies the three values. The displayed cubic is \(P_2+V_2\), which satisfies them as well. This editorial edge case does not affect the requested instance, for which \(c_N\ne1\).

For a separate example with no monic edge case, take \((N,p)=(2,11)\). The rational interpolant \((7x-x^2)/6\) reduces to \(3x+9x^2\). Adding the vanishing polynomial gives the following cubic; its values at \(0,1,4\) are respectively \(0,1,2\) modulo \(11\), and its leading coefficient is \(1\). This illustrates why the correction changes coefficients while preserving every required value.

\[\begin{aligned}V_2(x)&=x(x-1)(x-4)\equiv x^3+6x^2+4x\pmod{11},\\Q(x)&=P_2(x)+V_2(x)\equiv x^3+4x^2+7x\pmod{11}.\end{aligned}\]

Combining independently processed blocks

The C++ program can process consecutive index intervals in parallel. For a block \(B\), write \(g_k(x)=1-x/k^2\) and store, through degree \(9\),

\[G_B=\prod_{k\in B}g_k,\qquad H_B=\sum_{k\in B}\frac1{k(2k-1)} \prod_{\substack{j\in B\\j<k}}g_j.\]

If block \(L\) immediately precedes block \(R\), their summaries combine as

\[\boxed{G_{LR}=G_LG_R,\qquad H_{LR}=H_L+G_LH_R}.\]

This operation is associative, but order matters; the implementation merges blocks in increasing-index order. Their factorial products are multiplied at the same time. Truncated polynomial multiplication is sufficient because coefficients above degree \(9\) cannot contribute below that degree.

The block rule can also be read as an update of the incoming global state: \((S,A)\mapsto(S+AH_B,AG_B)\). Contributions from the right block need every factor in the left block before them, which accounts for \(G_LH_R\). For three consecutive blocks, either parenthesization produces exactly the same expression below, proving associativity. The empty block has summary \((G,H)=(1,0)\).

\[G_{LMR}=G_LG_MG_R,\qquad H_{LMR}=H_L+G_LH_M+G_LG_MH_R.\]

This also explains why a worker must retain all ten coefficients of \(H_B\), even though the final answer needs only degree nine. An earlier product mixes the degrees according to the convolution below. The merging code may keep only the ninth coefficient of the global sum, but not only the ninth coefficient of each worker's sum. Inversion batches and parallel blocks serve different purposes: a worker can process many inversion batches within its own consecutive interval.

\[[x^9](AH_B)=\sum_{r=0}^{9}([x^r]A)([x^{9-r}]H_B).\]

Correctness Argument

Lemma 1. The square nodes are distinct in \(\mathbb F_p\), so \(P_N\) exists and is the only interpolant of degree at most \(N\).

Lemma 2. Evaluating the divided-difference denominator and the alternating binomial sum gives the stated \(d_k\). Substitution into Newton's formula gives \(P_N/x=\sum A_{k-1}/(k(2k-1))\).

Lemma 3. The descending recurrence maintains \(a_d^{(k)}=[x^d]A_k\), while the accumulation maintains \(s_d^{(k)}=[x^d]\sum_{i=1}^{k}A_{i-1}/(i(2i-1))\). Degree-\(9\) truncation preserves these invariants exactly.

Lemma 4. Batch inversion returns the true inverse of each nonzero \(u_k\): \(U_{i-1}U_i^{-1}=1/u_i\). Hence it changes only how the same recurrence is evaluated.

Lemma 5. Every polynomial satisfying the data differs from \(P_N\) by a multiple of \(V_N\). The leading-coefficient test therefore selects \(P_N\) exactly when it is already monic and otherwise selects the monic polynomial \(P_N+V_N\) of the smallest possible degree.

Lemma 6. The ordered block identity inserts \(G_L\) before every contribution originating in \(R\); induction over merged blocks therefore reproduces the serial sum and product.

Combining the lemmas, the final expression is exactly the coefficient of \(x^{10}\) in the required lowest-degree monic functional inverse.

How the Code Works

In C++, Block.product[d] represents the truncated product coefficient and Block.sum[d] represents the corresponding weighted sum. The serial specialization of calculate_block keeps only sum[9]; parallel workers retain the full sum. Python and Java store the same product in coefficients and the ninth sum coefficient in interpolation_coefficient and interpolationCoefficient, respectively. The mathematical state is the same in all three implementations.

At the start of serial iteration \(k\), the array contains \(A_{k-1}\), the accumulator contains \(s_9^{(k-1)}\), and the factorial variable contains \((k-1)!\). Retrieve the batch weight, add its contribution to the accumulator, update the array downwards, and multiply the factorial variable by \(k\). These actions re-establish the same invariant for iteration \(k+1\). After the final iteration the state is exactly \((A_N,s_9^{(N)},N!)\), which is all the final formula needs.

The C++ default uses at most \(16\) worker threads, while its optional --threads argument can select another supported count. Workers process consecutive intervals and merge in their original order. Python and Java evaluate the serial recurrence. Binary modular exponentiation is provided by C++ power, Python's three-argument pow, and Java modPow. For this modulus, multiplying two reduced residues gives less than \((10^9+7)^2<2^{63}\), so C++ 64-bit arithmetic and Java long suffice; a 32-bit product would not suffice.

For validation, each implementation also includes a direct Lagrange interpolator that constructs the whole polynomial for \(1\le N\le40\). It uses a different construction from the optimized Newton recurrence. Comparisons with that reference include the already-monic case at \(N=1\), and check the degree-ten coefficient even when the small polynomial has no such term. The C++ self-tests additionally compare serial and parallel results around inversion-batch boundaries; Python and Java also test a deliberately small batch size.

Complexity Analysis

Let \(D=9\), let \(B=32768\) be the inversion-batch size, and let \(T\) be the number of workers, with \(T=1\) in a serial run. Each index requires \(D\) product-coefficient updates and a constant amount of other arithmetic; a full worker sum adds another \(D+1\) updates. Prefix and backward passes cost linear work in the batch size. The total number of batch inversions is at most \(\lceil N/B\rceil+T-1\), and testing the leading coefficient uses one additional inversion.

\[W(N,D,B,T,p)=O\!\left(ND+\left(\frac NB+T\right)\log p+TD^2\right).\]

The \(TD^2\) term accounts for the short polynomial convolutions at merges. With fixed \(D,B,T\) and this fixed modulus, total work is \(O(N)\). For the target \(N\), a serial run performs \(90{,}000{,}000\) product-coefficient updates and only \(306\) batch inversions, plus the final leading-coefficient inversion. Parallelism divides the main loop among workers, but does not change its total order of work; actual elapsed time also depends on scheduling and hardware.

Auxiliary storage is \(O(B+D)\) in the serial implementation and \(O(T(B+D))\) in the parallel implementation, excluding runtime thread stacks. There is no array indexed over all \(N\) nodes and no degree-\(N\) coefficient array. The same derivation can compute another positive degree \(r\) by retaining degrees \(0\) through \(r-1\); degree ten is what makes the present state only ten coefficients long.

Footnotes and References

1. Problem page: Project Euler 1008 - Functional Inverse
2. Polynomial interpolation: Wikipedia - Polynomial interpolation
3. Newton form and divided differences: Wikipedia - Newton polynomial
4. Finite fields: Wikipedia - Finite field
5. Modular multiplicative inverses: Wikipedia - Modular multiplicative inverse
6. Fermat's little theorem: Wikipedia - Fermat's little theorem