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\).
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.
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}\]
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)!}}.\]
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.
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.
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\).
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}\]
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).\]
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.
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.
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.
Für einen Primmodul \(M\) ist eine \((N,M)\)-funktionale Inverse von \(x^2\) ein normiertes Polynom \(Q(x)\) kleinsten Grades. Seine Koeffizienten werden durch \(\{0,\ldots,M-1\}\) dargestellt, und es erfüllt
\[Q(n^2)\equiv n\pmod M\qquad(0\le n\le N).\]
Hier sind \(N=10^7\) und \(M=p=10^9+7\); gesucht ist nur der Koeffizient von \(x^{10}\). Weder das Interpolationspolynom vom Grad \(N\) noch alle seine Koeffizienten müssen erzeugt werden: Die Rechnung lässt sich beweisbar exakt bei Grad \(9\) abschneiden.
Die Schreibweise \([x^r]F(x)\) bezeichnet den Koeffizienten von \(x^r\) in \(F\). Alle folgenden Gleichheiten gelten in \(\mathbb F_p\), sofern sie nicht ausdrücklich als rationale Identitäten angegeben sind. Ein Bruch \(a/b\) bedeutet \(a\,b^{-1}\pmod p\), keine Ganzzahl- oder Gleitkommadivision. Da \(p\) prim ist, besitzt jeder von null verschiedene Nenner eine eindeutige Inverse. Sein Repräsentant zwischen \(0\) und \(p-1\) liefert das in der Aufgabe geforderte Koeffizientenformat.
Gesucht ist keine Umkehrung des Quadrierens auf dem ganzen Körper: \(n\) und \(-n\) haben dasselbe Quadrat. Vorgeschrieben ist nur je ein Wert an den ausgewählten Quadratstellen. Die Schranke \(2N<p\) verhindert, dass zwei ausgewählte positive Indizes entgegengesetzte Quadratwurzeln derselben Stützstelle liefern. Deshalb ist Interpolation der passende Ausgangspunkt und nicht eine Potenzreihenentwicklung der reellen Quadratwurzel.
Wir rechnen im Körper \(\mathbb F_p\). Die Stützstellen \(0^2,1^2,\ldots,N^2\) sind verschieden: Aus \(i^2\equiv j^2\pmod p\) folgte \(p\mid(i-j)(i+j)\), doch für \(0\le i,j\le N\) gilt \(i+j\le2N<p\). Daher gibt es genau ein Polynom \(P_N\) vom Grad höchstens \(N\) mit \(P_N(n^2)=n\) an allen geforderten Stellen.
Seine Newton-Darstellung lautet
\[P_N(x)=\sum_{k=1}^{N}d_k\prod_{j=0}^{k-1}(x-j^2).\]
Der Term \(k=0\) ist null, weil der vorgeschriebene Wert bei \(0\) ebenfalls null ist.
Existenz und Eindeutigkeit lassen sich ausdrücklich begründen. Für festes \(k\) definiere die folgende Lagrange-Basis. Bei \(x=j^2\) ist ihr \(i\)-tes Element gleich \(1\), falls \(i=j\), und sonst gleich \(0\). Ihre gewichtete Summe nimmt daher die vorgeschriebenen Werte an. Hätten zwei Polynome vom Grad höchstens \(k\) diese Werte, besäße ihre Differenz \(k+1\) verschiedene Nullstellen; für ein von null verschiedenes Polynom dieses Grades ist das unmöglich.
\[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).\]
Die Newton-Form fügt die Stützstellen einzeln hinzu. Die Differenz \(P_k-P_{k-1}\) verschwindet an den ersten \(k\) Stellen und ist somit ein konstantes Vielfaches ihres normierten Produkts \(B_k\). Auswertung an der neuen Stelle bestimmt die Konstante. Weil \(P_{k-1}\) höchstens Grad \(k-1\) hat, ist \(d_k\) zugleich der Koeffizient von \(x^k\) in der Lagrange-Darstellung von \(P_k\). Das Ablesen dieses Leitkoeffizienten ergibt die folgende Summe dividierter Differenzen.
\[\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}\]
Die dividierte Differenz zu \(0^2,\ldots,k^2\) ist
\[d_k=\sum_{i=1}^{k}\frac{i} {\displaystyle\prod_{\substack{0\le j\le k\\j\ne i}}(i^2-j^2)}.\]
Zur Auswertung des Nenners behandeln wir die beiden linearen Produkte getrennt. Das erste enthält genau \(k-i\) negative Faktoren. Im zweiten muss aus dem fortlaufenden Produkt von \(i\) bis \(i+k\) der ausgelassene Faktor \(2i\) entfernt werden. Für \(1\le i\le k\) folgt
\[\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}\]
Beim Multiplizieren kürzt sich \(i!/(i-1)!=i\) gegen \(2i\). Alle auftretenden Fakultäten sind modulo \(p\) invertierbar, denn ihre Argumente sind höchstens \(2k\le2N<p\). Damit gilt
\[\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)!}.\]
Die verbleibende Summe lässt sich ohne unerklärte Kürzung auswerten. Multipliziere mit \((2k)!\), setze \(r=k-i\) und nenne das Ergebnis \(S_k\). Der Fakultätennenner wird dabei zu einem Binomialkoeffizienten:
\[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).\]
Für die alternierenden Summen wende Pascals Identität auf jeden Summanden an. Nach Verschiebung des zweiten Summationsindex um eins kürzen sich alle inneren Terme; nur der letzte bleibt übrig. Ein Binomialkoeffizient mit unterem Index \(-1\) sowie eine leere Summe werden dabei als null definiert.
\[\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}\]
Setze die Binomialdifferenz in \(S_k\) ein. Ihr Minuszeichen wird nach der Indexverschiebung zum Pluszeichen, da zugleich das alternierende Vorzeichen wechselt. Wende die Teilsummenidentität zweimal an und nutze dann das Verhältnis der benachbarten Binomialkoeffizienten. Mit der Konvention für leere Summen gilt die Rechnung auch für \(k=1\):
\[\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)!}}.\]
Definiere
\[A_k(x)=\prod_{j=1}^{k}\left(1-\frac{x}{j^2}\right), \qquad A_0(x)=1.\]
Wegen
\[\prod_{j=0}^{k-1}(x-j^2) =x(-1)^{k-1}((k-1)!)^2A_{k-1}(x),\]
heben sich Vorzeichen und Fakultäten gegen \(d_k\) weg. Es bleibt die einfache Identität
\[\boxed{\frac{P_N(x)}x= \sum_{k=1}^{N}\frac{A_{k-1}(x)}{k(2k-1)}}.\]
Somit ist \([x^{10}]P_N=[x^9](P_N/x)\). Setze \(a_d^{(k)}=[x^d]A_k\), und bezeichne mit \(s_d^{(k)}\) den nach den Termen \(1,\ldots,k\) akkumulierten Koeffizienten. Für \(1\le d\le9\) gilt schrittweise
\[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}.\]
Das Array \(a_d\) wird absteigend von Grad \(9\) bis \(1\) aktualisiert, sodass rechts stets ein Koeffizient von \(A_{k-1}\) steht. Das Weglassen höherer Grade ist exakt und keine Näherung: Spätere Multiplikationen haben nur nichtnegative Grade und können einen verworfenen Koeffizienten nie wieder auf Grad \(9\) absenken.
Ausdrücklich ist der Faktor vor \(A_{k-1}\) gleich \(d_k(-1)^{k-1}((k-1)!)^2=1/(k(2k-1))\). Initialisiere Koeffizientenzustand und Summe wie folgt. Der konstante Koeffizient von \(A_k\) bleibt stets \(1\) und muss daher nie aktualisiert werden.
\[a_0^{(0)}=1,\qquad a_d^{(0)}=0\ (1\le d\le9),\qquad s_d^{(0)}=0\ (0\le d\le9).\]
Zwei Reihenfolgen sind wesentlich. Zuerst wird der Beitrag von \(A_{k-1}\) aufsummiert, danach entsteht durch den neuen Faktor \(A_k\). Innerhalb dieser Multiplikation werden die Grade absteigend aktualisiert. Aufsteigend würde bereits das veränderte \(a_{d-1}\) verwendet. Ausgehend von \(1\) entstünde dann fälschlich ein Term \(x^2/k^4\), obwohl nur mit \(1-x/k^2\) multipliziert wurde und das Ergebnis Grad eins haben muss.
Formal rechnen wir im Quotientenring \(\mathbb F_p[x]/(x^{10})\). Die Reduktion modulo \(x^{10}\) ist mit Addition und Multiplikation verträglich. Gleichbedeutend hängt der Koeffizient vom Grad \(d\le9\) eines Produkts nur von Eingangskoeffizienten bis Grad \(d\) ab. Kein verworfener Term höheren Grades beeinflusst daher einen gespeicherten Koeffizienten:
\[[x^d](F(x)G(x))=\sum_{r=0}^{d}([x^r]F(x))([x^{d-r}]G(x)),\qquad 0\le d\le9.\]
Ein serieller Durchlauf speichert somit zehn Produktkoeffizienten, aber nur den Skalar \(s_9\). Die übrigen Summenkoeffizienten gehen nicht in \(s_9\) ein. Alle zehn Produktkoeffizienten sind dennoch nötig: Wiederholte Aktualisierungen übertragen Informationen vom konstanten Term über die Zwischenstufen bis Grad neun. Außerdem gilt \(a_9^{(k-1)}=0\) für \(k<10\); die ersten neun Summanden tragen deshalb nichts zum Interpolationsanteil des gesuchten Koeffizienten bei.
Wähle \(N=3\) und behalte zunächst rationale Brüche bei, damit die Algebra gut lesbar bleibt. Anschließend lassen sich diese Identitäten modulo jeder Primzahl \(p>6\) reduzieren. Die ersten normierten Produkte sind
\[A_0=1,\qquad A_1=1-x,\qquad A_2=(1-x)(1-x/4)=1-\frac54x+\frac14x^2.\]
Die Gewichte für \(k=1,2,3\) sind \(1,1/6,1/15\). Jedes Gewicht verwendet das Produkt, bevor sein eigener Faktor eingefügt wurde. Einsetzen und Zusammenfassen der Koeffizienten ergibt
\[\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.\]
Damit sind sowohl die vorgeschriebenen Werte als auch die Formel \(d_3=1/60\) überprüft. Zugleich wird der Unterschied zwischen Interpolation und Normiertheit sichtbar: Beim Zielmodul besitzt dieses kleine Interpolationspolynom den Leitkoeffizienten \(1/60\ne1\). Interpolation allein ist also noch nicht die vollständige Konstruktion. Die Korrektur zur Normiertheit wird unten hergeleitet.
Sei \(u_k=k(2k-1)\). Alle \(u_k\) sind in \(\mathbb F_p\) von null verschieden, denn \(1\le k\le N<p\) und \(2k-1<p\). Für einen Block \(u_1,\ldots,u_m\) bilde die Präfixprodukte \(U_0=1\), \(U_i=U_{i-1}u_i\). Nur \(U_m^{-1}=U_m^{p-2}\) wird potenziert; danach läuft man rückwärts:
\[u_i^{-1}=U_{i-1}U_i^{-1}, \qquad U_{i-1}^{-1}=U_i^{-1}u_i.\]
Eine einzige Fermat-Potenzierung liefert so alle Kehrwerte eines Blocks. Außerdem ist \(k^{-1}=(2k-1)u_k^{-1}\); sein Quadrat ergibt \(k^{-2}\) für die Produktaktualisierung. Die Implementierungen verwenden Blöcke mit \(32768\) Einträgen.
Die Indizes der Batch-Formel sind lokal; der erste globale Index des Batches muss nicht \(1\) sein. Im Rückwärtsdurchlauf ist die laufende Inverse stets \(U_i^{-1}\). Multiplikation mit dem noch unveränderten Präfix \(U_{i-1}\) liefert \(u_i^{-1}\), und Multiplikation der laufenden Inversen mit \(u_i\) bereitet den nächsten Schritt vor. Der Präfixplatz darf anschließend durch die berechnete Inverse überschrieben werden, da die weiteren Schritte nur frühere Plätze benötigen.
Als konkreten Test dieses allgemeinen Tricks invertiere \(2,3,5\) modulo \(11\). Die Präfixprodukte sind \(1,2,6,8\), und \(8^{-1}=7\). Rückwärts erhält man zuerst \(5^{-1}=6\cdot7=9\), dann \(3^{-1}=2\cdot2=4\) und schließlich \(2^{-1}=1\cdot6=6\), jeweils modulo \(11\). Nur eine Potenzierung war nötig. Im eigentlichen Algorithmus sind die Eingaben die von null verschiedenen Werte \(k(2k-1)\); jede Multiplikation wird modulo \(p\) reduziert.
Da jedes Newton-Basispolynom normiert ist, lautet der Leitkoeffizient von \(P_N\)
\[c_N=d_N= \frac{(-1)^{N-1}N}{(2N-1)(N!)^2}.\]
Das an allen Stützstellen verschwindende Polynom ist
\[V_N(x)=\prod_{j=0}^{N}(x-j^2) =(-1)^N(N!)^2xA_N(x).\]
Die Formel für \(d_N\) ist von null verschieden: Ihr Zähler ist ein Vorzeichen, und kein Nennerfaktor verschwindet. Somit hat \(P_N\) genau Grad \(N\), nicht nur höchstens Grad \(N\). Zur Beschreibung aller Lösungen sei \(F\) ein beliebiges Polynom mit den vorgeschriebenen Werten. Jeder lineare Faktor \(x-j^2\) teilt \(F-P_N\); wegen der verschiedenen Nullstellen teilt auch ihr Produkt diese Differenz. Also
\[F(x)=P_N(x)+V_N(x)R(x),\qquad R(x)\in\mathbb F_p[x].\]
Für \(R\ne0\) gilt \(\deg(V_NR)=N+1+\deg R>\deg P_N\); der höchste Term kann sich also nicht wegheben. Die einzige Lösung vom Grad höchstens \(N\) ist daher \(P_N\). Ist sie nicht normiert, muss für eine normierte Lösung des nächsthöheren Grades \(R\) konstant sein, und ihr Leitkoeffizient erzwingt genau die Konstante \(1\). Damit sind Minimalität und Eindeutigkeit bewiesen. Einfaches Dividieren von \(P_N\) durch \(c_N\) würde dagegen den geforderten Wert bei \(1\) von \(1\) auf \(c_N^{-1}\) ändern; Skalieren ist deshalb keine zulässige Korrektur.
Ist \(c_N=1\), so ist das Interpolationspolynom selbst bereits die eindeutige normierte Lösung kleinsten Grades. Andernfalls kann es keine andere Lösung vom Grad höchstens \(N\) geben; der kleinste mögliche normierte Grad ist dann \(N+1\), und das Addieren des normierten Polynoms \(V_N\) liefert
\[Q(x)=P_N(x)+V_N(x).\]
Für die Zielinstanz prüft der Code \(c_N\ne1\). Der gesuchte Koeffizient folgt daher, ohne eines der großen Polynome zu konstruieren:
\[\boxed{[x^{10}]Q=s_9^{(N)} +(-1)^N(N!)^2a_9^{(N)}}.\]
Hinweis zum veröffentlichten kleinen Beispiel. Wörtlich gelesen besitzt die Formulierung „kleinster Grad“ bei \((N,M)=(2,7)\) einen Sonderfall: Dort ist \(c_2=1\), also ist \(P_2(x)=x^2\) bereits normiert und erfüllt alle drei Werte. Das angegebene kubische Polynom ist \(P_2+V_2\) und erfüllt sie ebenfalls. Dieser redaktionelle Randfall berührt die Zielinstanz nicht, denn dort gilt \(c_N\ne1\).
Als separates Beispiel ohne Sonderfall der Normiertheit wähle \((N,p)=(2,11)\). Das rationale Interpolationspolynom \((7x-x^2)/6\) reduziert sich zu \(3x+9x^2\). Das Addieren des verschwindenden Polynoms ergibt das folgende kubische Polynom. Seine Werte bei \(0,1,4\) sind modulo \(11\) gleich \(0,1,2\), und sein Leitkoeffizient ist \(1\). So verändert die Korrektur die Koeffizienten, ohne einen vorgeschriebenen Wert zu ändern.
\[\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}\]
Das C++-Programm kann aufeinanderfolgende Indexintervalle parallel bearbeiten. Für einen Block \(B\) setze \(g_k(x)=1-x/k^2\) und speichere bis Grad \(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.\]
Steht Block \(L\) unmittelbar vor Block \(R\), so werden ihre Zusammenfassungen kombiniert durch
\[\boxed{G_{LR}=G_LG_R,\qquad H_{LR}=H_L+G_LH_R}.\]
Diese Operation ist assoziativ, aber nicht vertauschbar; der Code führt die Blöcke in aufsteigender Indexreihenfolge zusammen. Gleichzeitig werden ihre Fakultätsprodukte multipliziert. Abgeschnittene Polynomprodukte genügen, weil Grade oberhalb von \(9\) keinen Beitrag zu kleineren Graden leisten können.
Die Blockregel beschreibt auch die Aktualisierung des eingehenden globalen Zustands: \((S,A)\mapsto(S+AH_B,AG_B)\). Vor jedem Beitrag des rechten Blocks müssen alle Faktoren des linken stehen; daher entsteht \(G_LH_R\). Bei drei aufeinanderfolgenden Blöcken ergibt jede Klammerung genau den folgenden Ausdruck, womit die Assoziativität bewiesen ist. Der leere Block besitzt die Zusammenfassung \((G,H)=(1,0)\).
\[G_{LMR}=G_LG_MG_R,\qquad H_{LMR}=H_L+G_LH_M+G_LG_MH_R.\]
Deshalb muss ein Worker alle zehn Koeffizienten von \(H_B\) behalten, obwohl am Ende nur Grad neun benötigt wird. Ein vorangehendes Produkt mischt die Grade gemäß der folgenden Faltung. Der Merge-Code darf nur den neunten Koeffizienten der globalen Summe speichern, nicht aber nur den neunten Koeffizienten jeder Worker-Summe. Inversionsbatches und parallele Blöcke haben unterschiedliche Aufgaben: Ein Worker kann innerhalb seines zusammenhängenden Intervalls viele Inversionsbatches verarbeiten.
\[[x^9](AH_B)=\sum_{r=0}^{9}([x^r]A)([x^{9-r}]H_B).\]
Lemma 1. Die Quadratstellen sind in \(\mathbb F_p\) verschieden. Somit existiert \(P_N\) und ist das einzige Interpolationspolynom vom Grad höchstens \(N\).
Lemma 2. Die Auswertung des Nenners der dividierten Differenz und der alternierenden Binomialsumme ergibt das angegebene \(d_k\). Eingesetzt in Newtons Formel folgt \(P_N/x=\sum A_{k-1}/(k(2k-1))\).
Lemma 3. Die absteigende Rekursion erhält die Invariante \(a_d^{(k)}=[x^d]A_k\); die Akkumulation erhält \(s_d^{(k)}=[x^d]\sum_{i=1}^{k}A_{i-1}/(i(2i-1))\). Der Abbruch bei Grad \(9\) bewahrt beide Invarianten exakt.
Lemma 4. Die Batch-Inversion liefert den tatsächlichen Kehrwert jedes \(u_k\ne0\), denn \(U_{i-1}U_i^{-1}=1/u_i\). Sie ändert daher nur die Auswertungsweise derselben Rekursion.
Lemma 5. Jedes Polynom mit den vorgegebenen Werten unterscheidet sich von \(P_N\) um ein Vielfaches von \(V_N\). Der Leitkoeffiziententest wählt also genau dann \(P_N\), wenn es bereits normiert ist, und sonst das normierte Polynom \(P_N+V_N\) vom kleinsten möglichen Grad.
Lemma 6. Die geordnete Blockidentität setzt \(G_L\) vor jeden Beitrag aus \(R\). Induktion über die zusammengeführten Blöcke reproduziert daher exakt die serielle Summe und das serielle Produkt.
Zusammen zeigen die Lemmata, dass der Endausdruck genau der Koeffizient von \(x^{10}\) in der verlangten funktionalen Inversen kleinsten normierten Grades ist.
In C++ stellt Block.product[d] den Koeffizienten des abgeschnittenen Produkts und Block.sum[d] die entsprechende gewichtete Summe dar. Die serielle Spezialisierung von calculate_block behält nur sum[9]; parallele Worker speichern die ganze Summe. Python und Java halten dasselbe Produkt in coefficients und den neunten Summenkoeffizienten in interpolation_coefficient beziehungsweise interpolationCoefficient. Der mathematische Zustand ist in allen drei Implementierungen gleich.
Zu Beginn der seriellen Iteration \(k\) enthält das Array \(A_{k-1}\), der Akkumulator \(s_9^{(k-1)}\) und die Fakultätsvariable \((k-1)!\). Lies das Batch-Gewicht, addiere seinen Beitrag, aktualisiere das Array absteigend und multipliziere die Fakultätsvariable mit \(k\). Dadurch gilt dieselbe Invariante für Iteration \(k+1\). Nach der letzten Iteration ist der Zustand genau \((A_N,s_9^{(N)},N!)\), also alles, was die Endformel benötigt.
C++ verwendet standardmäßig höchstens \(16\) Worker-Threads; das optionale Argument --threads kann eine andere unterstützte Anzahl wählen. Die Worker bearbeiten zusammenhängende Intervalle und werden in ihrer ursprünglichen Reihenfolge zusammengeführt. Python und Java werten die serielle Rekurrenz aus. Binäre modulare Potenzierung erfolgt mit C++ power, Pythons dreistelligem pow und Java modPow. Bei diesem Modul ist das Produkt zweier reduzierter Reste kleiner als \((10^9+7)^2<2^{63}\); C++-64-Bit-Arithmetik und Java long genügen daher, ein 32-Bit-Produkt jedoch nicht.
Zur Überprüfung enthält jede Implementierung außerdem eine direkte Lagrange-Interpolation, die für \(1\le N\le40\) das vollständige Polynom konstruiert. Sie verwendet eine andere Konstruktion als die optimierte Newton-Rekurrenz. Die Vergleiche mit dieser Referenz umfassen den bereits normierten Fall \(N=1\) und prüfen den Koeffizienten vom Grad zehn auch dann, wenn das kleine Polynom keinen solchen Term hat. Die C++-Selbsttests vergleichen zusätzlich serielle und parallele Ergebnisse an Inversionsbatch-Grenzen; Python und Java testen auch eine bewusst kleine Batch-Größe.
Sei \(D=9\), \(B=32768\) die Größe eines Inversionsbatches und \(T\) die Worker-Anzahl, seriell also \(T=1\). Jeder Index benötigt \(D\) Aktualisierungen der Produktkoeffizienten und konstant viel weitere Arithmetik; eine vollständige Worker-Summe ergänzt \(D+1\) Aktualisierungen. Präfix- und Rückwärtsdurchläufe sind linear in der Batch-Größe. Insgesamt gibt es höchstens \(\lceil N/B\rceil+T-1\) Batch-Inversionen und eine weitere Inversion für den Leitkoeffiziententest.
\[W(N,D,B,T,p)=O\!\left(ND+\left(\frac NB+T\right)\log p+TD^2\right).\]
Der Term \(TD^2\) erfasst die kurzen Polynomfaltungen beim Zusammenführen. Für feste \(D,B,T\) und diesen festen Modul ist die gesamte Arbeit \(O(N)\). Beim Zielwert \(N\) führt ein serieller Durchlauf \(90{,}000{,}000\) Produktkoeffizienten-Aktualisierungen und nur \(306\) Batch-Inversionen aus, zuzüglich der letzten Inversion für den Leitkoeffizienten. Parallelisierung verteilt die Hauptschleife auf Worker, ändert aber nicht die Größenordnung der Gesamtarbeit; die tatsächliche Laufzeit hängt auch von Scheduling und Hardware ab.
Der zusätzliche Speicher beträgt seriell \(O(B+D)\) und parallel \(O(T(B+D))\), ohne die Thread-Stacks der Laufzeitumgebung. Es gibt weder ein Array für alle \(N\) Stützstellen noch ein Koeffizientenarray bis Grad \(N\). Dieselbe Herleitung kann einen anderen positiven Grad \(r\) berechnen, indem die Grade \(0\) bis \(r-1\) behalten werden. Der gesuchte Grad zehn begrenzt den vorliegenden Zustand auf zehn Koeffizienten.
Asal bir \(M\) modülü için \(x^2\)'nin \((N,M)\)-fonksiyonel tersi; katsayıları \(\{0,\ldots,M-1\}\) içinde temsil edilen, baş katsayısı \(1\) olan en düşük dereceli \(Q(x)\) polinomudur ve
\[Q(n^2)\equiv n\pmod M\qquad(0\le n\le N).\]
koşulunu sağlar. Bu soruda \(N=10^7\), \(M=p=10^9+7\) ve yalnızca \(x^{10}\)'un katsayısı isteniyor. Derecesi \(N\) olan enterpolasyon polinomunu kurmak, hatta tüm katsayılarını saklamak gerekmiyor: hesap tam doğrulukla \(9\). derecede kesilebilir.
\([x^r]F(x)\), \(F\) polinomundaki \(x^r\) katsayısını gösterir. Açıkça rasyonel özdeşlik olarak belirtilmedikçe aşağıdaki tüm eşitlikler \(\mathbb F_p\) içinde geçerlidir. \(a/b\) kesri, tam sayı bölmesi veya kayan noktalı bölme değil, \(a\,b^{-1}\pmod p\) anlamına gelir. \(p\) asal olduğundan sıfırdan farklı her paydanın tek bir tersi vardır. Sonucun \(0\) ile \(p-1\) arasındaki temsilcisini seçmek, sorunun istediği katsayı biçimini verir.
İstenen, kare alma işlemini bütün cisim üzerinde tersine çevirmek değildir: \(n\) ve \(-n\) aynı kareye sahiptir. Yalnızca seçilen kare düğümlerinin her birinde bir değer belirlenmiştir. \(2N<p\) sınırı, seçilen iki pozitif indisin aynı düğümün birbirinin negatifi olan kareköklerini vermesini engeller. Bu nedenle doğru başlangıç noktası, gerçek karekök fonksiyonunun kuvvet serisi değil, interpolasyondur.
\(\mathbb F_p\) cismi üzerinde çalışalım. \(0^2,1^2,\ldots,N^2\) düğümleri birbirinden farklıdır. Çünkü \(i^2\equiv j^2\pmod p\) olsaydı \(p\mid(i-j)(i+j)\) gerekirdi; oysa \(0\le i,j\le N\) ve \(i+j\le2N<p\). Dolayısıyla bütün düğümlerde \(P_N(n^2)=n\) sağlayan, derecesi en çok \(N\) olan tek bir \(P_N\) polinomu vardır.
Bu polinomu Newton biçiminde yazalım:
\[P_N(x)=\sum_{k=1}^{N}d_k\prod_{j=0}^{k-1}(x-j^2).\]
\(0\)'da istenen değer \(0\) olduğu için \(k=0\) terimi yoktur.
Varlık ve teklik açıkça şöyle gösterilir. Sabit bir \(k\) için aşağıdaki Lagrange bazını tanımlayalım. \(x=j^2\) noktasında \(i\). baz polinomu, \(i=j\) ise \(1\), aksi hâlde \(0\) olur. Dolayısıyla ağırlıklı toplamları istenen değerleri verir. Derecesi en fazla \(k\) olan iki polinom bu değerleri verseydi farklarının \(k+1\) farklı kökü olurdu; sıfır olmayan ve derecesi en fazla \(k\) olan bir polinom için bu imkânsızdır.
\[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 biçimi, düğümleri birer birer ekler. \(P_k-P_{k-1}\) farkı ilk \(k\) düğümde sıfır olduğundan bu düğümlerin monik çarpımı \(B_k\)'nın sabit katıdır. Yeni düğümde değerini hesaplamak bu sabiti belirler. \(P_{k-1}\)'in derecesi en fazla \(k-1\) olduğu için \(d_k\), aynı zamanda \(P_k\)'nın Lagrange ifadesindeki \(x^k\) katsayısıdır. Bu baş katsayıyı okumak, sıradaki bölünmüş fark toplamını verir.
\[\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}\]
\(0^2,\ldots,k^2\) üzerindeki bölünmüş fark
\[d_k=\sum_{i=1}^{k}\frac{i} {\displaystyle\prod_{\substack{0\le j\le k\\j\ne i}}(i^2-j^2)}.\]
Paydayı hesaplamak için iki doğrusal çarpımı ayrı ele alalım. İlk çarpımda tam \(k-i\) tane negatif çarpan vardır. İkincisinde ise \(i\)'den \(i+k\)'ya kadar uzanan ardışık çarpımdan, atlanan \(2i\) çarpanı çıkarılmalıdır. \(1\le i\le k\) için böylece
\[\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}\]
Bu ifadeleri çarparken \(i!/(i-1)!=i\) ile \(2i\) sadeleşir. Ortaya çıkan tüm faktöriyeller mod \(p\) terslenebilir; çünkü argümanları en fazla \(2k\le2N<p\)'dir. Dolayısıyla
\[\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)!}.\]
Kalan toplamı, sadeleşmenin hiçbir adımını atlamadan hesaplayabiliriz. \((2k)!\) ile çarpıp \(r=k-i\) yazalım ve sonucu \(S_k\) olarak adlandıralım. Faktöriyel paydası böylece bir binom katsayısına dönüşür:
\[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).\]
Değişken işaretli toplamlar için her terime Pascal özdeşliğini uygulayalım. İkinci toplamın indisini bir kaydırınca iç terimlerin hepsi birbirini götürür ve yalnızca son terim kalır. Alt indisi \(-1\) olan binom katsayısını ve boş toplamı sıfır kabul ediyoruz.
\[\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}\]
Binom katsayılarının farkını \(S_k\)'da yerine koyalım. İndis kaydırılırken değişken işaret de değiştiğinden farktaki eksi işareti artıya dönüşür. Kısmi toplam özdeşliğini iki kez uygulayıp ardından komşu iki binom katsayısının oranını kullanalım. Boş toplam kuralıyla hesap \(k=1\) için de geçerlidir:
\[\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)!}}.\]
Şu polinomları tanımlayalım:
\[A_k(x)=\prod_{j=1}^{k}\left(1-\frac{x}{j^2}\right), \qquad A_0(x)=1.\]
Çünkü
\[\prod_{j=0}^{k-1}(x-j^2) =x(-1)^{k-1}((k-1)!)^2A_{k-1}(x),\]
işaretler ve faktöriyeller \(d_k\) ile sadeleşir ve şu yalın özdeşlik kalır:
\[\boxed{\frac{P_N(x)}x= \sum_{k=1}^{N}\frac{A_{k-1}(x)}{k(2k-1)}}.\]
Bu nedenle \([x^{10}]P_N=[x^9](P_N/x)\). \(a_d^{(k)}=[x^d]A_k\) olsun; \(s_d^{(k)}\) de \(1,\ldots,k\) terimlerinden sonra biriken katsayıyı göstersin. \(1\le d\le9\) için her adım
\[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}.\]
biçimindedir. \(a_d\) dizisi \(9\)'dan \(1\)'e doğru güncellenir; böylece sağ taraftaki katsayı hâlâ \(A_{k-1}\)'e aittir. \(9\)'un üzerindeki dereceleri atmak bir yaklaşım değil, tam bir işlemdir: sonraki çarpanların dereceleri negatif olamaz ve atılmış bir terim tekrar \(9\). dereceye inemez.
Daha açık yazarsak \(A_{k-1}\)'in önündeki çarpan \(d_k(-1)^{k-1}((k-1)!)^2=1/(k(2k-1))\) olur. Katsayı durumunu ve toplamı aşağıdaki gibi başlatırız. \(A_k\)'nın sabit katsayısı her zaman \(1\) kaldığından bu katsayının güncellenmesi gerekmez.
\[a_0^{(0)}=1,\qquad a_d^{(0)}=0\ (1\le d\le9),\qquad s_d^{(0)}=0\ (0\le d\le9).\]
İki ayrı işlem sırası önemlidir. Önce \(A_{k-1}\)'in katkısı toplama eklenir, ardından yeni çarpanla çarpılarak \(A_k\) elde edilir. Bu çarpım içinde dereceler azalan sırada güncellenir. Artan sırada güncelleme, değiştirilmiş \(a_{d-1}\)'i kullanırdı. Örneğin \(1\)'den başlarken yalnızca \(1-x/k^2\) ile çarpmak, sonuç birinci dereceden olması gerekirken yanlışlıkla \(x^2/k^4\) terimi üretirdi.
Kesmenin güvenli olmasının biçimsel nedeni, hesabı \(\mathbb F_p[x]/(x^{10})\) bölüm halkasında yapmamızdır. \(x^{10}\)'a göre indirgeme, toplama ve çarpmayla uyumludur. Eşdeğer olarak bir çarpımın \(d\le9\) dereceli katsayısı, yalnızca girdilerin en fazla \(d\) dereceli katsayılarını kullanır. Atılan hiçbir yüksek dereceli terim saklanan katsayıları etkileyemez:
\[[x^d](F(x)G(x))=\sum_{r=0}^{d}([x^r]F(x))([x^{d-r}]G(x)),\qquad 0\le d\le9.\]
Bu nedenle seri hesap, çarpımın on katsayısını fakat toplamın yalnızca \(s_9\) skalerini saklar. Toplamın diğer katsayıları \(s_9\)'u beslemez. Buna karşılık çarpımın on katsayısı da gereklidir: ardışık güncellemeler, bilgiyi sabit terimden ara dereceler üzerinden dokuzuncu dereceye taşır. Ayrıca \(k<10\) için \(a_9^{(k-1)}=0\) olduğundan ilk dokuz terim, aranan katsayının interpolasyon kısmına katkı yapmaz.
\(N=3\) alalım ve cebirsel işlemleri açık görmek için geçici olarak rasyonel kesirlerle çalışalım. Bu özdeşlikler daha sonra herhangi bir \(p>6\) asalı için modüler olarak indirgenebilir. İlk normalize edilmiş çarpımlar şunlardır:
\[A_0=1,\qquad A_1=1-x,\qquad A_2=(1-x)(1-x/4)=1-\frac54x+\frac14x^2.\]
\(k=1,2,3\) için ağırlıklar \(1,1/6,1/15\)'tir. Her ağırlık, kendi çarpanı eklenmeden önceki çarpımı kullanır. Yerine koyup katsayıları toplarsak
\[\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.\]
Bu hesap hem istenen değerleri hem de \(d_3=1/60\) katsayı formülünü doğrular. Ayrıca interpolasyon ile monik koşulunun ayrı olduğunu gösterir: hedef modül altında bu küçük interpolasyon polinomunun baş katsayısı \(1/60\ne1\)'dir. Dolayısıyla yalnızca interpolasyon yapmak, yapıyı tamamlamaz. Monik düzeltmeyi aşağıda türeteceğiz.
\(u_k=k(2k-1)\) diyelim. \(1\le k\le N<p\) ve \(2k-1<p\) olduğundan her \(u_k\), \(\mathbb F_p\)'de sıfırdan farklıdır. Bir \(u_1,\ldots,u_m\) partisi için \(U_0=1\), \(U_i=U_{i-1}u_i\) önek çarpımları kurulur. Yalnızca \(U_m^{-1}=U_m^{p-2}\) hesaplanır ve dizi geriye doğru taranır:
\[u_i^{-1}=U_{i-1}U_i^{-1}, \qquad U_{i-1}^{-1}=U_i^{-1}u_i.\]
Böylece tek bir Fermat üs alması bütün partinin terslerini verir. Ayrıca \(k^{-1}=(2k-1)u_k^{-1}\) olduğundan bunun karesi, çarpım güncellemesinde gereken \(k^{-2}\)'yi verir. Uygulamalar \(32768\) elemanlı partiler kullanır.
Toplu ters alma formülündeki indisler o gruba aittir; grubun ilk genel indisi \(1\) olmak zorunda değildir. Geriye doğru geçişte, taşınan tersin \(U_i^{-1}\) olmasını değişmez olarak koruruz. Henüz değiştirilmemiş \(U_{i-1}\) önekiyle çarpmak \(u_i^{-1}\)'i verir; taşınan tersi \(u_i\) ile çarpmak ise sonraki adıma hazırlar. Artık önek hücresinin üzerine bulunan ters yazılabilir, çünkü sonraki adımlar yalnızca daha önceki hücreleri kullanacaktır.
Bu genel yöntemi somut olarak görmek için \(2,3,5\) sayılarının mod \(11\) terslerini bulalım. Önek çarpımları \(1,2,6,8\)'dir ve \(8^{-1}=7\) olur. Geriye doğru geçiş önce \(5^{-1}=6\cdot7=9\), sonra \(3^{-1}=2\cdot2=4\), son olarak \(2^{-1}=1\cdot6=6\) verir; bütün eşitlikler mod \(11\)'dedir. Yalnızca bir modüler üs alma yeterli olmuştur. Asıl algoritmada girdiler sıfırdan farklı \(k(2k-1)\) değerleridir ve her çarpım mod \(p\) indirgenir.
Her Newton taban polinomu monik olduğundan \(P_N\)'in baş katsayısı
\[c_N=d_N= \frac{(-1)^{N-1}N}{(2N-1)(N!)^2}.\]
olur. Bütün enterpolasyon düğümlerinde sıfırlanan polinom ise
\[V_N(x)=\prod_{j=0}^{N}(x-j^2) =(-1)^N(N!)^2xA_N(x).\]
\(d_N\) formülü sıfırdan farklıdır: payı bir işarettir ve paydadaki hiçbir çarpan sıfır olmaz. Bu nedenle \(P_N\)'in derecesi yalnızca en fazla \(N\) değil, tam olarak \(N\)'dir. Tüm çözümleri belirlemek için istenen değerleri veren herhangi bir \(F\) polinomu alalım. Her \(x-j^2\) doğrusal çarpanı \(F-P_N\)'i böler; kökler farklı olduğundan bunların çarpımı da bu farkı böler. Yani
\[F(x)=P_N(x)+V_N(x)R(x),\qquad R(x)\in\mathbb F_p[x].\]
\(R\ne0\) ise \(\deg(V_NR)=N+1+\deg R>\deg P_N\) olur; en yüksek dereceli terim yok olamaz. Dolayısıyla derecesi en fazla \(N\) olan tek çözüm \(P_N\)'dir. Bu çözüm monik değilse bir sonraki derecede monik bir çözüm elde etmek için \(R\) sabit olmalı; baş katsayı koşulu da bu sabiti tam olarak \(1\) yapmalıdır. Böylece hem en küçük derece hem de teklik kanıtlanır. \(P_N\)'i doğrudan \(c_N\)'e bölmek ise \(1\) noktasındaki gerekli değeri \(1\)'den \(c_N^{-1}\)'e değiştirir; bu yüzden ölçekleme geçerli bir düzeltme değildir.
biçimindedir. \(c_N=1\) ise enterpolant zaten en düşük dereceli tek monik çözümdür. Değilse, derecesi en çok \(N\) olan başka bir çözüm bulunamaz; mümkün olan en küçük monik derece \(N+1\)'dir ve monik \(V_N\)'i eklemek
\[Q(x)=P_N(x)+V_N(x).\]
çözümünü verir. Hedef örnek için kod \(c_N\ne1\) olduğunu denetler. Dolayısıyla aranan katsayı iki büyük polinom da oluşturulmadan hesaplanır:
\[\boxed{[x^{10}]Q=s_9^{(N)} +(-1)^N(N!)^2a_9^{(N)}}.\]
Yayımlanmış küçük örnek hakkında not. “En düşük derece” ifadesi kelimesi kelimesine okunduğunda \((N,M)=(2,7)\) için özel bir kenar durum vardır: \(c_2=1\) olduğundan \(P_2(x)=x^2\) zaten moniktir ve üç değeri de sağlar. Gösterilen kübik polinom \(P_2+V_2\)'dir ve o da değerleri sağlar. Bu editoryal kenar durum hedef hesabı etkilemez; hedefte \(c_N\ne1\)'dir.
Monik olma istisnasının bulunmadığı ayrı bir örnek için \((N,p)=(2,11)\) alalım. Rasyonel interpolasyon polinomu \((7x-x^2)/6\), \(3x+9x^2\)'ye indirgenir. Sıfırlayan polinomu ekleyince aşağıdaki kübik polinom elde edilir. Bu polinomun \(0,1,4\) noktalarındaki değerleri mod \(11\) sırasıyla \(0,1,2\)'dir ve baş katsayısı \(1\)'dir. Bu örnek, düzeltmenin katsayıları değiştirirken gerekli bütün değerleri nasıl koruduğunu gösterir.
\[\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}\]
C++ programı ardışık indis aralıklarını paralel işleyebilir. Bir \(B\) bloğu için \(g_k(x)=1-x/k^2\) yazıp \(9\). dereceye kadar
\[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.\]
özetlerini saklayalım. \(L\) bloğunun hemen ardından \(R\) bloğu geliyorsa birleşim
\[\boxed{G_{LR}=G_LG_R,\qquad H_{LR}=H_L+G_LH_R}.\]
özdeşliğiyle yapılır. İşlem birleşmelidir fakat değişmeli değildir; uygulama blokları artan indis sırasında birleştirir. Faktöriyel çarpımları da aynı anda çarpılır. \(9\). derecede kesilmiş polinom çarpımı yeterlidir, çünkü daha yüksek dereceler daha düşük derecelere katkı yapamaz.
Blok kuralı, gelen genel durumun \((S,A)\mapsto(S+AH_B,AG_B)\) biçiminde güncellenmesi olarak da okunabilir. Sağ bloktan gelen her katkının önünde sol bloğun bütün çarpanları bulunmalıdır; \(G_LH_R\) terimi buradan gelir. Ardışık üç blokta iki parantezleme de aşağıdaki aynı ifadeyi verir; bu, birleşme özelliğini kanıtlar. Boş bloğun özeti \((G,H)=(1,0)\)'dır.
\[G_{LMR}=G_LG_MG_R,\qquad H_{LMR}=H_L+G_LH_M+G_LG_MH_R.\]
Bu durum, sonuçta yalnızca dokuzuncu derece gerekmesine rağmen bir iş parçacığının neden \(H_B\)'nin on katsayısını da tutması gerektiğini açıklar. Önceki çarpım, aşağıdaki konvolüsyonla dereceleri karıştırır. Birleştirme kodu genel toplamın yalnızca dokuzuncu katsayısını saklayabilir, ancak her iş parçacığının toplamı için aynısını yapamaz. Toplu ters alma grupları ile paralel bloklar farklı amaçlara hizmet eder: bir iş parçacığı kendi ardışık aralığı içinde birçok ters alma grubu işleyebilir.
\[[x^9](AH_B)=\sum_{r=0}^{9}([x^r]A)([x^{9-r}]H_B).\]
Lemma 1. Kare düğümler \(\mathbb F_p\)'de farklıdır; dolayısıyla \(P_N\) vardır ve derecesi en çok \(N\) olan tek enterpolanttır.
Lemma 2. Bölünmüş fark paydasının ve almaşık binom toplamının değerlendirilmesi belirtilen \(d_k\)'yi verir. Bu sonuç Newton formuna konunca \(P_N/x=\sum A_{k-1}/(k(2k-1))\) elde edilir.
Lemma 3. Azalan derece sırasındaki yineleme \(a_d^{(k)}=[x^d]A_k\) değişmezini; toplama adımı da \(s_d^{(k)}=[x^d]\sum_{i=1}^{k}A_{i-1}/(i(2i-1))\) değişmezini korur. \(9\). derecede kesme bunları tam olarak korur.
Lemma 4. \(U_{i-1}U_i^{-1}=1/u_i\) olduğundan toplu ters alma, sıfırdan farklı her \(u_k\)'nin gerçek tersini üretir. Yalnızca aynı yinelemenin değerlendirilme biçimi değişir.
Lemma 5. Verilen değerleri sağlayan her polinomun \(P_N\)'den farkı \(V_N\)'in bir katıdır. Baş katsayı testi bu nedenle, yalnızca zaten monikse \(P_N\)'i; aksi hâlde mümkün olan en düşük dereceli monik \(P_N+V_N\)'i seçer.
Lemma 6. Sıralı blok özdeşliği, \(R\)'den gelen her katkının önüne \(G_L\)'yi yerleştirir. Birleştirilen bloklar üzerinde tümevarım, seri toplamı ve çarpımı aynen yeniden üretir.
Bu lemmalar birlikte, son ifadenin istenen en düşük dereceli monik fonksiyonel tersin \(x^{10}\) katsayısı olduğunu kanıtlar.
C++ kodunda Block.product[d], kesilmiş çarpımın katsayısını; Block.sum[d] ise ilgili ağırlıklı toplamı temsil eder. calculate_block'un seri özelleştirmesi yalnızca sum[9]'u tutar; paralel çalışan iş parçacıkları toplamın tamamını saklar. Python ve Java aynı çarpımı coefficients dizisinde, toplamın dokuzuncu katsayısını ise sırasıyla interpolation_coefficient ve interpolationCoefficient değişkenlerinde tutar. Üç uygulamada da matematiksel durum aynıdır.
Seri döngünün \(k\). adımı başlarken dizi \(A_{k-1}\)'i, biriktirici \(s_9^{(k-1)}\)'i, faktöriyel değişkeni ise \((k-1)!\)'i içerir. Gruptan ağırlık alınır, katkısı biriktiriciye eklenir, dizi azalan sırada güncellenir ve faktöriyel değişkeni \(k\) ile çarpılır. Bu işlemler aynı değişmezi \(k+1\). adım için yeniden sağlar. Son adımın ardından durum tam olarak \((A_N,s_9^{(N)},N!)\)'dir; son formülün ihtiyacı olan bütün bilgi budur.
C++ varsayılan olarak en fazla \(16\) iş parçacığı kullanır; isteğe bağlı --threads argümanı desteklenen başka bir sayıyı seçebilir. İş parçacıkları ardışık aralıkları işler ve başlangıç sıralarında birleştirilir. Python ve Java seri bağıntıyı uygular. İkili modüler üs alma C++'ta power, Python'da üç argümanlı pow, Java'da modPow ile yapılır. Bu modülde iki indirgenmiş kalanın çarpımı \((10^9+7)^2<2^{63}\)'ten küçüktür; bu nedenle C++'ın 64 bit aritmetiği ve Java'nın long türü yeterlidir, fakat 32 bit çarpım yeterli olmaz.
Doğrulama için her uygulama ayrıca \(1\le N\le40\) aralığında polinomun tamamını kuran doğrudan bir Lagrange interpolasyonu içerir. Bu yöntem, optimize Newton bağıntısından farklı bir yapı kullanır. Referansla karşılaştırmalar \(N=1\)'de zaten monik olan durumu da kapsar ve küçük polinomda onuncu derece terimi bulunmadığında bile bu katsayıyı kontrol eder. C++ öz testleri ayrıca ters alma gruplarının sınırları çevresinde seri ve paralel sonuçları karşılaştırır; Python ve Java da özellikle küçük seçilmiş bir grup boyutunu sınar.
\(D=9\), ters alma grubu boyutu \(B=32768\) ve iş parçacığı sayısı \(T\) olsun; seri hesapta \(T=1\)'dir. Her indis, \(D\) çarpım katsayısı güncellemesi ve sabit miktarda ek aritmetik gerektirir; tam iş parçacığı toplamı için buna \(D+1\) güncelleme daha eklenir. Önek ve geriye doğru geçişlerin maliyeti grup boyutuyla doğrusaldır. Toplam grup tersi sayısı en fazla \(\lceil N/B\rceil+T-1\)'dir; baş katsayı testi bir ek ters alma kullanır.
\[W(N,D,B,T,p)=O\!\left(ND+\left(\frac NB+T\right)\log p+TD^2\right).\]
\(TD^2\) terimi, birleştirmelerdeki kısa polinom konvolüsyonlarını temsil eder. \(D,B,T\) ve bu modül sabitken toplam işlem miktarı \(O(N)\)'dir. Hedef \(N\) için seri hesap \(90{,}000{,}000\) çarpım katsayısı güncellemesi ve yalnızca \(306\) grup tersi hesaplar; son baş katsayı için de bir ters daha alınır. Paralellik ana döngüyü iş parçacıklarına böler, ancak toplam işlem miktarının mertebesini değiştirmez; gerçek süre zamanlamaya ve donanıma da bağlıdır.
Çalışma zamanının iş parçacığı yığınları hariç ek bellek, seri uygulamada \(O(B+D)\), paralel uygulamada \(O(T(B+D))\)'dir. Bütün \(N\) düğümü indisleyen bir dizi veya derece \(N\)'ye kadar katsayı dizisi bulunmaz. Aynı türetimle başka bir pozitif \(r\) derecesi de \(0\)'dan \(r-1\)'e kadar dereceler saklanarak hesaplanabilir. Bu soruda onuncu derecenin istenmesi, durumu yalnızca on katsayıyla sınırlamamızı sağlar.
Para un módulo primo \(M\), una inversa funcional \((N,M)\) de \(x^2\) es un polinomio mónico \(Q(x)\) de grado mínimo, con coeficientes representados en \(\{0,\ldots,M-1\}\), que cumple
\[Q(n^2)\equiv n\pmod M\qquad(0\le n\le N).\]
En este caso \(N=10^7\), \(M=p=10^9+7\), y solo se pide el coeficiente de \(x^{10}\). No hace falta construir el polinomio interpolador de grado \(N\), ni siquiera guardar todos sus coeficientes: el cálculo puede truncarse rigurosamente en grado \(9\).
La notación \([x^r]F(x)\) indica el coeficiente de \(x^r\) en \(F\). Todas las igualdades siguientes se interpretan en \(\mathbb F_p\), salvo que se presenten explícitamente como identidades racionales. Una fracción \(a/b\) significa \(a\,b^{-1}\pmod p\), no división entera ni de coma flotante. Como \(p\) es primo, cada denominador no nulo tiene una inversa única. Elegir el representante entre \(0\) y \(p-1\) proporciona el formato de coeficientes exigido.
No se pide invertir la función cuadrado en todo el cuerpo: \(n\) y \(-n\) tienen el mismo cuadrado. Solo se prescribe un valor en cada nodo cuadrado seleccionado. La cota \(2N<p\) impide que dos índices positivos seleccionados representen raíces cuadradas opuestas del mismo nodo. Por eso el punto de partida adecuado es la interpolación, no una serie de potencias de la raíz cuadrada real.
Trabajamos en el cuerpo \(\mathbb F_p\). Los nodos \(0^2,1^2,\ldots,N^2\) son distintos: si \(i^2\equiv j^2\pmod p\), entonces \(p\mid(i-j)(i+j)\), pero \(0\le i,j\le N\) e \(i+j\le2N<p\). Por tanto existe un único polinomio \(P_N\), de grado a lo sumo \(N\), tal que \(P_N(n^2)=n\) en todos los nodos requeridos.
Lo escribimos en forma de Newton:
\[P_N(x)=\sum_{k=1}^{N}d_k\prod_{j=0}^{k-1}(x-j^2).\]
El término \(k=0\) es nulo porque el valor prescrito en \(0\) es \(0\).
La existencia y la unicidad se justifican de forma explícita. Para un \(k\) fijo, definimos la siguiente base de Lagrange. En \(x=j^2\), su miembro \(i\) vale \(1\) si \(i=j\) y \(0\) en caso contrario. Su suma ponderada tiene, por tanto, los valores prescritos. Si dos polinomios de grado a lo sumo \(k\) tuvieran esos valores, su diferencia tendría \(k+1\) raíces distintas, algo imposible para un polinomio no nulo de ese grado.
\[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).\]
La forma de Newton añade los nodos uno a uno. La diferencia \(P_k-P_{k-1}\) se anula en los primeros \(k\) nodos, así que es un múltiplo constante de su producto mónico \(B_k\). Evaluar en el nuevo nodo determina esa constante. Como \(P_{k-1}\) tiene grado a lo sumo \(k-1\), \(d_k\) también es el coeficiente de \(x^k\) en la expresión de Lagrange de \(P_k\). Leer ese coeficiente principal proporciona la suma de diferencias divididas que sigue.
\[\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}\]
La diferencia dividida sobre \(0^2,\ldots,k^2\) es
\[d_k=\sum_{i=1}^{k}\frac{i} {\displaystyle\prod_{\substack{0\le j\le k\\j\ne i}}(i^2-j^2)}.\]
Para evaluar el denominador, separamos los dos productos lineales. El primero contiene exactamente \(k-i\) factores negativos. En el segundo hay que quitar el factor omitido \(2i\) del producto consecutivo desde \(i\) hasta \(i+k\). Para \(1\le i\le k\), obtenemos
\[\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}\]
Al multiplicar, \(i!/(i-1)!=i\) se cancela con el factor \(2i\). Todos los factoriales que aparecen son invertibles módulo \(p\), pues sus argumentos no superan \(2k\le2N<p\). Por tanto
\[\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)!}.\]
Podemos evaluar la suma restante sin omitir ninguna cancelación. Multiplicamos por \((2k)!\), ponemos \(r=k-i\) y llamamos \(S_k\) al resultado. El denominador factorial se convierte en un coeficiente binomial:
\[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).\]
Para las sumas alternadas, aplicamos la identidad de Pascal a cada término. Al desplazar un índice la segunda suma, se cancelan todos los términos interiores y solo queda el último. Adoptamos la convención de que un coeficiente binomial con índice inferior \(-1\) y una suma vacía valen cero.
\[\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}\]
Sustituimos la diferencia binomial en \(S_k\). Su signo menos se vuelve más al desplazar el índice, porque también cambia el signo alternado. Aplicamos dos veces la identidad de sumas parciales y después usamos el cociente de los dos coeficientes binomiales adyacentes. Con la convención de suma vacía, el cálculo también vale para \(k=1\):
\[\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)!}}.\]
Definimos
\[A_k(x)=\prod_{j=1}^{k}\left(1-\frac{x}{j^2}\right), \qquad A_0(x)=1.\]
Como
\[\prod_{j=0}^{k-1}(x-j^2) =x(-1)^{k-1}((k-1)!)^2A_{k-1}(x),\]
los signos y factoriales se cancelan con \(d_k\), dejando la identidad sencilla
\[\boxed{\frac{P_N(x)}x= \sum_{k=1}^{N}\frac{A_{k-1}(x)}{k(2k-1)}}.\]
Así, \([x^{10}]P_N=[x^9](P_N/x)\). Sea \(a_d^{(k)}=[x^d]A_k\), y sea \(s_d^{(k)}\) el coeficiente acumulado después de los términos \(1,\ldots,k\). Para \(1\le d\le9\), cada paso es
\[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}.\]
El arreglo \(a_d\) se actualiza desde el grado \(9\) hasta el \(1\), de modo que cada miembro derecho todavía pertenece a \(A_{k-1}\). Descartar los grados superiores a \(9\) es exacto, no una aproximación: las multiplicaciones posteriores tienen grado no negativo y nunca pueden bajar un coeficiente descartado al grado \(9\).
De forma explícita, el factor de \(A_{k-1}\) es \(d_k(-1)^{k-1}((k-1)!)^2=1/(k(2k-1))\). Inicializamos los coeficientes y la suma como sigue. El término constante de \(A_k\) permanece igual a \(1\), por lo que nunca hay que actualizarlo.
\[a_0^{(0)}=1,\qquad a_d^{(0)}=0\ (1\le d\le9),\qquad s_d^{(0)}=0\ (0\le d\le9).\]
Hay dos requisitos distintos sobre el orden. Primero se acumula la contribución de \(A_{k-1}\) y después se multiplica por el nuevo factor para obtener \(A_k\). Dentro de esa multiplicación, los grados se actualizan en orden descendente. En orden ascendente se usaría un \(a_{d-1}\) ya modificado. Por ejemplo, partiendo de \(1\), aparecería incorrectamente un término \(x^2/k^4\) tras multiplicar solo por \(1-x/k^2\), aunque el producto tiene grado uno.
La justificación formal del truncamiento es que calculamos en el anillo cociente \(\mathbb F_p[x]/(x^{10})\). La reducción módulo \(x^{10}\) es compatible con la suma y el producto. Equivalentemente, el coeficiente de grado \(d\le9\) de un producto solo usa coeficientes de entrada de grados a lo sumo \(d\). Ningún término de grado alto descartado puede afectar a los coeficientes almacenados:
\[[x^d](F(x)G(x))=\sum_{r=0}^{d}([x^r]F(x))([x^{d-r}]G(x)),\qquad 0\le d\le9.\]
Una ejecución secuencial guarda, por tanto, diez coeficientes del producto, pero solo el escalar \(s_9\). Los demás coeficientes de la suma no alimentan a \(s_9\). Los diez coeficientes del producto sí son necesarios: las actualizaciones sucesivas trasladan información desde el término constante, pasando por los grados intermedios, hasta el grado nueve. Además, \(a_9^{(k-1)}=0\) para \(k<10\), así que los primeros nueve sumandos no contribuyen a la parte de interpolación del coeficiente buscado.
Tomemos \(N=3\) y mantengamos fracciones racionales de momento para que el álgebra sea legible. Después podemos reducir estas identidades módulo cualquier primo \(p>6\). Los primeros productos normalizados son
\[A_0=1,\qquad A_1=1-x,\qquad A_2=(1-x)(1-x/4)=1-\frac54x+\frac14x^2.\]
Los pesos para \(k=1,2,3\) son \(1,1/6,1/15\). Cada peso utiliza el producto anterior a la inserción de su propio factor. Al sustituir y agrupar coeficientes obtenemos
\[\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.\]
Esto comprueba tanto los valores prescritos como la fórmula \(d_3=1/60\). También muestra la diferencia entre interpolar y cumplir la condición mónica: con el módulo objetivo, este pequeño interpolante tiene coeficiente principal \(1/60\ne1\). Por tanto, la interpolación por sí sola no completa la construcción. La corrección mónica se deriva más abajo.
Sea \(u_k=k(2k-1)\). Todos los \(u_k\) son no nulos en \(\mathbb F_p\), pues \(1\le k\le N<p\) y \(2k-1<p\). En un lote \(u_1,\ldots,u_m\), se forman los productos prefijo \(U_0=1\), \(U_i=U_{i-1}u_i\). Solo se calcula \(U_m^{-1}=U_m^{p-2}\), y luego se recorre hacia atrás:
\[u_i^{-1}=U_{i-1}U_i^{-1}, \qquad U_{i-1}^{-1}=U_i^{-1}u_i.\]
Una sola exponenciación de Fermat proporciona todos los inversos del lote. Además, \(k^{-1}=(2k-1)u_k^{-1}\), cuyo cuadrado da el \(k^{-2}\) requerido por la actualización del producto. Las implementaciones usan lotes de \(32768\) elementos.
Los índices de la fórmula de inversión por lotes son locales al lote; su primer índice global no tiene por qué ser \(1\). Durante el recorrido hacia atrás mantenemos que la inversa acumulada es \(U_i^{-1}\). Al multiplicarla por el prefijo aún intacto \(U_{i-1}\), obtenemos \(u_i^{-1}\); al multiplicarla por \(u_i\), preparamos la siguiente iteración. Después podemos sobrescribir la posición del prefijo con la inversa obtenida, porque las iteraciones siguientes solo necesitan posiciones anteriores.
Como comprobación concreta de este recurso general, invertimos \(2,3,5\) módulo \(11\). Los productos prefijos son \(1,2,6,8\), y \(8^{-1}=7\). El recorrido inverso da primero \(5^{-1}=6\cdot7=9\), después \(3^{-1}=2\cdot2=4\) y finalmente \(2^{-1}=1\cdot6=6\), todo módulo \(11\). Solo hizo falta una exponenciación. En el algoritmo real, las entradas son los valores no nulos \(k(2k-1)\), y cada multiplicación se reduce módulo \(p\).
Como cada polinomio base de Newton es mónico, el coeficiente principal de \(P_N\) es
\[c_N=d_N= \frac{(-1)^{N-1}N}{(2N-1)(N!)^2}.\]
El polinomio que se anula en todos los nodos es
\[V_N(x)=\prod_{j=0}^{N}(x-j^2) =(-1)^N(N!)^2xA_N(x).\]
La fórmula de \(d_N\) no es nula: su numerador es un signo y ningún factor del denominador se anula. Por tanto, \(P_N\) tiene grado exactamente \(N\), no solo a lo sumo \(N\). Para describir todas las soluciones, sea \(F\) cualquier polinomio con los valores prescritos. Cada factor lineal \(x-j^2\) divide a \(F-P_N\); como las raíces son distintas, también lo divide su producto. Así
\[F(x)=P_N(x)+V_N(x)R(x),\qquad R(x)\in\mathbb F_p[x].\]
Si \(R\ne0\), entonces \(\deg(V_NR)=N+1+\deg R>\deg P_N\), de modo que el término de mayor grado no puede cancelarse. La única solución de grado a lo sumo \(N\) es, por tanto, \(P_N\). Si no es mónica, una solución mónica del grado siguiente exige que \(R\) sea constante, y el coeficiente principal obliga a que esa constante sea exactamente \(1\). Esto demuestra tanto la minimalidad como la unicidad. Dividir simplemente \(P_N\) por \(c_N\) cambiaría el valor exigido en \(1\), de \(1\) a \(c_N^{-1}\), así que escalar no es una corrección válida.
Si \(c_N=1\), el interpolador ya es la única solución mónica de grado mínimo. En caso contrario no puede existir otra solución de grado a lo sumo \(N\); el menor grado mónico posible es \(N+1\), y al sumar el polinomio mónico \(V_N\) se obtiene
\[Q(x)=P_N(x)+V_N(x).\]
Para la instancia objetivo, el código verifica \(c_N\ne1\). Por tanto, el coeficiente pedido se obtiene sin construir ninguno de los dos polinomios grandes:
\[\boxed{[x^{10}]Q=s_9^{(N)} +(-1)^N(N!)^2a_9^{(N)}}.\]
Nota sobre el ejemplo pequeño publicado. Leída literalmente, la expresión «grado mínimo» presenta un caso límite en \((N,M)=(2,7)\): allí \(c_2=1\), por lo que \(P_2(x)=x^2\) ya es mónico y satisface los tres valores. El cúbico mostrado es \(P_2+V_2\), que también los satisface. Este detalle editorial no afecta a la instancia pedida, donde \(c_N\ne1\).
Como ejemplo separado sin el caso especial de monicidad, tomemos \((N,p)=(2,11)\). El interpolante racional \((7x-x^2)/6\) se reduce a \(3x+9x^2\). Al añadir el polinomio que se anula en los nodos obtenemos el siguiente cúbico. Sus valores en \(0,1,4\) son respectivamente \(0,1,2\) módulo \(11\), y su coeficiente principal es \(1\). Esto ilustra cómo la corrección cambia los coeficientes mientras conserva todos los valores exigidos.
\[\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}\]
El programa C++ puede procesar en paralelo intervalos consecutivos de índices. Para un bloque \(B\), escribimos \(g_k(x)=1-x/k^2\) y guardamos, hasta grado \(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.\]
Si el bloque \(L\) precede inmediatamente al bloque \(R\), sus resúmenes se combinan mediante
\[\boxed{G_{LR}=G_LG_R,\qquad H_{LR}=H_L+G_LH_R}.\]
La operación es asociativa, pero el orden importa; la implementación fusiona los bloques en orden creciente de índices. Sus productos factoriales se multiplican al mismo tiempo. Basta la multiplicación polinómica truncada, porque los coeficientes por encima del grado \(9\) no pueden contribuir a grados inferiores.
La regla de bloques también actualiza el estado global de entrada: \((S,A)\mapsto(S+AH_B,AG_B)\). Las contribuciones del bloque derecho necesitan antes todos los factores del izquierdo, lo que explica \(G_LH_R\). Para tres bloques consecutivos, ambas parentizaciones dan exactamente la misma expresión siguiente, demostrando la asociatividad. El bloque vacío tiene resumen \((G,H)=(1,0)\).
\[G_{LMR}=G_LG_MG_R,\qquad H_{LMR}=H_L+G_LH_M+G_LG_MH_R.\]
Esto explica por qué cada trabajador debe conservar los diez coeficientes de \(H_B\), aunque la respuesta final solo requiera el grado nueve. Un producto anterior mezcla los grados mediante la convolución siguiente. El código de combinación puede guardar solo el noveno coeficiente de la suma global, pero no únicamente el noveno de la suma de cada trabajador. Los lotes de inversión y los bloques paralelos tienen funciones distintas: un trabajador puede procesar muchos lotes de inversión dentro de su intervalo consecutivo.
\[[x^9](AH_B)=\sum_{r=0}^{9}([x^r]A)([x^{9-r}]H_B).\]
Lema 1. Los nodos cuadrados son distintos en \(\mathbb F_p\), por lo que \(P_N\) existe y es el único interpolador de grado a lo sumo \(N\).
Lema 2. Evaluar el denominador de la diferencia dividida y la suma binomial alternada produce el \(d_k\) indicado. Al sustituirlo en la fórmula de Newton se obtiene \(P_N/x=\sum A_{k-1}/(k(2k-1))\).
Lema 3. La recurrencia descendente mantiene \(a_d^{(k)}=[x^d]A_k\), mientras que la acumulación mantiene \(s_d^{(k)}=[x^d]\sum_{i=1}^{k}A_{i-1}/(i(2i-1))\). El truncamiento en grado \(9\) conserva exactamente ambos invariantes.
Lema 4. La inversión por lotes devuelve el inverso verdadero de cada \(u_k\ne0\), ya que \(U_{i-1}U_i^{-1}=1/u_i\). Solo cambia la manera de evaluar la misma recurrencia.
Lema 5. Todo polinomio que cumple los datos difiere de \(P_N\) en un múltiplo de \(V_N\). Por ello, la prueba del coeficiente principal elige \(P_N\) exactamente cuando ya es mónico; en caso contrario elige \(P_N+V_N\), mónico y del menor grado posible.
Lema 6. La identidad ordenada de bloques antepone \(G_L\) a toda contribución procedente de \(R\). Una inducción sobre los bloques fusionados reproduce exactamente la suma y el producto seriales.
En conjunto, los lemas demuestran que la expresión final es precisamente el coeficiente de \(x^{10}\) de la inversa funcional mónica de grado mínimo solicitada.
En C++, Block.product[d] representa el coeficiente del producto truncado y Block.sum[d] la suma ponderada correspondiente. La especialización secuencial de calculate_block conserva solo sum[9]; los trabajadores paralelos guardan la suma completa. Python y Java almacenan el mismo producto en coefficients y el noveno coeficiente de la suma en interpolation_coefficient e interpolationCoefficient, respectivamente. El estado matemático es idéntico en las tres implementaciones.
Al comenzar la iteración secuencial \(k\), el arreglo contiene \(A_{k-1}\), el acumulador contiene \(s_9^{(k-1)}\) y la variable factorial contiene \((k-1)!\). Se recupera el peso del lote, se añade su contribución, se actualiza el arreglo en orden descendente y se multiplica la variable factorial por \(k\). Así se restablece el mismo invariante para la iteración \(k+1\). Al terminar, el estado es exactamente \((A_N,s_9^{(N)},N!)\), toda la información que necesita la fórmula final.
C++ usa por defecto como máximo \(16\) hilos trabajadores; el argumento opcional --threads permite elegir otra cantidad admitida. Los trabajadores procesan intervalos consecutivos y se combinan en su orden original. Python y Java evalúan la recurrencia secuencial. La exponenciación modular binaria se realiza mediante power en C++, pow con tres argumentos en Python y modPow en Java. Con este módulo, el producto de dos residuos reducidos es menor que \((10^9+7)^2<2^{63}\), por lo que bastan la aritmética de 64 bits de C++ y long de Java; un producto de 32 bits no basta.
Para validar, cada implementación incluye además un interpolador directo de Lagrange que construye el polinomio completo para \(1\le N\le40\). Usa una construcción distinta de la recurrencia optimizada de Newton. Las comparaciones con esa referencia incluyen el caso ya mónico \(N=1\) y comprueban el coeficiente de grado diez incluso cuando el polinomio pequeño no tiene ese término. Las autopruebas de C++ también comparan resultados secuenciales y paralelos cerca de los límites de los lotes de inversión; Python y Java prueban asimismo un tamaño de lote deliberadamente pequeño.
Sea \(D=9\), sea \(B=32768\) el tamaño de los lotes de inversión y sea \(T\) el número de trabajadores, con \(T=1\) en la ejecución secuencial. Cada índice requiere \(D\) actualizaciones de coeficientes del producto y una cantidad constante de aritmética adicional; una suma completa de trabajador añade \(D+1\) actualizaciones. Los recorridos de prefijos y hacia atrás son lineales en el tamaño del lote. El número total de inversiones por lotes no supera \(\lceil N/B\rceil+T-1\), y la prueba del coeficiente principal usa una inversión adicional.
\[W(N,D,B,T,p)=O\!\left(ND+\left(\frac NB+T\right)\log p+TD^2\right).\]
El término \(TD^2\) corresponde a las convoluciones cortas de polinomios durante las combinaciones. Con \(D,B,T\) y este módulo fijos, el trabajo total es \(O(N)\). Para el \(N\) objetivo, una ejecución secuencial realiza \(90{,}000{,}000\) actualizaciones de coeficientes del producto y solo \(306\) inversiones por lotes, más la inversión final del coeficiente principal. El paralelismo reparte el bucle principal entre trabajadores sin cambiar el orden del trabajo total; el tiempo real también depende de la planificación y del hardware.
El almacenamiento auxiliar es \(O(B+D)\) en la implementación secuencial y \(O(T(B+D))\) en la paralela, sin contar las pilas de los hilos del entorno de ejecución. No hay un arreglo para los \(N\) nodos ni un arreglo de coeficientes hasta grado \(N\). La misma derivación permite calcular otro grado positivo \(r\) conservando los grados de \(0\) a \(r-1\); que aquí se pida el grado diez limita el estado a diez coeficientes.
对于素数模数 \(M\)，\(x^2\) 的 \((N,M)\)-函数逆是次数最低的首一多项式 \(Q(x)\)。其系数以 \(\{0,\ldots,M-1\}\) 中的数表示，并满足
\[Q(n^2)\equiv n\pmod M\qquad(0\le n\le N).\]
本题中 \(N=10^7\)、\(M=p=10^9+7\)，只需要 \(x^{10}\) 的系数。无须构造次数为 \(N\) 的插值多项式，甚至无须保存它的全部系数；计算可以严格而准确地截断到 \(9\) 次。
记号 \([x^r]F(x)\) 表示 \(F\) 中 \(x^r\) 的系数。除非明确说明是有理数恒等式，下文的等式都在 \(\mathbb F_p\) 中成立。分数 \(a/b\) 表示 \(a\,b^{-1}\pmod p\)，不是整数除法或浮点除法。由于 \(p\) 是素数，每个非零分母都有唯一的逆元。将结果取为 \(0\) 到 \(p-1\) 之间的代表元，就得到题目要求的系数形式。
题目并不要求在整个有限域上逆转平方运算，因为 \(n\) 与 \(-n\) 的平方相同。它只在选定的平方节点上各指定一个值。条件 \(2N<p\) 保证两个选定的正下标不会成为同一节点的互为相反数的平方根。因此，应从插值入手，而不是使用实数平方根函数的幂级数展开。
在有限域 \(\mathbb F_p\) 上计算。节点 \(0^2,1^2,\ldots,N^2\) 两两不同：若 \(i^2\equiv j^2\pmod p\)，则 \(p\mid(i-j)(i+j)\)；然而 \(0\le i,j\le N\) 且 \(i+j\le2N<p\)。因此恰有一个次数不超过 \(N\) 的多项式 \(P_N\)，在所有指定节点上满足 \(P_N(n^2)=n\)。
把它写成 Newton 形式：
\[P_N(x)=\sum_{k=1}^{N}d_k\prod_{j=0}^{k-1}(x-j^2).\]
由于 \(0\) 点的指定值也是 \(0\)，所以 \(k=0\) 项为零。
存在性和唯一性可以直接证明。固定 \(k\)，定义下面的拉格朗日基。在 \(x=j^2\) 处，第 \(i\) 个基多项式在 \(i=j\) 时等于 \(1\)，否则等于 \(0\)，因此它们的加权和满足所有指定值。若两个次数至多为 \(k\) 的多项式都满足这些值，则其差有 \(k+1\) 个不同的根；次数至多为 \(k\) 的非零多项式不可能如此。
\[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).\]
牛顿形式逐个加入节点。差 \(P_k-P_{k-1}\) 在前 \(k\) 个节点上为零，所以它是这些节点对应的首一乘积 \(B_k\) 的常数倍。在新节点处求值便可确定该常数。由于 \(P_{k-1}\) 的次数至多为 \(k-1\)，\(d_k\) 也等于 \(P_k\) 的拉格朗日表达式中 \(x^k\) 的系数。读出这个首项系数，就得到下一节的差商求和公式。
\[\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}\]
节点 \(0^2,\ldots,k^2\) 上的差商为
\[d_k=\sum_{i=1}^{k}\frac{i} {\displaystyle\prod_{\substack{0\le j\le k\\j\ne i}}(i^2-j^2)}.\]
计算分母时，将两个线性因子乘积分开处理。第一个乘积恰有 \(k-i\) 个负因子。第二个乘积对应从 \(i\) 到 \(i+k\) 的连续整数乘积，但必须除去遗漏的因子 \(2i\)。因此，对于 \(1\le i\le k\)，有
\[\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}\]
将两式相乘，\(i!/(i-1)!=i\) 与分母中的 \(2i\) 约去。所有阶乘的参数都至多为 \(2k\le2N<p\)，所以这些阶乘模 \(p\) 都可逆。于是
\[\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)!}.\]
下面逐步化简剩余的和，不略去抵消过程。两边乘以 \((2k)!\)，令 \(r=k-i\)，并将结果记为 \(S_k\)。阶乘分母随之变为二项式系数：
\[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).\]
对交错和逐项使用帕斯卡恒等式，再将第二个和的下标平移一位，所有内部项都会抵消，只留下末项。约定下标为 \(-1\) 的二项式系数以及空和都等于零。
\[\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}\]
将二项式系数之差代入 \(S_k\)。下标平移时交错符号也改变，因此差式中的减号变成加号。两次应用部分和恒等式，再利用相邻二项式系数的比值，即可得到下面的化简。采用空和约定后，\(k=1\) 时也成立：
\[\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)!}}.\]
定义
\[A_k(x)=\prod_{j=1}^{k}\left(1-\frac{x}{j^2}\right), \qquad A_0(x)=1.\]
由于
\[\prod_{j=0}^{k-1}(x-j^2) =x(-1)^{k-1}((k-1)!)^2A_{k-1}(x),\]
符号和阶乘与 \(d_k\) 抵消，得到简洁恒等式
\[\boxed{\frac{P_N(x)}x= \sum_{k=1}^{N}\frac{A_{k-1}(x)}{k(2k-1)}}.\]
所以 \([x^{10}]P_N=[x^9](P_N/x)\)。令 \(a_d^{(k)}=[x^d]A_k\)，并以 \(s_d^{(k)}\) 表示累加完第 \(1,\ldots,k\) 项后的系数。对 \(1\le d\le9\)，每一步为
\[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}.\]
\(a_d\) 数组按次数从 \(9\) 到 \(1\) 逆序更新，因此右端始终还是 \(A_{k-1}\) 的系数。丢弃 \(9\) 次以上的项是精确操作而非近似：后续乘法只含非负次数，不可能把已丢弃的高次系数降回 \(9\) 次。
具体地说，\(A_{k-1}\) 前的标量为 \(d_k(-1)^{k-1}((k-1)!)^2=1/(k(2k-1))\)。系数状态与累加和按下式初始化。\(A_k\) 的常数项始终为 \(1\)，因此无需更新。
\[a_0^{(0)}=1,\qquad a_d^{(0)}=0\ (1\le d\le9),\qquad s_d^{(0)}=0\ (0\le d\le9).\]
这里有两个独立的顺序要求：先累加 \(A_{k-1}\) 的贡献，再乘入新因子得到 \(A_k\)；而乘法内部必须按次数递减更新。若按递增顺序更新，就会用到已经改变的 \(a_{d-1}\)。例如，从 \(1\) 出发，仅乘以 \(1-x/k^2\) 就会错误地产生 \(x^2/k^4\) 项，而正确结果只有一次。
截断的严格依据是我们在商环 \(\mathbb F_p[x]/(x^{10})\) 中计算。模 \(x^{10}\) 的约化与加法、乘法相容。等价地，乘积中 \(d\le9\) 次项的系数仅依赖输入中次数不超过 \(d\) 的系数。因此，任何被舍弃的高次项都不会影响保留的系数：
\[[x^d](F(x)G(x))=\sum_{r=0}^{d}([x^r]F(x))([x^{d-r}]G(x)),\qquad 0\le d\le9.\]
因此，串行计算保存乘积的十个系数，但累加和只需保存标量 \(s_9\)，因为其他累加系数不会传入 \(s_9\)。不过，乘积的十个系数都不可少：连续更新会把常数项的信息逐级传递到九次项。另外，当 \(k<10\) 时，\(a_9^{(k-1)}=0\)，所以前九个求和项对目标系数的插值部分没有贡献。
取 \(N=3\)，暂时保留有理数分数，以便看清代数过程。之后可将这些恒等式按任意素数 \(p>6\) 取模。最初的几个归一化乘积为
\[A_0=1,\qquad A_1=1-x,\qquad A_2=(1-x)(1-x/4)=1-\frac54x+\frac14x^2.\]
\(k=1,2,3\) 对应的权重为 \(1,1/6,1/15\)。每个权重使用的都是尚未乘入自身因子之前的乘积。代入并合并同类项，得到
\[\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.\]
这同时验证了指定值与系数公式 \(d_3=1/60\)。它也说明插值条件和首一条件是两回事：在目标模数下，这个小插值多项式的首项系数为 \(1/60\ne1\)，所以仅完成插值还不够。下面将推导首一修正。
令 \(u_k=k(2k-1)\)。因为 \(1\le k\le N<p\) 且 \(2k-1<p\)，所有 \(u_k\) 在 \(\mathbb F_p\) 中都非零。对一批 \(u_1,\ldots,u_m\)，构造前缀积 \(U_0=1\)、\(U_i=U_{i-1}u_i\)。只计算一次 \(U_m^{-1}=U_m^{p-2}\)，然后反向扫描：
\[u_i^{-1}=U_{i-1}U_i^{-1}, \qquad U_{i-1}^{-1}=U_i^{-1}u_i.\]
这样，一次 Fermat 幂运算就得到整批元素的逆元。此外 \(k^{-1}=(2k-1)u_k^{-1}\)，平方后即得乘积更新需要的 \(k^{-2}\)。实现采用每批 \(32768\) 个元素。
批量求逆公式中的下标是批内下标，该批的起始全局下标不必为 \(1\)。反向遍历时，保持当前逆元等于 \(U_i^{-1}\) 这一不变式。将它乘以尚未修改的前缀 \(U_{i-1}\)，便得到 \(u_i^{-1}\)；再将当前逆元乘以 \(u_i\)，即可进入下一次迭代。随后可以用求得的逆元覆盖该前缀位置，因为后续迭代只会访问更早的位置。
用一个具体例子检查这个通用技巧：在模 \(11\) 下求 \(2,3,5\) 的逆元。前缀乘积为 \(1,2,6,8\)，且 \(8^{-1}=7\)。反向遍历先得到 \(5^{-1}=6\cdot7=9\)，再得到 \(3^{-1}=2\cdot2=4\)，最后得到 \(2^{-1}=1\cdot6=6\)，所有等式均按模 \(11\) 理解。整个过程只需一次幂运算。实际算法的输入是非零的 \(k(2k-1)\)，每次乘法都模 \(p\) 约化。
每个 Newton 基多项式都是首一的，因此 \(P_N\) 的首项系数为
\[c_N=d_N= \frac{(-1)^{N-1}N}{(2N-1)(N!)^2}.\]
在所有插值节点上为零的多项式是
\[V_N(x)=\prod_{j=0}^{N}(x-j^2) =(-1)^N(N!)^2xA_N(x).\]
\(d_N\) 的公式不为零：分子只是正负号，分母中的因子也都非零。因此，\(P_N\) 的次数恰好是 \(N\)，而不只是至多为 \(N\)。为了描述所有解，设 \(F\) 为任意满足指定值的多项式。每个线性因子 \(x-j^2\) 都整除 \(F-P_N\)；由于这些根互不相同，它们的乘积也整除该差。因此
\[F(x)=P_N(x)+V_N(x)R(x),\qquad R(x)\in\mathbb F_p[x].\]
若 \(R\ne0\)，则 \(\deg(V_NR)=N+1+\deg R>\deg P_N\)，最高次项不可能抵消。所以次数至多为 \(N\) 的解只有 \(P_N\)。如果它不是首一的，那么要在下一次数得到首一解，\(R\) 必须是常数，且首项系数条件迫使这个常数恰好等于 \(1\)。这同时证明了最小次数和唯一性。直接将 \(P_N\) 除以 \(c_N\) 会使 \(1\) 处要求的值从 \(1\) 变成 \(c_N^{-1}\)，因此缩放不能作为合法修正。
若 \(c_N=1\)，插值多项式本身已经是唯一的最低次数首一解。否则，次数不超过 \(N\) 的其他解不存在，最低可能的首一次数为 \(N+1\)；加上首一多项式 \(V_N\) 即得
\[Q(x)=P_N(x)+V_N(x).\]
代码对目标实例验证了 \(c_N\ne1\)。因此不用构造任一巨大多项式便可得到所求系数：
\[\boxed{[x^{10}]Q=s_9^{(N)} +(-1)^N(N!)^2a_9^{(N)}}.\]
关于题面小例子的说明。 若严格按字面理解“最低次数”，\((N,M)=(2,7)\) 是一个特殊边界情形：此时 \(c_2=1\)，所以 \(P_2(x)=x^2\) 已经首一并满足三个给定值。题面展示的三次式是 \(P_2+V_2\)，也满足这些值。这一编辑层面的边界情形不影响目标实例，因为目标中 \(c_N\ne1\)。
再看一个没有首一特殊情形的例子，取 \((N,p)=(2,11)\)。有理数插值多项式 \((7x-x^2)/6\) 模约化后为 \(3x+9x^2\)。加上节点消失多项式后，得到下面的三次多项式；它在 \(0,1,4\) 处的值模 \(11\) 分别为 \(0,1,2\)，且首项系数为 \(1\)。这说明修正虽然改变系数，却保留了全部指定值。
\[\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}\]
C++ 程序可以并行处理连续的下标区间。对区块 \(B\)，记 \(g_k(x)=1-x/k^2\)，并把下列多项式保留到 \(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.\]
若区块 \(L\) 紧接在区块 \(R\) 之前，则摘要按下式合并：
\[\boxed{G_{LR}=G_LG_R,\qquad H_{LR}=H_L+G_LH_R}.\]
该运算满足结合律，但次序不可交换；实现按下标递增顺序合并区块，同时相乘各区块的阶乘积。截断的多项式乘法已经足够，因为 \(9\) 次以上的系数无法影响更低次数。
分块规则也可以看作对输入全局状态的更新：\((S,A)\mapsto(S+AH_B,AG_B)\)。右块的每项贡献之前都必须乘上左块的全部因子，因此出现 \(G_LH_R\)。对三个连续块，无论怎样加括号都会得到下面同一表达式，从而证明结合律。空块的摘要为 \((G,H)=(1,0)\)。
\[G_{LMR}=G_LG_MG_R,\qquad H_{LMR}=H_L+G_LH_M+G_LG_MH_R.\]
这也解释了为什么每个工作线程必须保存 \(H_B\) 的十个系数，尽管最终只需要九次项。先前的乘积会按下面的卷积混合各次数。合并代码可以只保存全局累加和的九次项，却不能只保存各线程局部累加和的九次项。求逆批次与并行块用途不同：一个线程可以在自己的连续区间内处理多个求逆批次。
\[[x^9](AH_B)=\sum_{r=0}^{9}([x^r]A)([x^{9-r}]H_B).\]
引理 1。平方节点在 \(\mathbb F_p\) 中互不相同，所以 \(P_N\) 存在，并且是次数不超过 \(N\) 的唯一插值多项式。
引理 2。计算差商分母和交错二项式和可得上述 \(d_k\)。代入 Newton 公式便有 \(P_N/x=\sum A_{k-1}/(k(2k-1))\)。
引理 3。逆序递推保持 \(a_d^{(k)}=[x^d]A_k\)，累加则保持 \(s_d^{(k)}=[x^d]\sum_{i=1}^{k}A_{i-1}/(i(2i-1))\)。截断到 \(9\) 次完整保留这两个不变量。
引理 4。由于 \(U_{i-1}U_i^{-1}=1/u_i\)，批量求逆确实返回每个非零 \(u_k\) 的逆元；它只改变同一递推的求值方式。
引理 5。任何满足给定数据的多项式与 \(P_N\) 之差都是 \(V_N\) 的倍数。因此首项系数检查在 \(P_N\) 已首一时选择它，否则选择次数尽可能低的首一多项式 \(P_N+V_N\)。
引理 6。有序区块恒等式把 \(G_L\) 放在来自 \(R\) 的每个贡献之前。对已合并区块作归纳，即可精确恢复串行的和与积。
综合这些引理，最终表达式正是所求最低次数首一函数逆中 \(x^{10}\) 的系数。
在 C++ 中，Block.product[d] 表示截断乘积的系数，Block.sum[d] 表示对应的加权和。calculate_block 的串行特化只保留 sum[9]，而并行工作线程保留完整的和。Python 和 Java 都将乘积存入 coefficients，将累加和的九次项系数分别存入 interpolation_coefficient 和 interpolationCoefficient。三种实现维护的是同一数学状态。
串行第 \(k\) 次迭代开始时，数组保存 \(A_{k-1}\)，累加器保存 \(s_9^{(k-1)}\)，阶乘变量保存 \((k-1)!\)。取出批次中的权重并累加其贡献，按次数递减更新数组，再将阶乘变量乘以 \(k\)。这些操作为第 \(k+1\) 次迭代重新建立同一不变式。最后一次迭代结束后，状态恰好为 \((A_N,s_9^{(N)},N!)\)，其中包含最终公式所需的全部信息。
C++ 默认最多使用 \(16\) 个工作线程，可选参数 --threads 允许选择其他受支持的数量。线程处理连续区间，并按原顺序合并。Python 和 Java 则执行串行递推。二进制模幂运算分别由 C++ 的 power、Python 的三参数 pow 和 Java 的 modPow 完成。对本题模数，两个已约化剩余类代表元的乘积小于 \((10^9+7)^2<2^{63}\)，因此 C++ 的 64 位运算和 Java 的 long 足够，但 32 位乘法不够。
为进行验证，每种实现还包含直接拉格朗日插值，在 \(1\le N\le40\) 时构造完整多项式，其构造方式不同于优化后的牛顿递推。与该参考方法的比较包含 \(N=1\) 时已经首一的情况，也会检查小多项式没有十次项时的对应系数。C++ 自测还会在求逆批次边界附近比较串行与并行结果；Python 和 Java 也测试了刻意设小的批次大小。
令 \(D=9\)，求逆批次大小 \(B=32768\)，工作线程数为 \(T\)，串行时 \(T=1\)。每个下标需要更新 \(D\) 个乘积系数，并进行常数次其他运算；保存完整线程累加和时，还要增加 \(D+1\) 次更新。前缀与反向遍历的工作量均与批次大小成正比。批次求逆总次数至多为 \(\lceil N/B\rceil+T-1\)，首项系数判断另需一次求逆。
\[W(N,D,B,T,p)=O\!\left(ND+\left(\frac NB+T\right)\log p+TD^2\right).\]
\(TD^2\) 项对应合并时的短多项式卷积。在 \(D,B,T\) 及本题模数固定时，总工作量为 \(O(N)\)。对目标 \(N\)，串行计算进行 \(90{,}000{,}000\) 次乘积系数更新，仅需 \(306\) 次批次求逆，外加最后一次首项系数求逆。并行化把主循环分配给多个线程，但不改变总工作量的阶；实际耗时还取决于调度和硬件。
不计运行时线程栈，串行实现的辅助空间为 \(O(B+D)\)，并行实现为 \(O(T(B+D))\)。算法不建立覆盖全部 \(N\) 个节点的数组，也不建立次数达到 \(N\) 的系数数组。同样的推导还能计算任意其他正次数 \(r\) 的系数，只需保留 \(0\) 到 \(r-1\) 次项；本题要求十次项，因此状态中只需十个系数。
Для простого модуля \(M\) \((N,M)\)-функциональной обратной к \(x^2\) называется унитарный многочлен \(Q(x)\) наименьшей степени. Его коэффициенты представлены числами из \(\{0,\ldots,M-1\}\), и он удовлетворяет условиям
\[Q(n^2)\equiv n\pmod M\qquad(0\le n\le N).\]
Здесь \(N=10^7\), \(M=p=10^9+7\), а требуется лишь коэффициент при \(x^{10}\). Строить интерполяционный многочлен степени \(N\) или хранить все его коэффициенты не нужно: вычисление можно строго и без потери точности усечь на степени \(9\).
Обозначение \([x^r]F(x)\) означает коэффициент при \(x^r\) в \(F\). Все равенства ниже рассматриваются в \(\mathbb F_p\), если явно не указано, что это рациональные тождества. Дробь \(a/b\) означает \(a\,b^{-1}\pmod p\), а не целочисленное деление или вычисление с плавающей точкой. Поскольку \(p\) простое, каждый ненулевой знаменатель имеет единственный обратный элемент. Выбор представителя от \(0\) до \(p-1\) даёт требуемый формат коэффициентов.
Задача не требует обратить возведение в квадрат на всём поле: \(n\) и \(-n\) имеют одинаковые квадраты. Значение задано лишь в каждом из выбранных квадратных узлов. Условие \(2N<p\) не позволяет двум выбранным положительным индексам оказаться противоположными квадратными корнями одного узла. Поэтому подходящей отправной точкой служит интерполяция, а не разложение вещественного квадратного корня в степенной ряд.
Будем работать в поле \(\mathbb F_p\). Узлы \(0^2,1^2,\ldots,N^2\) попарно различны: из \(i^2\equiv j^2\pmod p\) следовало бы \(p\mid(i-j)(i+j)\), но \(0\le i,j\le N\) и \(i+j\le2N<p\). Поэтому существует единственный многочлен \(P_N\) степени не выше \(N\), для которого \(P_N(n^2)=n\) во всех требуемых узлах.
Запишем его в форме Ньютона:
\[P_N(x)=\sum_{k=1}^{N}d_k\prod_{j=0}^{k-1}(x-j^2).\]
Член \(k=0\) равен нулю, поскольку заданное значение в нуле также равно нулю.
Дадим явное доказательство существования и единственности. При фиксированном \(k\) определим следующий базис Лагранжа. В точке \(x=j^2\) его \(i\)-й элемент равен \(1\) при \(i=j\) и \(0\) иначе. Поэтому взвешенная сумма принимает заданные значения. Если бы два многочлена степени не выше \(k\) имели эти значения, их разность обладала бы \(k+1\) различными корнями, что невозможно для ненулевого многочлена такой степени.
\[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).\]
Форма Ньютона добавляет узлы по одному. Разность \(P_k-P_{k-1}\) обращается в нуль в первых \(k\) узлах, поэтому она является постоянным кратным их приведённого произведения \(B_k\). Подстановка нового узла определяет эту константу. Поскольку степень \(P_{k-1}\) не выше \(k-1\), число \(d_k\) также является коэффициентом при \(x^k\) в формуле Лагранжа для \(P_k\). Выделяя этот старший коэффициент, получаем следующую сумму разделённых разностей.
\[\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}\]
Разделённая разность на узлах \(0^2,\ldots,k^2\) равна
\[d_k=\sum_{i=1}^{k}\frac{i} {\displaystyle\prod_{\substack{0\le j\le k\\j\ne i}}(i^2-j^2)}.\]
Чтобы вычислить знаменатель, рассмотрим два произведения линейных множителей отдельно. В первом ровно \(k-i\) отрицательных множителей. Во втором из произведения последовательных чисел от \(i\) до \(i+k\) нужно убрать пропущенный множитель \(2i\). При \(1\le i\le k\) получаем
\[\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}\]
При умножении выражений \(i!/(i-1)!=i\) сокращается с множителем \(2i\). Все встречающиеся факториалы обратимы по модулю \(p\), поскольку их аргументы не превосходят \(2k\le2N<p\). Следовательно
\[\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)!}.\]
Оставшуюся сумму можно вычислить, не пропуская сокращений. Умножим на \((2k)!\), положим \(r=k-i\) и обозначим результат через \(S_k\). Факториальный знаменатель превращается в биномиальный коэффициент:
\[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).\]
Для знакопеременных сумм применим тождество Паскаля к каждому слагаемому. После сдвига индекса второй суммы на единицу все внутренние члены сокращаются и остаётся только последний. Биномиальный коэффициент с нижним индексом \(-1\) и пустую сумму считаем равными нулю.
\[\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}\]
Подставим разность биномиальных коэффициентов в \(S_k\). После сдвига индекса минус перед разностью становится плюсом, поскольку меняется и знак слагаемого. Дважды применим тождество для частичных сумм, затем воспользуемся отношением соседних биномиальных коэффициентов. При принятом соглашении о пустой сумме вычисление верно и для \(k=1\):
\[\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)!}}.\]
Определим
\[A_k(x)=\prod_{j=1}^{k}\left(1-\frac{x}{j^2}\right), \qquad A_0(x)=1.\]
Так как
\[\prod_{j=0}^{k-1}(x-j^2) =x(-1)^{k-1}((k-1)!)^2A_{k-1}(x),\]
знаки и факториалы сокращаются с \(d_k\), и остаётся простое тождество
\[\boxed{\frac{P_N(x)}x= \sum_{k=1}^{N}\frac{A_{k-1}(x)}{k(2k-1)}}.\]
Значит, \([x^{10}]P_N=[x^9](P_N/x)\). Обозначим \(a_d^{(k)}=[x^d]A_k\), а через \(s_d^{(k)}\) — накопленный коэффициент после членов \(1,\ldots,k\). Для \(1\le d\le9\) один шаг имеет вид
\[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}.\]
Массив \(a_d\) обновляется по убыванию степени от \(9\) до \(1\), поэтому справа остаются коэффициенты \(A_{k-1}\). Отбрасывание степеней выше \(9\) является точным, а не приближённым: последующие множители имеют только неотрицательные степени и не могут вернуть отброшенный коэффициент к степени \(9\).
В явном виде множитель перед \(A_{k-1}\) равен \(d_k(-1)^{k-1}((k-1)!)^2=1/(k(2k-1))\). Начальные значения коэффициентов и суммы задаются ниже. Свободный член \(A_k\) всегда равен \(1\), поэтому обновлять его не требуется.
\[a_0^{(0)}=1,\qquad a_d^{(0)}=0\ (1\le d\le9),\qquad s_d^{(0)}=0\ (0\le d\le9).\]
Важны два независимых требования к порядку действий. Сначала прибавляем вклад \(A_{k-1}\), затем умножаем на новый множитель и получаем \(A_k\). При этом коэффициенты обновляются по убыванию степени. При возрастании использовалось бы уже изменённое \(a_{d-1}\). Например, начиная с \(1\), мы ошибочно получили бы член \(x^2/k^4\) после умножения только на \(1-x/k^2\), хотя произведение имеет степень один.
Формальное обоснование усечения состоит в вычислении в факторкольце \(\mathbb F_p[x]/(x^{10})\). Приведение по модулю \(x^{10}\) согласуется со сложением и умножением. Иначе говоря, коэффициент степени \(d\le9\) в произведении зависит только от входных коэффициентов степеней не выше \(d\). Ни один отброшенный член высокой степени не влияет на сохранённые коэффициенты:
\[[x^d](F(x)G(x))=\sum_{r=0}^{d}([x^r]F(x))([x^{d-r}]G(x)),\qquad 0\le d\le9.\]
Поэтому последовательное вычисление хранит десять коэффициентов произведения, но лишь скаляр \(s_9\) для суммы. Остальные коэффициенты суммы не участвуют в обновлении \(s_9\). Однако все десять коэффициентов произведения необходимы: последовательные обновления передают информацию от свободного члена через промежуточные степени к девятой. Кроме того, при \(k<10\) имеем \(a_9^{(k-1)}=0\), поэтому первые девять слагаемых не дают вклада в интерполяционную часть искомого коэффициента.
Возьмём \(N=3\) и временно сохраним рациональные дроби, чтобы алгебраические действия были наглядными. Затем эти тождества можно привести по модулю любого простого \(p>6\). Первые нормированные произведения имеют вид
\[A_0=1,\qquad A_1=1-x,\qquad A_2=(1-x)(1-x/4)=1-\frac54x+\frac14x^2.\]
Веса для \(k=1,2,3\) равны \(1,1/6,1/15\). Каждый вес умножается на произведение до включения соответствующего ему множителя. Подстановка и сбор коэффициентов дают
\[\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.\]
Так проверяются и заданные значения, и формула \(d_3=1/60\). Пример также показывает различие между интерполяцией и условием приведённости: при целевом модуле старший коэффициент этого малого интерполянта равен \(1/60\ne1\). Поэтому одной интерполяции недостаточно для полной конструкции. Поправка, обеспечивающая приведённость, выводится ниже.
Пусть \(u_k=k(2k-1)\). Все \(u_k\) ненулевые в \(\mathbb F_p\), поскольку \(1\le k\le N<p\) и \(2k-1<p\). Для пакета \(u_1,\ldots,u_m\) строятся префиксные произведения \(U_0=1\), \(U_i=U_{i-1}u_i\). Вычисляется лишь \(U_m^{-1}=U_m^{p-2}\), после чего выполняется обратный проход:
\[u_i^{-1}=U_{i-1}U_i^{-1}, \qquad U_{i-1}^{-1}=U_i^{-1}u_i.\]
Одно возведение в степень по теореме Ферма даёт все обратные элементы пакета. Кроме того, \(k^{-1}=(2k-1)u_k^{-1}\), а его квадрат есть \(k^{-2}\), нужный при обновлении произведения. В реализациях размер пакета равен \(32768\).
Индексы в формуле пакетного обращения локальны для пакета; его первый глобальный индекс не обязан быть \(1\). При обратном проходе поддерживается инвариант: текущий обратный элемент равен \(U_i^{-1}\). Умножение на ещё не изменённый префикс \(U_{i-1}\) даёт \(u_i^{-1}\), а умножение текущего обратного элемента на \(u_i\) подготавливает следующую итерацию. Затем ячейку префикса можно заменить найденным обратным элементом, поскольку далее потребуются только более ранние ячейки.
Для конкретной проверки этого общего приёма найдём обратные к \(2,3,5\) по модулю \(11\). Префиксные произведения равны \(1,2,6,8\), а \(8^{-1}=7\). Обратный проход сначала даёт \(5^{-1}=6\cdot7=9\), затем \(3^{-1}=2\cdot2=4\), наконец \(2^{-1}=1\cdot6=6\), всё по модулю \(11\). Потребовалось лишь одно возведение в степень. В основном алгоритме входами служат ненулевые значения \(k(2k-1)\), а каждое произведение приводится по модулю \(p\).
Каждый базисный многочлен Ньютона унитарен, поэтому старший коэффициент \(P_N\) равен
\[c_N=d_N= \frac{(-1)^{N-1}N}{(2N-1)(N!)^2}.\]
Многочлен, обращающийся в ноль во всех узлах интерполяции, имеет вид
\[V_N(x)=\prod_{j=0}^{N}(x-j^2) =(-1)^N(N!)^2xA_N(x).\]
Формула для \(d_N\) даёт ненулевой элемент: числитель равен знаку, а ни один множитель знаменателя не обращается в нуль. Значит, степень \(P_N\) равна в точности \(N\), а не просто не превосходит \(N\). Для описания всех решений пусть \(F\) — любой многочлен с заданными значениями. Каждый линейный множитель \(x-j^2\) делит \(F-P_N\); поскольку корни различны, их произведение также делит эту разность. Поэтому
\[F(x)=P_N(x)+V_N(x)R(x),\qquad R(x)\in\mathbb F_p[x].\]
Если \(R\ne0\), то \(\deg(V_NR)=N+1+\deg R>\deg P_N\), поэтому старший член не может сократиться. Следовательно, единственное решение степени не выше \(N\) — это \(P_N\). Если оно не приведённое, то для приведённого решения следующей степени \(R\) обязано быть константой, а условие на старший коэффициент заставляет эту константу равняться \(1\). Это доказывает и минимальность, и единственность. Простое деление \(P_N\) на \(c_N\) изменило бы требуемое значение в точке \(1\) с \(1\) на \(c_N^{-1}\), поэтому масштабирование не является допустимой поправкой.
Если \(c_N=1\), интерполянт уже является единственным унитарным решением наименьшей степени. Иначе другого решения степени не выше \(N\) нет, наименьшая возможная степень унитарного решения равна \(N+1\), и добавление унитарного \(V_N\) даёт
\[Q(x)=P_N(x)+V_N(x).\]
Для целевого случая программа проверяет \(c_N\ne1\). Поэтому искомый коэффициент находится без построения обоих больших многочленов:
\[\boxed{[x^{10}]Q=s_9^{(N)} +(-1)^N(N!)^2a_9^{(N)}}.\]
Замечание об опубликованном малом примере. При буквальном чтении слов «наименьшей степени» случай \((N,M)=(2,7)\) оказывается особым: здесь \(c_2=1\), поэтому \(P_2(x)=x^2\) уже унитарен и удовлетворяет всем трём значениям. Приведённый кубический многочлен равен \(P_2+V_2\) и также им удовлетворяет. Эта редакционная особенность не влияет на целевой случай, где \(c_N\ne1\).
В качестве отдельного примера без особого случая приведённости возьмём \((N,p)=(2,11)\). Рациональный интерполянт \((7x-x^2)/6\) приводится к \(3x+9x^2\). Добавление многочлена, исчезающего в узлах, даёт следующий кубический многочлен. Его значения в \(0,1,4\) равны соответственно \(0,1,2\) по модулю \(11\), а старший коэффициент равен \(1\). Это показывает, как поправка меняет коэффициенты, сохраняя все заданные значения.
\[\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}\]
Программа на C++ может параллельно обрабатывать последовательные интервалы индексов. Для блока \(B\) положим \(g_k(x)=1-x/k^2\) и сохраним до степени \(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.\]
Если блок \(L\) непосредственно предшествует блоку \(R\), их сводки объединяются по формулам
\[\boxed{G_{LR}=G_LG_R,\qquad H_{LR}=H_L+G_LH_R}.\]
Операция ассоциативна, но порядок существенен; реализация объединяет блоки по возрастанию индексов. Одновременно перемножаются их факториальные произведения. Достаточно усеченного умножения многочленов, поскольку коэффициенты степеней выше \(9\) не влияют на меньшие степени.
Правило объединения блоков можно понимать и как обновление входного глобального состояния: \((S,A)\mapsto(S+AH_B,AG_B)\). Перед вкладом правого блока должны стоять все множители левого, откуда возникает \(G_LH_R\). Для трёх последовательных блоков обе расстановки скобок дают одно и то же выражение ниже, что доказывает ассоциативность. Пустому блоку соответствует сводка \((G,H)=(1,0)\).
\[G_{LMR}=G_LG_MG_R,\qquad H_{LMR}=H_L+G_LH_M+G_LG_MH_R.\]
Это объясняет, почему рабочий поток должен хранить все десять коэффициентов \(H_B\), хотя в ответе нужна только девятая степень. Предыдущее произведение смешивает степени согласно свёртке ниже. При объединении можно хранить лишь девятый коэффициент глобальной суммы, но не только девятый коэффициент суммы каждого потока. Пакеты обращения и параллельные блоки решают разные задачи: один поток может обработать много пакетов обращения внутри своего последовательного интервала.
\[[x^9](AH_B)=\sum_{r=0}^{9}([x^r]A)([x^{9-r}]H_B).\]
Лемма 1. Квадратные узлы различны в \(\mathbb F_p\), поэтому \(P_N\) существует и является единственным интерполянтом степени не выше \(N\).
Лемма 2. Вычисление знаменателя разделённой разности и знакопеременной биномиальной суммы даёт указанное \(d_k\). Подстановка в формулу Ньютона даёт \(P_N/x=\sum A_{k-1}/(k(2k-1))\).
Лемма 3. Рекурсия по убыванию степени сохраняет \(a_d^{(k)}=[x^d]A_k\), а накопление сохраняет \(s_d^{(k)}=[x^d]\sum_{i=1}^{k}A_{i-1}/(i(2i-1))\). Усечение на степени \(9\) сохраняет обе инварианты точно.
Лемма 4. Пакетное обращение возвращает истинный обратный элемент каждого \(u_k\ne0\), так как \(U_{i-1}U_i^{-1}=1/u_i\). Оно меняет лишь способ вычисления той же рекурсии.
Лемма 5. Любой многочлен с заданными значениями отличается от \(P_N\) на кратное \(V_N\). Поэтому проверка старшего коэффициента выбирает \(P_N\) именно тогда, когда он уже унитарен, а иначе — унитарный многочлен \(P_N+V_N\) наименьшей возможной степени.
Лемма 6. Упорядоченное блочное тождество ставит \(G_L\) перед каждым вкладом из \(R\). Индукция по объединяемым блокам точно воспроизводит последовательные сумму и произведение.
Следовательно, итоговое выражение точно равно коэффициенту при \(x^{10}\) в требуемой унитарной функциональной обратной наименьшей степени.
В C++ Block.product[d] представляет коэффициент усечённого произведения, а Block.sum[d] — соответствующую взвешенную сумму. Последовательная специализация calculate_block сохраняет только sum[9]; параллельные потоки хранят всю сумму. Python и Java сохраняют то же произведение в coefficients, а девятый коэффициент суммы — соответственно в interpolation_coefficient и interpolationCoefficient. Во всех трёх реализациях математическое состояние одинаково.
В начале последовательной итерации \(k\) массив содержит \(A_{k-1}\), накопитель — \(s_9^{(k-1)}\), а переменная факториала — \((k-1)!\). Берём вес из пакета, прибавляем его вклад, обновляем массив по убыванию степени и умножаем факториал на \(k\). Эти действия восстанавливают тот же инвариант для итерации \(k+1\). После последней итерации состояние в точности равно \((A_N,s_9^{(N)},N!)\), то есть содержит всё необходимое для итоговой формулы.
По умолчанию C++ использует не более \(16\) рабочих потоков; необязательный аргумент --threads позволяет выбрать другое поддерживаемое число. Потоки обрабатывают последовательные интервалы и объединяются в исходном порядке. Python и Java вычисляют последовательную рекуррентную формулу. Бинарное модульное возведение в степень выполняют power в C++, трёхаргументный pow в Python и modPow в Java. Для данного модуля произведение двух приведённых остатков меньше \((10^9+7)^2<2^{63}\), поэтому 64-битной арифметики C++ и типа long в Java достаточно; 32-битного произведения недостаточно.
Для проверки каждая реализация также содержит прямую интерполяцию Лагранжа, строящую весь многочлен при \(1\le N\le40\). Это другая конструкция, не использующая оптимизированную рекуррентную формулу Ньютона. Сравнения с эталоном включают уже приведённый случай \(N=1\) и проверяют коэффициент десятой степени даже тогда, когда у малого многочлена такого члена нет. Самопроверки C++ дополнительно сравнивают последовательные и параллельные результаты около границ пакетов обращения; Python и Java также проверяют намеренно малый размер пакета.
Пусть \(D=9\), размер пакета обращения \(B=32768\), а число рабочих потоков равно \(T\); при последовательном вычислении \(T=1\). Каждый индекс требует обновления \(D\) коэффициентов произведения и постоянного числа других операций; полная сумма рабочего потока добавляет ещё \(D+1\) обновлений. Префиксный и обратный проходы линейны по размеру пакета. Общее число пакетных обращений не превосходит \(\lceil N/B\rceil+T-1\), а проверка старшего коэффициента требует ещё одного обращения.
\[W(N,D,B,T,p)=O\!\left(ND+\left(\frac NB+T\right)\log p+TD^2\right).\]
Член \(TD^2\) учитывает короткие полиномиальные свёртки при объединении. При фиксированных \(D,B,T\) и данном фиксированном модуле общая работа составляет \(O(N)\). Для целевого \(N\) последовательный запуск выполняет \(90{,}000{,}000\) обновлений коэффициентов произведения и лишь \(306\) пакетных обращений плюс итоговое обращение для старшего коэффициента. Параллелизм делит основной цикл между потоками, не меняя порядка общей работы; реальное время зависит также от планирования и оборудования.
Дополнительная память составляет \(O(B+D)\) в последовательной реализации и \(O(T(B+D))\) в параллельной, без учёта стеков потоков среды выполнения. Нет ни массива на все \(N\) узлов, ни массива коэффициентов до степени \(N\). Тот же вывод позволяет вычислить коэффициент другой положительной степени \(r\), сохраняя степени от \(0\) до \(r-1\); требование десятой степени ограничивает текущий массив десятью коэффициентами.
لمعيار أولي \(M\)، يُقصد بالمعكوس الدالي \((N,M)\) للدالة \(x^2\) كثير الحدود الأحادي \(Q(x)\) ذو أصغر درجة، على أن تُمثَّل معاملاته بعناصر \(\{0,\ldots,M-1\}\) وأن يحقق
\[Q(n^2)\equiv n\pmod M\qquad(0\le n\le N).\]
في هذه المسألة \(N=10^7\) و\(M=p=10^9+7\)، والمطلوب هو معامل \(x^{10}\) فقط. لا حاجة إلى بناء كثير حدود الاستيفاء ذي الدرجة \(N\)، ولا حتى إلى حفظ جميع معاملاته؛ إذ يمكن قطع الحساب عند الدرجة \(9\) قطعًا دقيقًا ومبرهنًا.
يرمز \([x^r]F(x)\) إلى معامل \(x^r\) في \(F\). تُفهم جميع المساويات التالية داخل \(\mathbb F_p\)، إلا إذا ذُكر صراحة أنها متطابقات على الأعداد النسبية. الكسر \(a/b\) يعني \(a\,b^{-1}\pmod p\)، وليس قسمة صحيحة أو قسمة بالفاصلة العائمة. بما أن \(p\) أولي، فلكل مقام غير صفري معكوس وحيد. واختيار ممثل النتيجة بين \(0\) و\(p-1\) يعطي صيغة المعاملات المطلوبة في المسألة.
لا تطلب المسألة عكس التربيع على الحقل كله، لأن \(n\) و\(-n\) لهما المربع نفسه. إنها تحدد قيمة واحدة عند كل عقدة مربعة مختارة فقط. ويضمن الشرط \(2N<p\) ألّا يكون فهرسان موجبان مختاران جذرين تربيعيين متعاكسين للعقدة نفسها. لذلك تكون نقطة البداية المناسبة هي الاستيفاء، لا متسلسلة قوى الجذر التربيعي الحقيقي.
نعمل في الحقل \(\mathbb F_p\). العقد \(0^2,1^2,\ldots,N^2\) متمايزة: فلو كان \(i^2\equiv j^2\pmod p\)، للزم أن \(p\mid(i-j)(i+j)\)، بينما \(0\le i,j\le N\) و\(i+j\le2N<p\). لذلك يوجد كثير حدود وحيد \(P_N\) درجته لا تتجاوز \(N\) ويحقق \(P_N(n^2)=n\) عند كل العقد المطلوبة.
نكتبه بصيغة نيوتن:
\[P_N(x)=\sum_{k=1}^{N}d_k\prod_{j=0}^{k-1}(x-j^2).\]
حد \(k=0\) يساوي صفرًا لأن القيمة المفروضة عند الصفر هي صفر أيضًا.
يمكن إثبات الوجود والوحدانية صراحة. عند تثبيت \(k\)، نعرّف أساس لاغرانج التالي. عند \(x=j^2\)، يساوي عنصره ذو الفهرس \(i\) العدد \(1\) إذا كان \(i=j\)، ويساوي \(0\) خلاف ذلك. لذلك يحقق مجموعها الموزون القيم المطلوبة. ولو حققها كثيرا حدود درجتهما لا تتجاوز \(k\)، لكان لفرقهما \(k+1\) جذور مختلفة، وهذا مستحيل لكثير حدود غير صفري من هذه الدرجة.
\[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).\]
تضيف صيغة نيوتن العقد واحدة تلو الأخرى. ينعدم الفرق \(P_k-P_{k-1}\) عند العقد \(k\) الأولى، ولذلك فهو مضاعف ثابت لحاصل ضربها الأحادي \(B_k\). ويحدد التعويض بالعقدة الجديدة هذا الثابت. ولأن درجة \(P_{k-1}\) لا تتجاوز \(k-1\)، فإن \(d_k\) هو أيضا معامل \(x^k\) في صيغة لاغرانج لـ\(P_k\). واستخراج هذا المعامل الرئيس يعطي مجموع الفروق المقسومة المستخدم فيما يلي.
\[\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}\]
فرق القسمة على العقد \(0^2,\ldots,k^2\) هو
\[d_k=\sum_{i=1}^{k}\frac{i} {\displaystyle\prod_{\substack{0\le j\le k\\j\ne i}}(i^2-j^2)}.\]
لحساب المقام، نفصل حاصلي ضرب العوامل الخطية. يحتوي الأول على \(k-i\) عوامل سالبة بالضبط. وفي الثاني يجب حذف العامل المستبعد \(2i\) من حاصل ضرب الأعداد المتتالية من \(i\) إلى \(i+k\). لذلك، عندما \(1\le i\le k\)، نحصل على
\[\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}\]
عند ضرب التعبيرين، يختزل \(i!/(i-1)!=i\) مع العامل \(2i\). وجميع المضروبات الظاهرة قابلة للعكس بترديد \(p\)، لأن مدخلاتها لا تتجاوز \(2k\le2N<p\). ومن ثم
\[\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)!}.\]
يمكن حساب المجموع المتبقي دون إغفال خطوات الاختزال. نضرب في \((2k)!\)، ونضع \(r=k-i\)، ونسمي النتيجة \(S_k\). عندئذ يتحول المقام المحتوي على المضروبات إلى معامل ثنائي الحدين:
\[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).\]
لمعالجة المجاميع متناوبة الإشارة، نطبق متطابقة باسكال على كل حد. وبعد إزاحة فهرس المجموع الثاني بمقدار واحد، تتلاشى جميع الحدود الداخلية ولا يبقى إلا الحد الأخير. ونعتمد أن معامل ثنائي الحدين ذي الفهرس السفلي \(-1\)، وكذلك المجموع الفارغ، يساويان صفرا.
\[\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}\]
نعوض فرق معاملي ثنائي الحدين في \(S_k\). تتحول إشارة الطرح إلى جمع بعد إزاحة الفهرس، لأن الإشارة المتناوبة تتغير أيضا. نطبق متطابقة المجموع الجزئي مرتين، ثم نستخدم نسبة معاملي ثنائي الحدين المتجاورين. وبفضل اتفاقية المجموع الفارغ، يصح الحساب أيضا عندما \(k=1\):
\[\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)!}}.\]
نعرّف
\[A_k(x)=\prod_{j=1}^{k}\left(1-\frac{x}{j^2}\right), \qquad A_0(x)=1.\]
وبما أن
\[\prod_{j=0}^{k-1}(x-j^2) =x(-1)^{k-1}((k-1)!)^2A_{k-1}(x),\]
تُختزل الإشارات والمضروبات مع \(d_k\)، فنحصل على المتطابقة البسيطة
\[\boxed{\frac{P_N(x)}x= \sum_{k=1}^{N}\frac{A_{k-1}(x)}{k(2k-1)}}.\]
إذن \([x^{10}]P_N=[x^9](P_N/x)\). لنضع \(a_d^{(k)}=[x^d]A_k\)، وليرمز \(s_d^{(k)}\) إلى المعامل المتراكم بعد الحدود \(1,\ldots,k\). لكل \(1\le d\le9\)، تكون الخطوة
\[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}.\]
تُحدَّث مصفوفة \(a_d\) تنازليًا من الدرجة \(9\) إلى \(1\)، ولذلك يبقى الطرف الأيمن معاملًا من \(A_{k-1}\). حذف الدرجات الأعلى من \(9\) دقيق وليس تقريبًا: فالضرب اللاحق لا يحوي درجات سالبة، ولا يستطيع إعادة معامل محذوف إلى الدرجة \(9\).
بصورة صريحة، العامل المضروب في \(A_{k-1}\) هو \(d_k(-1)^{k-1}((k-1)!)^2=1/(k(2k-1))\). نهيئ حالة المعاملات والمجموع كما يلي. يبقى الحد الثابت لـ\(A_k\) مساويا \(1\) دائما، فلا حاجة إلى تحديثه.
\[a_0^{(0)}=1,\qquad a_d^{(0)}=0\ (1\le d\le9),\qquad s_d^{(0)}=0\ (0\le d\le9).\]
هناك شرطان مستقلان لترتيب العمليات. نجمع أولا مساهمة \(A_{k-1}\)، ثم نضرب بالعامل الجديد لنحصل على \(A_k\). وداخل هذا الضرب نحدّث الدرجات تنازليا. أما التحديث التصاعدي فيستخدم \(a_{d-1}\) بعد تغييره. فمثلا، عند البدء من \(1\)، سيولد خطأ الحد \(x^2/k^4\) بعد الضرب في \(1-x/k^2\) وحده، رغم أن حاصل الضرب من الدرجة الأولى.
التبرير الشكلي لصحة القطع هو أننا نحسب في حلقة القسمة \(\mathbb F_p[x]/(x^{10})\). الاختزال بترديد \(x^{10}\) متوافق مع الجمع والضرب. وبصيغة مكافئة، يعتمد معامل الدرجة \(d\le9\) في حاصل الضرب على معاملات المدخلات ذات الدرجات التي لا تتجاوز \(d\) فقط. لذلك لا يؤثر أي حد محذوف ذي درجة عالية في المعاملات المخزنة:
\[[x^d](F(x)G(x))=\sum_{r=0}^{d}([x^r]F(x))([x^{d-r}]G(x)),\qquad 0\le d\le9.\]
لذلك يخزن الحساب التسلسلي عشرة معاملات لحاصل الضرب، لكنه يحتاج إلى العدد \(s_9\) وحده للمجموع. فمعاملات المجموع الأخرى لا تدخل في تحديث \(s_9\). ومع ذلك تبقى معاملات حاصل الضرب العشرة كلها ضرورية، لأن التحديثات المتتابعة تنقل المعلومات من الحد الثابت عبر الدرجات الوسيطة إلى الدرجة التاسعة. كذلك يكون \(a_9^{(k-1)}=0\) عندما \(k<10\)، فلا تسهم الحدود التسعة الأولى في جزء الاستيفاء من المعامل المطلوب.
لنأخذ \(N=3\)، ونحتفظ مؤقتا بالكسور النسبية لتوضيح الجبر. يمكن بعد ذلك اختزال هذه المتطابقات بترديد أي عدد أولي \(p>6\). حواصل الضرب المطَبَّعة الأولى هي
\[A_0=1,\qquad A_1=1-x,\qquad A_2=(1-x)(1-x/4)=1-\frac54x+\frac14x^2.\]
الأوزان عند \(k=1,2,3\) هي \(1,1/6,1/15\). يستخدم كل وزن حاصل الضرب قبل إدراج عامله الخاص. بالتعويض وتجميع المعاملات نحصل على
\[\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.\]
يتحقق بذلك كل من القيم المطلوبة وصيغة المعامل \(d_3=1/60\). ويوضح المثال أيضا الفرق بين الاستيفاء وشرط أن يكون المعامل الرئيس واحدا: عند استخدام معامل الترديد المستهدف، يكون المعامل الرئيس لهذا المثال الصغير \(1/60\ne1\). ولذلك لا يكفي الاستيفاء وحده لإتمام البناء. سنشتق التصحيح الذي يجعل كثير الحدود أحاديا فيما يلي.
لنضع \(u_k=k(2k-1)\). كل \(u_k\) غير صفري في \(\mathbb F_p\)، لأن \(1\le k\le N<p\) و\(2k-1<p\). في دفعة \(u_1,\ldots,u_m\)، نبني نواتج البادئة \(U_0=1\) و\(U_i=U_{i-1}u_i\). نحسب \(U_m^{-1}=U_m^{p-2}\) مرة واحدة فقط، ثم نسير إلى الخلف:
\[u_i^{-1}=U_{i-1}U_i^{-1}, \qquad U_{i-1}^{-1}=U_i^{-1}u_i.\]
وهكذا تعطي عملية أس واحدة بمبرهنة فيرما كل مقلوبات الدفعة. كذلك \(k^{-1}=(2k-1)u_k^{-1}\)، وتربيعه يعطي \(k^{-2}\) اللازم لتحديث الناتج. تستخدم التطبيقات دفعات من \(32768\) عنصرًا.
الفهارس في صيغة العكس الدفعي محلية للدفعة، وليس ضروريا أن يكون فهرسها العام الأول \(1\). أثناء المرور العكسي نحافظ على أن المعكوس الجاري يساوي \(U_i^{-1}\). ضربه في البادئة \(U_{i-1}\) التي لم تُعدَّل بعد يعطي \(u_i^{-1}\)، وضرب المعكوس الجاري في \(u_i\) يهيئ التكرار التالي. ويمكن عندئذ استبدال خانة البادئة بالمعكوس المحسوب، لأن التكرارات التالية لا تحتاج إلا إلى الخانات السابقة.
للتحقق من هذه الحيلة العامة بمثال ملموس، نحسب معكوسات \(2,3,5\) بترديد \(11\). حواصل ضرب البادئات هي \(1,2,6,8\)، ولدينا \(8^{-1}=7\). يعطي المرور العكسي أولا \(5^{-1}=6\cdot7=9\)، ثم \(3^{-1}=2\cdot2=4\)، وأخيرا \(2^{-1}=1\cdot6=6\)، وكل ذلك بترديد \(11\). لم نحتج إلا إلى عملية رفع إلى قوة واحدة. وفي الخوارزمية الفعلية تكون المدخلات هي القيم غير الصفرية \(k(2k-1)\)، ويُختزل كل حاصل ضرب بترديد \(p\).
كل كثير حدود أساسي في صيغة نيوتن أحادي، لذا فإن المعامل الرئيس لـ\(P_N\) هو
\[c_N=d_N= \frac{(-1)^{N-1}N}{(2N-1)(N!)^2}.\]
وكثير الحدود المنعدم عند جميع عقد الاستيفاء هو
\[V_N(x)=\prod_{j=0}^{N}(x-j^2) =(-1)^N(N!)^2xA_N(x).\]
قيمة \(d_N\) غير صفرية: بسطها مجرد إشارة، ولا ينعدم أي عامل في مقامها. لذلك فإن درجة \(P_N\) تساوي \(N\) بالضبط، وليست فقط أقل منها أو مساوية لها. لوصف جميع الحلول، ليكن \(F\) أي كثير حدود يحقق القيم المطلوبة. يقسم كل عامل خطي \(x-j^2\) الفرق \(F-P_N\)، وبما أن الجذور مختلفة، فإن حاصل ضرب هذه العوامل يقسمه أيضا. ومن ثم
\[F(x)=P_N(x)+V_N(x)R(x),\qquad R(x)\in\mathbb F_p[x].\]
إذا كان \(R\ne0\)، فإن \(\deg(V_NR)=N+1+\deg R>\deg P_N\)، ولذلك لا يمكن إلغاء الحد الأعلى درجة. إذن الحل الوحيد الذي لا تتجاوز درجته \(N\) هو \(P_N\). وإذا لم يكن أحاديا، فإن الحصول على حل أحادي من الدرجة التالية يقتضي أن يكون \(R\) ثابتا، ويفرض المعامل الرئيس أن يساوي هذا الثابت \(1\) بالضبط. وهذا يثبت أصغرية الدرجة والوحدانية معا. أما قسمة \(P_N\) مباشرة على \(c_N\) فتغير القيمة المطلوبة عند \(1\) من \(1\) إلى \(c_N^{-1}\)، لذا لا يصلح التحجيم تصحيحا.
إذا كان \(c_N=1\)، فإن كثير حدود الاستيفاء نفسه هو الحل الأحادي الوحيد ذو أصغر درجة. وإلا فلا يوجد حل آخر درجته لا تتجاوز \(N\)، وتكون أصغر درجة أحادية ممكنة \(N+1\)؛ وبإضافة كثير الحدود الأحادي \(V_N\) نحصل على
\[Q(x)=P_N(x)+V_N(x).\]
يتحقق البرنامج في الحالة المطلوبة من أن \(c_N\ne1\). ومن ثم يُستخرج المعامل المطلوب من دون بناء أي من كثيري الحدود الكبيرين:
\[\boxed{[x^{10}]Q=s_9^{(N)} +(-1)^N(N!)^2a_9^{(N)}}.\]
ملاحظة حول المثال الصغير المنشور. عند قراءة عبارة «أصغر درجة» حرفيًا، تظهر حالة حدية خاصة عند \((N,M)=(2,7)\): ففيها \(c_2=1\)، ولذلك يكون \(P_2(x)=x^2\) أحاديًا بالفعل ويحقق القيم الثلاث. أما كثير الحدود التكعيبي المعروض فهو \(P_2+V_2\)، وهو يحققها أيضًا. هذه الملاحظة التحريرية لا تؤثر في الحالة المطلوبة، حيث \(c_N\ne1\).
لمثال مستقل لا تظهر فيه الحالة الخاصة للمعامل الرئيس، نأخذ \((N,p)=(2,11)\). يُختزل كثير حدود الاستيفاء النسبي \((7x-x^2)/6\) إلى \(3x+9x^2\). وبإضافة كثير الحدود المنعدم عند العقد نحصل على كثير الحدود التكعيبي التالي. قيمه عند \(0,1,4\) هي على الترتيب \(0,1,2\) بترديد \(11\)، ومعامله الرئيس \(1\). يوضح ذلك كيف يغير التصحيح المعاملات مع الحفاظ على جميع القيم المطلوبة.
\[\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}\]
يستطيع برنامج C++ معالجة فترات متتالية من الفهارس بالتوازي. لكتلة \(B\)، نكتب \(g_k(x)=1-x/k^2\) ونحفظ حتى الدرجة \(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.\]
إذا سبقت الكتلة \(L\) الكتلة \(R\) مباشرة، فإن ملخصيهما يندمجان وفق
\[\boxed{G_{LR}=G_LG_R,\qquad H_{LR}=H_L+G_LH_R}.\]
هذه العملية تجميعية، لكن ترتيبها مهم؛ لذلك يدمج التطبيق الكتل حسب تزايد الفهارس. وتُضرب نواتج المضروبات في الوقت نفسه. يكفي ضرب كثيرات الحدود المقطوع، لأن المعاملات فوق الدرجة \(9\) لا يمكنها الإسهام في درجة أدنى.
يمكن قراءة قاعدة الكتل أيضا كتحديث للحالة العامة الداخلة: \((S,A)\mapsto(S+AH_B,AG_B)\). تحتاج كل مساهمة من الكتلة اليمنى إلى جميع عوامل الكتلة اليسرى قبلها، وهذا يفسر الحد \(G_LH_R\). ولثلاث كتل متتالية، يعطي كلا ترتيبي الأقواس التعبير نفسه أدناه، مما يثبت التجميعية. وملخص الكتلة الفارغة هو \((G,H)=(1,0)\).
\[G_{LMR}=G_LG_MG_R,\qquad H_{LMR}=H_L+G_LH_M+G_LG_MH_R.\]
يفسر هذا أيضا ضرورة احتفاظ كل خيط تنفيذ بالمعاملات العشرة لـ\(H_B\)، رغم أن الجواب النهائي يحتاج إلى الدرجة التاسعة فقط. فحاصل الضرب السابق يمزج الدرجات وفق الالتفاف أدناه. ويمكن لشفرة الدمج الاحتفاظ بالمعامل التاسع وحده للمجموع العام، ولكن ليس بالمعامل التاسع وحده لمجموع كل خيط. دفعات حساب المعكوس والكتل المتوازية تؤدي أغراضا مختلفة: يستطيع الخيط معالجة دفعات عكس عديدة داخل مجاله المتتالي.
\[[x^9](AH_B)=\sum_{r=0}^{9}([x^r]A)([x^{9-r}]H_B).\]
اللمّة 1. العقد المربعة متمايزة في \(\mathbb F_p\)، لذلك يوجد \(P_N\) وهو كثير حدود الاستيفاء الوحيد الذي لا تتجاوز درجته \(N\).
اللمّة 2. يعطي حساب مقام فرق القسمة ومجموع ثنائي الحدين المتناوب القيمة المذكورة لـ\(d_k\). وبالتعويض في صيغة نيوتن نحصل على \(P_N/x=\sum A_{k-1}/(k(2k-1))\).
اللمّة 3. يحافظ التكرار التنازلي على \(a_d^{(k)}=[x^d]A_k\)، بينما يحافظ التجميع على \(s_d^{(k)}=[x^d]\sum_{i=1}^{k}A_{i-1}/(i(2i-1))\). ويحفظ القطع عند الدرجة \(9\) هذين الثابتين بدقة.
اللمّة 4. يعيد العكس على دفعات المقلوب الصحيح لكل \(u_k\ne0\)، لأن \(U_{i-1}U_i^{-1}=1/u_i\). فهو لا يغير إلا طريقة تقييم التكرار نفسه.
اللمّة 5. يختلف كل كثير حدود يحقق البيانات عن \(P_N\) بمضاعف لـ\(V_N\). لذلك يختار اختبار المعامل الرئيس \(P_N\) بالضبط حين يكون أحاديًا أصلًا، وإلا فيختار \(P_N+V_N\)، وهو أحادي وذو أصغر درجة ممكنة.
اللمّة 6. تضع متطابقة الكتل المرتبة \(G_L\) قبل كل مساهمة آتية من \(R\). ومن ثم يعيد الاستقراء على الكتل المدمجة تمامًا المجموع والناتج التسلسليين.
وبجمع هذه اللمّات، يكون التعبير النهائي بالضبط معامل \(x^{10}\) في المعكوس الدالي الأحادي المطلوب ذي أصغر درجة.
في C++ يمثل Block.product[d] معامل حاصل الضرب المقطوع، ويمثل Block.sum[d] المجموع الموزون المقابل. يحتفظ التخصص التسلسلي للدالة calculate_block بـsum[9] فقط، بينما تحتفظ الخيوط المتوازية بالمجموع كاملا. يخزن Python وJava حاصل الضرب نفسه في coefficients، ومعامل الدرجة التاسعة للمجموع في interpolation_coefficient وinterpolationCoefficient على الترتيب. الحالة الرياضية واحدة في التطبيقات الثلاثة.
عند بداية التكرار التسلسلي \(k\)، تحتوي المصفوفة على \(A_{k-1}\)، والمجمّع على \(s_9^{(k-1)}\)، ومتغير المضروب على \((k-1)!\). نستخرج وزن الدفعة، ونضيف مساهمته إلى المجمّع، ونحدّث المصفوفة تنازليا، ثم نضرب متغير المضروب في \(k\). تعيد هذه العمليات تحقيق الثابت نفسه للتكرار \(k+1\). وبعد التكرار الأخير تصبح الحالة بالضبط \((A_N,s_9^{(N)},N!)\)، وهي كل ما تحتاج إليه الصيغة النهائية.
يستخدم C++ افتراضيا ما لا يزيد على \(16\) خيط تنفيذ، ويمكن للوسيط الاختياري --threads اختيار عدد آخر مدعوم. تعالج الخيوط مجالات متتالية، ثم تُدمج بترتيبها الأصلي. أما Python وJava فينفذان العلاقة التكرارية تسلسليا. تنفذ دوال power في C++، وpow ذات الوسائط الثلاثة في Python، وmodPow في Java الرفع الثنائي إلى قوة بترديد. وللمعامل المستخدم هنا، يكون حاصل ضرب باقيين مختزلين أصغر من \((10^9+7)^2<2^{63}\)، لذا تكفي حسابات C++ ذات 64 بت ونوع long في Java، بينما لا يكفي ضرب 32 بت.
للتحقق، يحتوي كل تطبيق أيضا على استيفاء مباشر بطريقة لاغرانج يبني كثير الحدود كاملا عندما \(1\le N\le40\). وهو يستخدم بناء مختلفا عن علاقة نيوتن المحسنة. تشمل المقارنات مع هذا المرجع الحالة الأحادية أصلا عند \(N=1\)، وتفحص معامل الدرجة العاشرة حتى عندما لا يحتوي كثير الحدود الصغير على هذا الحد. وتقارن اختبارات C++ الذاتية أيضا النتائج التسلسلية والمتوازية قرب حدود دفعات العكس، بينما يختبر Python وJava حجما صغيرا للدفعة اختير عمدا.
ليكن \(D=9\)، وحجم دفعة حساب المعكوس \(B=32768\)، وعدد خيوط التنفيذ \(T\)، حيث \(T=1\) في الحساب التسلسلي. يحتاج كل فهرس إلى تحديث \(D\) معاملات لحاصل الضرب وكمية ثابتة من الحسابات الأخرى، ويضيف المجموع الكامل للخيط \(D+1\) تحديثات. وتكلفة مروري البادئات والرجوع خطية في حجم الدفعة. لا يتجاوز العدد الكلي لعمليات عكس الدفعات \(\lceil N/B\rceil+T-1\)، ويستخدم اختبار المعامل الرئيس عملية عكس إضافية.
\[W(N,D,B,T,p)=O\!\left(ND+\left(\frac NB+T\right)\log p+TD^2\right).\]
يمثل الحد \(TD^2\) التفافات كثيرات الحدود القصيرة عند الدمج. وعند تثبيت \(D,B,T\) ومعامل الترديد هذا، تكون كمية العمل الكلية \(O(N)\). للقيمة المستهدفة من \(N\)، ينفذ الحساب التسلسلي \(90{,}000{,}000\) تحديث لمعاملات حاصل الضرب، و\(306\) عمليات عكس دفعات فقط، إضافة إلى عملية العكس الأخيرة للمعامل الرئيس. توزع الموازاة الحلقة الرئيسية على الخيوط دون تغيير رتبة العمل الكلي، كما يعتمد الزمن الفعلي على الجدولة والعتاد.
تبلغ الذاكرة الإضافية \(O(B+D)\) في التطبيق التسلسلي و\(O(T(B+D))\) في المتوازي، باستثناء مكدسات خيوط بيئة التشغيل. لا توجد مصفوفة مفهرسة على جميع العقد \(N\)، ولا مصفوفة معاملات حتى الدرجة \(N\). ويمكن للاشتقاق نفسه حساب معامل أي درجة موجبة أخرى \(r\) بالاحتفاظ بالدرجات من \(0\) إلى \(r-1\)؛ وطلب الدرجة العاشرة هنا هو ما يجعل الحالة مقتصرة على عشرة معاملات.