vi. Multipoint Evaluation and Interpolation

Converting between coefficients and values: evaluation at $n$ arbitrary points and interpolation in $O(n{\log }^{2}n)$ with a subproduct tree, $O(n)$ Lagrange interpolation at consecutive points (sums of powers, “the answer is a polynomial in $n$”), shifting sample points, $n!\mathrm{mod}p$ in $\tilde{O}(\sqrt{n})$, and the Newton basis. Uses Poly, polyMul, polyDivMod, derivative, evalPoly from NTT and Power Series, Comb from Modular Arithmetic.

Task	Naive	Fast
evaluate at 1 point	$O(n)$ Horner
evaluate at $n$ consecutive points / shift samples $f(\mathrm{0..}n-1)\to f(c..c+m-1)$	$O(n^2)$	$O((n+m)\log )$, one convolution
evaluate at $n$ points in geometric progression	$O(n^2)$	$O(n\log n)$ chirp-z (Convolution Variants)
evaluate at $n$ arbitrary points	$O(n^2)$	$O(n{\log }^{2}n)$ subproduct tree
interpolate from $n$ arbitrary points (coefficients)	$O(n^2)$	$O(n{\log }^{2}n)$
value at one new point from consecutive samples	$O(n^2)$	$O(n)$ Lagrange
shift $P(x)\to P(x+c)$	$O(n^2)$	$O(n\log n)$ (Taylor Shift)

Check the point structure first

Consecutive points give $O(n)$ interpolation and $O(n\log n)$ evaluation; geometric points give chirp-z; only genuinely arbitrary points need the $O(n{\log }^{2}n)$ tree, whose constant (several NTTs per node) makes it slower than $O(n^2)$ Horner below $n\approx 5000$. In contests the points are usually $0,1,2,\dots$.

vector<Poly> subproductTree(const vector<mint>& xs) {       // tree[1] = prod (x - x_i); node p covers a contiguous range
    int n = sz(xs); vector<Poly> t(4 * max(n, 1));
    auto build = [&](auto&& self, int p, int l, int r) -> void {
        if (r - l == 1) { t[p] = {-xs[l], 1}; return; }
        int m = (l + r) / 2; self(self, 2 * p, l, m); self(self, 2 * p + 1, m, r);
        t[p] = polyMul(t[2 * p], t[2 * p + 1]);
    };
    if (n) build(build, 1, 0, n);
    return t;
}
vector<mint> multipointEval(const Poly& f, const vector<mint>& xs) {   // O(n log^2 n): remainders down the tree
    int n = sz(xs); if (n == 0) return {};
    auto t = subproductTree(xs); vector<mint> res(n);
    auto down = [&](auto&& self, int p, int l, int r, Poly g) -> void {
        g = polyDivMod(g, t[p]).second;                     // g = f mod prod_{i in [l,r)} (x - x_i)
        if (r - l <= 32) { for (int i = l; i < r; i++) res[i] = evalPoly(g, xs[i]); return; }
        int m = (l + r) / 2; self(self, 2 * p, l, m, g); self(self, 2 * p + 1, m, r, g);
    };
    down(down, 1, 0, n, f);
    return res;
}
Poly interpolate(const vector<mint>& xs, const vector<mint>& ys) {   // distinct xs, O(n log^2 n)
    int n = sz(xs); if (n == 0) return {};
    auto t = subproductTree(xs);
    vector<mint> d = multipointEval(derivative(t[1]), xs);  // M'(x_i) = prod_{j != i} (x_i - x_j)
    auto up = [&](auto&& self, int p, int l, int r) -> Poly {
        if (r - l == 1) return {ys[l] / d[l]};
        int m = (l + r) / 2;
        Poly a = polyMul(self(self, 2 * p, l, m), t[2 * p + 1]), b = polyMul(self(self, 2 * p + 1, m, r), t[2 * p]);
        if (sz(a) < sz(b)) swap(a, b); for (int i = 0; i < sz(b); i++) a[i] += b[i];
        return a;
    };
    Poly res = up(up, 1, 0, n); res.resize(n); return res;
}
 
mint lagrangeGeneral(const vector<mint>& xs, const vector<mint>& ys, mint t) {   // O(n^2), one inverse
    int n = sz(xs); mint num = 0, den = 1; vector<mint> dens(n);
    for (int i = 0; i < n; i++) { if (xs[i] == t) return ys[i]; mint d = t - xs[i]; for (int j = 0; j < n; j++) if (j != i) d *= xs[i] - xs[j]; dens[i] = d; }
    vector<mint> pre(n + 1, 1); for (int i = 0; i < n; i++) pre[i + 1] = pre[i] * dens[i];   // batch inverse
    mint inv = pre[n].inv(), full = 1; for (int i = 0; i < n; i++) full *= t - xs[i];
    for (int i = n - 1; i >= 0; i--) { num += ys[i] * pre[i] * inv; inv *= dens[i]; }
    return num * full;
}
mint lagrangeConsecutive(const vector<mint>& y, ll k, const Comb& cb) {   // samples at x = 0..n-1, value at k (any ll), O(n)
    int n = sz(y); if (k >= 0 && k < n) return y[k];
    vector<mint> pre(n + 1, 1), suf(n + 1, 1);
    for (int i = 0; i < n; i++) pre[i + 1] = pre[i] * mint(k - i);
    for (int i = n - 1; i >= 0; i--) suf[i] = suf[i + 1] * mint(k - i);
    mint res = 0;
    for (int i = 0; i < n; i++) { mint term = y[i] * pre[i] * suf[i + 1] * cb.fi[i] * cb.fi[n - 1 - i]; res += (n - 1 - i) % 2 ? -term : term; }
    return res;
}
mint powerSum(ll n, int k, const Comb& cb) {               // sum_{i=1}^{n} i^k, O(k log) for any n
    vector<mint> y(k + 2, 0);
    for (int i = 1; i <= k + 1; i++) y[i] = y[i - 1] + mint(i).pow(k);
    return lagrangeConsecutive(y, n, cb);
}
Poly shiftSamples(const Poly& y, mint c, int m, const Comb& cb) {   // f(0..n-1) -> f(c), ..., f(c+m-1) for deg f < n; needs c + i not in [0, n) and cb >= n
    int n = sz(y); Poly a(n), b(n + m - 1);
    for (int i = 0; i < n; i++) { a[i] = y[i] * cb.fi[i] * cb.fi[n - 1 - i]; if ((n - 1 - i) % 2) a[i] = -a[i]; }
    vector<mint> v(n + m - 1); for (int j = 0; j < n + m - 1; j++) v[j] = c - mint(n - 1) + mint(j);   // c - (n-1) + j
    vector<mint> pre(n + m, 1); for (int j = 0; j < n + m - 1; j++) pre[j + 1] = pre[j] * v[j];
    mint inv = pre[n + m - 1].inv();
    for (int j = n + m - 2; j >= 0; j--) { b[j] = inv * pre[j]; inv *= v[j]; }   // b_j = 1 / (c - (n-1) + j)
    Poly conv = polyMul(a, b), res(m);
    mint prod = 1; for (int j = 0; j < n; j++) prod *= c - mint(j);   // prod_{j<n} (c + 0 - j)
    for (int k = 0; k < m; k++) {
        res[k] = conv[n - 1 + k] * prod;                    // f(c+k) = prod_{j<n}(c+k-j) * sum_i a_i / (c+k-i)
        prod *= (c + mint(k + 1)) / (c + mint(k + 1) - mint(n));
    }
    return res;
}
mint factorialMod(ll n) {                                   // n! mod 998244353, O(sqrt n log^2 n); n < p
    if (n <= 1) return 1;
    ll v = max(1LL, (ll)sqrtl((long double)n));
    vector<mint> rs(v); for (ll i = 0; i < v; i++) rs[i] = -mint(i + 1);   // f(x) = (x+1)...(x+v)
    Poly f = subproductTree(rs)[1];
    vector<mint> xs(v); for (ll i = 0; i < v; i++) xs[i] = mint(i * v);
    mint res = 1; for (auto val : multipointEval(f, xs)) res *= val;   // covers 1..v^2
    for (ll j = v * v + 1; j <= n; j++) res *= mint(j);
    return res;
}

Proof: subproduct tree, interpolation, Lagrange, sample shifting, factorial

Evaluation. By the remainder theorem $f(x_i)=f\mathrm{mod}(x-x_i)$, and for $S^{\prime }\subseteq S$, $(f\mathrm{mod}{M}_S)\mathrm{mod}{M}_{S^{\prime }}=f\mathrm{mod}{M}_{S^{\prime }}$ since $M_{S^{\prime }}\mid M_S$. So descending the tree with remainders yields $f\mathrm{mod}(x-x_i)$ at the leaves (small nodes are finished by Horner). Each level does divisions of total degree $O(n)$, $O(n\log n)$ per level, $O(n{\log }^{2}n)$ in total. Interpolation. Lagrange: $P={\sum }_i\frac{y_i}{M^{\prime }(x_i)}{\prod }_{j\ne i}(x-x_j)$ because ${\prod }_{j\ne i}(x_i-x_j)=M^{\prime }(x_i)$ (product rule on $M=\prod (x-x_j)$). For a node with children $L,R$, ${\sum }_{i\in L\cup R}{c}_i{\prod }_{j\in L\cup R,j\ne i}(x-x_j)=M_R{\sum }_{i\in L}{c}_i{\prod }_{j\in L\setminus i}+M_L{\sum }_{i\in R}{c}_i{\prod }_{j\in R\setminus i}$, which is the bottom-up combination; uniqueness of the degree $<n$ interpolant needs distinct $x_i$.

Consecutive points. With $x_j=j$: ${\prod }_{j\ne i}\frac{k-j}{i-j}=\frac{{\mathrm{pre}}_i{\mathrm{suf}}_{i+1}}{i!(n-1-i)!(-1{)}^{n-1-i}}$. Power sums: $S_k(n)={\sum }_{i\le n}{i}^k$ satisfies $S_k(n)-S_k(n-1)=n^k$, a polynomial of degree $k$, so $S_k$ is a polynomial of degree $k+1$ (finite differences), determined by $k+2$ samples. lagrangeGeneral computes $\prod (t-x_j){\sum }_i\frac{y_i}{(t-x_i){\prod }_{j\ne i}(x_i-x_j)}$ with all denominators inverted by one batch inverse.

Sample shifting. For $\mathrm{deg}f<n$ and $t=c+k\notin \{\mathrm{0..}n-1\}$: $f(t)={\prod }_{j<n}(t-j){\sum }_i\frac{a_i}{t-i}$ with $a_i=\frac{f(i)(-1{)}^{n-1-i}}{i!(n-1-i)!}$ (Lagrange), and ${\sum }_i{a}_i\frac{1}{c+k-i}={\sum }_i{a}_i{b}_{k+n-1-i}$ with $b_j=\frac{1}{c-(n-1)+j}$: coefficient $n-1+k$ of $a\ast b$. The prefactor updates by $\frac{{\prod }_j(c+k+1-j)}{{\prod }_j(c+k-j)}=\frac{c+k+1}{c+k+1-n}$. All the $c-(n-1)+j$ must be non-zero modulo $p$. Factorial: $n!={\prod }_{i<v}f(iv)\cdot {\prod }_{j=v^2+1}^{n}j$ with $f(x)=(x+1)\cdots (x+v)$ ($f(iv)$ is the product of the block $iv+\mathrm{1..}iv+v$), and $v=\lfloor \sqrt{n}\rfloor$ makes both parts $O(\sqrt{n})$ terms; evaluating $f$ at $v$ points is $O(\sqrt{n}{\log }^{2}n)$ (Min_25’s sample shifting gives $O(\sqrt{n}\log n)$). All verified against direct evaluation, brute interpolation, direct power sums and a direct factorial loop up to $10^6$.

What the tree enables
coefficient $\leftrightarrow$ value representations	both directions
fast polynomial gcd (half-gcd), factorisation, Chinese remaindering of polynomials	the same tree shape
$n!\mathrm{mod}p$ for $n\approx 10^{12}$ (with a large NTT-friendly or CRT-split modulus)	$\tilde{O}(\sqrt{n})$
integer CRT from many moduli, batch modular inverse (linear special case), segment trees of products	integer analogues

Recognising “the answer is a polynomial in $n$”

Signals: the answer is a sum of a polynomial over a range (degree goes up by 1); it counts lattice points in a region scaling linearly with $n$ (Ehrhart polynomials); a DP whose answer is eventually polynomial; huge $n$ with a structure of small degree. Verification: compute the first $2d+4$ values by brute force, interpolate from the first $d+2$, and check the predictions on the rest. Pitfalls: repeated $x_i$ (deduplicate), a wrong degree guess (silently wrong answers: verify with extra points), a target equal to a sample point (return it directly).

Newton basis and divided differences

$P(x)={\sum }_i{c}_i{\prod }_{j<i}(x-x_j)$; the $c_i=f[x_0,\dots ,x_i]$ are divided differences (with repeated nodes: confluent / Hermite data). The table costs $O(n^2)$, but adding a point is $O(n)$, useful when the sample set grows. Fast conversion: split the nodes in halves with $P_L={\prod }_{\text{left}}(x-x_j)$; then $f=\mathrm{rem}+P_{L}q$ with $\mathrm{deg}\mathrm{rem}<|L|$, the left coefficients are those of rem and the right ones those of $q$ on the right nodes; backwards, $f=\mathrm{left}+P_L\cdot \mathrm{right}$ (Library Checker “Conversion from Monomial Basis to Newton Basis”).

struct NewtonInterpolator {                                 // add points one at a time, O(n) each; distinct x
    vector<mint> xs, c, diag;                              // diag[j] = f[x_j..x_last]
    void add(mint x, mint y) {
        int n = sz(xs); xs.push_back(x); diag.push_back(y);
        for (int j = n - 1; j >= 0; j--) diag[j] = (diag[j + 1] - diag[j]) / (x - xs[j]);
        c.push_back(diag[0]);
    }
    mint eval(mint t) const { mint r = 0; for (int i = sz(c) - 1; i >= 0; i--) r = r * (t - xs[i]) + c[i]; return r; }
};
vector<mint> toNewtonBasis(const Poly& f, const vector<mint>& xs) {   // deg f < n = |xs|, O(n log^2 n)
    int n = sz(xs); if (n == 0) return {};
    auto t = subproductTree(xs); vector<mint> c(n);
    auto rec = [&](auto&& self, int p, int l, int r, Poly g) -> void {
        g.resize(r - l);
        if (r - l == 1) { c[l] = g[0]; return; }
        int m = (l + r) / 2; auto [q, rem] = polyDivMod(g, t[2 * p]);   // g = rem + P_left q
        self(self, 2 * p, l, m, rem); self(self, 2 * p + 1, m, r, q);
    };
    rec(rec, 1, 0, n, f);
    return c;
}
Poly fromNewtonBasis(const vector<mint>& c, const vector<mint>& xs) {
    int n = sz(xs); if (n == 0) return {};
    auto t = subproductTree(xs);
    auto rec = [&](auto&& self, int p, int l, int r) -> Poly {
        if (r - l == 1) return {c[l]};
        int m = (l + r) / 2; Poly a = self(self, 2 * p, l, m), b = polyMul(t[2 * p], self(self, 2 * p + 1, m, r));
        if (sz(a) < sz(b)) swap(a, b); for (int i = 0; i < sz(b); i++) a[i] += b[i];
        return a;
    };
    Poly res = rec(rec, 1, 0, n); res.resize(n); return res;
}

Proof: divided differences and the Newton basis conversion

Divided differences. Define $f[x_j..x_k]=\frac{f[x_{j+1}..x_k]-f[x_j..x_{k-1}]}{x_k-x_j}$; then $c_i=f[x_0..x_i]$ are the Newton coefficients (induction: the Newton form through $x_0..x_k$ minus the one through $x_0..x_{k-1}$ vanishes at $x_0..x_{k-1}$ and is $c_k{\prod }_{j<k}(x-x_j)$; its value at $x_k$, compared with the form through $x_1..x_k$, gives the recurrence). When point $x_k$ arrives, diag holds $f[x_j..x_{k-1}]$ for all $j$, and the update computes $f[x_j..x_k]$ from $j=k-1$ down, so $c_k=f[x_0..x_k]$ in $O(k)$. Fast conversion: write $f=\mathrm{rem}+P_{L}q$; for $i<|L|$ the basis polynomials ${\prod }_{j<i}(x-x_j)$ have degree $<|L|$ and span the polynomials of degree $<|L|$, while for $i\ge |L|$ they equal $P_L{\prod }_{|L|\le j<i}(x-x_j)$. So the Newton expansion of $f$ splits into that of rem on the left nodes plus $P_L$ times that of $q$ on the right nodes, and both directions recurse; cost $O(n{\log }^{2}n)$. Verified by round trips and against the $O(n^2)$ divided-difference table.

Related notes

• Composition and Taylor Shift
• Factorials mod p
• Counting with Power Series (polynomial sums)