有限域上的O(n log n)时间拉格朗日插值算法

嗨，我来帮你搞定这个问题！针对你在GF(p)上有连续点（0到2e6）的多项式值、需要O(n logn)求系数的需求，其实有几个非常实用的方案，完全不需要受Pan那本书里复数域条件的限制：

一、针对连续点的优化版快速插值（基于NTT） #

因为你的插值点是0到n-1（n≈2e6）这种连续整数，这刚好可以利用阶乘、逆阶乘的性质大幅简化计算，再结合快速数论变换（NTT）实现O(n logn)复杂度，步骤如下：

1. 预处理阶乘与逆阶乘
   先在GF(p)上计算两个数组：

   • fact[i]：i的阶乘模p，递推公式：fact[0] = 1，fact[i] = fact[i-1] * i % p
   • inv_fact[i]：i阶乘的逆元模p，先算inv_fact[n-1] = pow(fact[n-1], p-2, p)（费马小定理），再倒推：inv_fact[i-1] = inv_fact[i] * i % p
     这一步是O(n)时间，完全没问题。

2. 计算每个点的权重
   对于每个i（0≤i<n），拉格朗日基对应的权重w_i可以通过公式直接计算：

   sign = p-1 if (n-1 - i) % 2 else 1  # 等价于(-1)的幂次模p
   w_i = y_i * inv_fact[i] % p
   w_i = w_i * inv_fact[n-1 - i] % p
   w_i = w_i * sign % p
   

   这一步也是O(n)时间。

3. 用NTT做卷积计算系数
   构造两个多项式：

   • A(x) = Σ(w_i * x^i) （i从0到n-1）
   • B(x) = Σ(fact[i] * ((p-1) if i%2 else 1) * x^i) （i从0到n-1）
     用NTT计算它们的卷积C(x) = A(x) * B(x)，这一步是O(n logn)时间。

4. 提取最终系数
   插值多项式P(x)的k次项系数c_k就是：

   c_k = C[n-1 + k] * inv_fact[k] % p
   

   这里的C[m]是卷积结果多项式的m次项系数，k从0到n-1。

二、通用任意点集的O(n logn)插值（如果后续点不连续） #

如果之后你的插值点不是连续整数了，也可以用基于NTT的任意点快速插值算法，核心思路是把拉格朗日插值转化为多项式乘法与除法的组合：

• 先构造所有插值点的乘积多项式L(x) = Π(x - x_i)
• 计算每个点对应的L'(x_i)（多项式L在x_i处的导数），这可以通过预处理点的前缀积、后缀积快速得到
• 构造辅助多项式Q(x) = Σ(y_i / L'(x_i) * 1/(x - x_i))，再通过NTT将这个有理函数转化为多项式，最终P(x) = L(x) * Q(x)

三、处理不支持NTT的有限域 #

如果你的素数p不满足NTT的条件（比如p-1没有足够大的2的幂次来覆盖n的长度），可以用这两个办法：

• Bluestein算法：实现任意长度的NTT，复杂度还是O(n logn)，只是常数会大一点
• 多模插值：先在几个满足NTT条件的素数域上计算出多项式系数，再用中国剩余定理合并到GF(p)上，适合p比较大的情况

实际实现小提示 #

• 2e6规模的NTT在C++里跑起来非常快，Python的话可以用优化过的库（比如自己实现快速卷积逻辑），但要注意内存和速度的平衡
• 所有模运算都要严格处理，避免负数结果（比如把负数加上p再取模）

希望这些方法能帮你解决问题！

备注：内容来源于stack exchange，提问作者Cnoob