引言
在数值计算和算法竞赛中，我们经常遇到这样的问题：已知一个多项式在某些点的取值，如何快速求出该多项式在任意点的值？或者反过来，如何根据若干离散数据点构造出通过这些点的多项式？拉格朗日插值（Lagrange Interpolation）正是解决这类问题的经典方法。它以简洁优美的公式，将插值问题转化为基函数的线性组合。本文将用 Java 完整实现拉格朗日插值，并展示其在数列通项公式求解、前缀和快速计算等场景中的应用。
插值问题定义
给定 n+1 个互不相同的点 (x₀, y₀), (x₁, y₁), …, (xₙ, yₙ)，其中 yᵢ = f(xᵢ)。拉格朗日插值的目标是找到一个次数不超过 n 的多项式 P(x)，使得 P(xᵢ) = yᵢ 对所有 i 成立。
唯一性定理
满足上述条件的插值多项式是唯一的。这意味着无论使用何种方法（拉格朗日、牛顿、范德蒙德），得到的插值多项式在代数上是等价的。
拉格朗日插值公式
基函数构造
拉格朗日插值的核心思想是构造一组基函数 Lᵢ(x)，每个基函数满足：
Lᵢ(xⱼ) = 1 当 i = j
Lᵢ(xⱼ) = 0 当 i ≠ j
这样的基函数可以构造为：
(x - x₀)(x - x₁)...(x - xᵢ₋₁)(x - xᵢ₊₁)...(x - xₙ)
Lᵢ(x) = ───────────────────────────────────────────────────
(xᵢ - x₀)(xᵢ - x₁)...(xᵢ - xᵢ₋₁)(xᵢ - xᵢ₊₁)...(xᵢ - xₙ)
即分子是所有 (x – xⱼ) 的乘积（j ≠ i），分母是将 x 替换为 xᵢ 后的值。
插值多项式
最终的插值多项式为这些基函数的线性组合：
P(x) = y₀·L₀(x) + y₁·L₁(x) + ... + yₙ·Lₙ(x)
算法核心思想
O(n²) 插值法
直接按照公式计算，每次求值需要 O(n²) 时间。
O(n) 快速求值（连续整数点优化）
当插值点 x₀, x₁, …, xₙ 为连续整数（如 0, 1, 2, …, n）时，可以进行大量预计算，将单次求值优化到 O(n)：
- 预计算前缀积和后缀积
- 预计算阶乘和阶乘逆元
- 利用阶乘表示分母，避免重复计算
对于连续整数点 0, 1, …, n，分母可简化为：
(xᵢ - x₀)...(xᵢ - xᵢ₋₁)(xᵢ - xᵢ₊₁)...(xᵢ - xₙ)
= i! · (-1)^(n-i) · (n-i)!
Java 完整实现
import java.util.*;
/**
* 拉格朗日插值（Lagrange Interpolation）
* 支持一般点和连续整数点的多项式插值
*/
public class LagrangeInterpolation {
private static final long MOD = 998244353L; // 常用模数（质数）
private long[] x; // 插值点横坐标
private long[] y; // 插值点纵坐标
private int n; // 点的数量（多项式次数为 n-1）
public LagrangeInterpolation(long[] x, long[] y) {
if (x.length != y.length) {
throw new IllegalArgumentException("x and y must have same length");
}
this.n = x.length;
this.x = Arrays.copyOf(x, n);
this.y = Arrays.copyOf(y, n);
}
/**
* 快速幂：计算 a^b mod MOD
*/
private static long power(long a, long b) {
long result = 1;
a %= MOD;
while (b > 0) {
if ((b & 1) == 1) {
result = (result * a) % MOD;
}
a = (a * a) % MOD;
b >>= 1;
}
return result;
}
/**
* 模逆元：计算 a^(-1) mod MOD
*/
private static long modInverse(long a) {
return power(a, MOD - 2);
}
/**
* 通用拉格朗日插值：O(n²) 求值
* 适用于任意插值点
* @param targetX 要求值的点
* @return P(targetX)
*/
public double interpolate(double targetX) {
double result = 0.0;
for (int i = 0; i < n; i++) {
double li = 1.0;
// 计算基函数 L_i(targetX)
for (int j = 0; j < n; j++) {
if (i != j) {
li *= (targetX - x[j]) / (x[i] - x[j]);
}
}
result += y[i] * li;
}
return result;
}
/**
* 模意义下的通用拉格朗日插值
* @param targetX 要求值的点
* @return P(targetX) mod MOD
*/
public long interpolateMod(long targetX) {
long result = 0;
for (int i = 0; i < n; i++) {
long numerator = 1;
long denominator = 1;
for (int j = 0; j < n; j++) {
if (i != j) {
numerator = (numerator * ((targetX - x[j] + MOD) % MOD)) % MOD;
denominator = (denominator * ((x[i] - x[j] + MOD) % MOD)) % MOD;
}
}
long li = (numerator * modInverse(denominator)) % MOD;
result = (result + (y[i] % MOD + MOD) % MOD * li) % MOD;
}
return result;
}
/**
* 连续整数点（0, 1, ..., n-1）的 O(n) 快速求值
* 需要预计算阶乘
* @param y 在 0, 1, ..., n-1 处的函数值
* @param targetX 要求值的点
* @return P(targetX)
*/
public static long interpolateConsecutive(long[] y, long targetX) {
int n = y.length;
// 预计算阶乘
long[] fact = new long[n + 1];
long[] invFact = new long[n + 1];
fact[0] = 1;
for (int i = 1; i <= n; i++) {
fact[i] = (fact[i - 1] * i) % MOD;
}
invFact[n] = modInverse(fact[n]);
for (int i = n - 1; i >= 0; i--) {
invFact[i] = (invFact[i + 1] * (i + 1)) % MOD;
}
// 预计算前缀积和后缀积
long[] pre = new long[n];
long[] suf = new long[n];
// pre[i] = (targetX - 0) * (targetX - 1) * ... * (targetX - (i-1))
pre[0] = 1;
for (int i = 1; i < n; i++) {
pre[i] = (pre[i - 1] * ((targetX - (i - 1) + MOD) % MOD)) % MOD;
}
// suf[i] = (targetX - (n-1)) * (targetX - (n-2)) * ... * (targetX - (i+1))
suf[n - 1] = 1;
for (int i = n - 2; i >= 0; i--) {
suf[i] = (suf[i + 1] * ((targetX - (i + 1) + MOD) % MOD)) % MOD;
}
long result = 0;
for (int i = 0; i < n; i++) {
// 分子 = pre[i] * suf[i]
long numerator = (pre[i] * suf[i]) % MOD;
// 分母 = i! * (-1)^(n-1-i) * (n-1-i)!
long sign = ((n - 1 - i) % 2 == 0) ? 1 : MOD - 1;
long denominator = (fact[i] * fact[n - 1 - i]) % MOD;
denominator = (denominator * sign) % MOD;
long li = (numerator * modInverse(denominator)) % MOD;
result = (result + ((y[i] % MOD + MOD) % MOD) * li) % MOD;
}
return result;
}
/**
* 求数列的通项公式值（已知前几项，求第k项）
* 假设数列由一个多项式生成
* @param firstTerms 数列的前几项
* @param k 要求的项数（从0开始计数）
* @return 第k项的值
*/
public static long findKthTerm(long[] firstTerms, long k) {
return interpolateConsecutive(firstTerms, k);
}
/**
* 多项式前缀和快速计算
* 若 f(k) = 1^k + 2^k + ... + n^k 是一个关于 n 的 k+1 次多项式
* 已知 f(0), f(1), ..., f(k+1)，求 f(n)
* @param prefixValues f(0), f(1), ..., f(k+1)
* @param n 要求的前缀和上限
* @return f(n)
*/
public static long polynomialPrefixSum(long[] prefixValues, long n) {
return interpolateConsecutive(prefixValues, n);
}
/**
* 打印插值多项式的基函数展开（用于展示）
*/
public void printBasisFunctions() {
System.out.println("=== Lagrange Basis Functions ===");
for (int i = 0; i < n; i++) {
System.out.print("L_" + i + "(x) = ");
boolean first = true;
for (int j = 0; j < n; j++) {
if (i != j) {
if (!first) System.out.print(" * ");
System.out.print("(x - " + x[j] + ")/(" + x[i] + " - " + x[j] + ")");
first = false;
}
}
System.out.println();
}
}
// ==================== 应用示例 ====================
/**
* 示例1：通用插值（非连续点）
*/
public static void demoGeneralInterpolation() {
System.out.println("=== Demo 1: General Interpolation ===");
// 已知 f(1)=2, f(2)=3, f(4)=5，求 f(3)
long[] x = {1, 2, 4};
long[] y = {2, 3, 5};
LagrangeInterpolation li = new LagrangeInterpolation(x, y);
li.printBasisFunctions();
double result = li.interpolate(3);
System.out.println("P(3) = " + result); // 应接近 4
// 验证已知点
for (int i = 0; i < x.length; i++) {
System.out.println("P(" + x[i] + ") = " + li.interpolate(x[i]) + " (expected: " + y[i] + ")");
}
}
/**
* 示例2：连续整数点的 O(n) 快速求值
*/
public static void demoConsecutiveInterpolation() {
System.out.println("\n=== Demo 2: Consecutive Points (O(n)) ===");
// 已知 f(0)=0, f(1)=1, f(2)=4, f(3)=9，求 f(100)
// 这对应于 f(x) = x^2
long[] y = {0, 1, 4, 9};
long result = LagrangeInterpolation.interpolateConsecutive(y, 100);
System.out.println("f(100) = " + result + " (expected: 10000)");
// 验证中间点
for (int i = 0; i < y.length; i++) {
long val = LagrangeInterpolation.interpolateConsecutive(y, i);
System.out.println("f(" + i + ") = " + val + " (expected: " + y[i] + ")");
}
}
/**
* 示例3：求数列通项公式
*/
public static void demoSequenceTerm() {
System.out.println("\n=== Demo 3: Sequence General Term ===");
// 数列：0, 1, 5, 14, 30... （三角锥数 / 四面体数）
// 通项公式：n(n+1)(n+2)/6
long[] terms = {0, 1, 5, 14, 30};
System.out.println("Sequence: 0, 1, 5, 14, 30...");
long term10 = LagrangeInterpolation.findKthTerm(terms, 10);
System.out.println("10th term = " + term10 + " (expected: 220)");
long term20 = LagrangeInterpolation.findKthTerm(terms, 20);
System.out.println("20th term = " + term20 + " (expected: 1540)");
}
/**
* 示例4：多项式前缀和
* 求 1^2 + 2^2 + ... + n^2 = n(n+1)(2n+1)/6
*/
public static void demoPrefixSum() {
System.out.println("\n=== Demo 4: Polynomial Prefix Sum ===");
// f(n) = 1^2 + 2^2 + ... + n^2，这是关于 n 的3次多项式
// f(0)=0, f(1)=1, f(2)=5, f(3)=14
long[] prefixValues = {0, 1, 5, 14};
long sum100 = LagrangeInterpolation.polynomialPrefixSum(prefixValues, 100);
System.out.println("1^2 + 2^2 + ... + 100^2 = " + sum100 + " (expected: 338350)");
long sum1000 = LagrangeInterpolation.polynomialPrefixSum(prefixValues, 1000);
System.out.println("1^2 + 2^2 + ... + 1000^2 = " + sum1000 + " (expected: 333833500)");
}
/**
* 示例5：三次多项式拟合
*/
public static void demoCubicPolynomial() {
System.out.println("\n=== Demo 5: Cubic Polynomial Fitting ===");
// f(x) = 2x^3 - 3x^2 + x - 1
// f(0)=-1, f(1)=-1, f(2)=5, f(3)=29
long[] y = {-1, -1, 5, 29};
long[] testPoints = {4, 5, 10};
for (long x : testPoints) {
long expected = 2*x*x*x - 3*x*x + x - 1;
long result = LagrangeInterpolation.interpolateConsecutive(y, x);
System.out.println("f(" + x + ") = " + result + " (expected: " + expected + ")");
}
}
public static void main(String[] args) {
demoGeneralInterpolation();
demoConsecutiveInterpolation();
demoSequenceTerm();
demoPrefixSum();
demoCubicPolynomial();
}
}
算法复杂度分析
| 方法 | 预处理 | 单次求值 | 适用场景 |
|---|---|---|---|
| 通用 O(n²) 插值 | O(1) | O(n²) | 任意插值点 |
| 连续整数 O(n) 优化 | O(n) | O(n) | xᵢ = 0, 1, …, n-1 |
| 牛顿插值 | O(n²) | O(n) | 需要多次求值 |
| 范德蒙德矩阵 | O(n³) | O(n²) | 理论分析 |
对于连续整数点的 O(n) 优化：
– 预计算阶乘和阶乘逆元：O(n)
– 预计算前缀积和后缀积：O(n)
– 单次求值：O(n)
– 总空间复杂度：O(n)
拉格朗日插值 vs 其他插值方法
| 方法 | 优点 | 缺点 |
|---|---|---|
| 拉格朗日 | 公式直观、易于理解 | 添加新点需重新计算 |
| 牛顿 | 添加新点只需 O(n) | 公式不够直观 |
| 分段线性 | 简单、避免龙格现象 | 光滑性差 |
| 样条插值 | 光滑性好 | 计算复杂 |
拉格朗日插值最适合于：插值点固定、需要理论分析或单次/少量求值的场景。
应用场景
拉格朗日插值在以下场景中非常实用：
- 数列通项求解：已知多项式型数列的前几项，快速求任意项
- 前缀和快速计算：将求和转化为多项式求值，如自然数幂和
- 组合数恒等式证明：通过验证有限个点成立来证明多项式恒等
- 有限域计算：在模质数的域中进行插值，避免浮点精度问题
- 信号处理：离散数据的连续化重建
注意事项
- 龙格现象：等距节点的高次插值可能在区间端点附近产生剧烈振荡
- 数值稳定性：大量乘法运算可能导致浮点误差累积，模运算场景更稳定
- 唯一性保证：n+1 个点确定唯一的 n 次插值多项式
- 计算优化：连续整数点是竞赛中最常见的情况，务必掌握 O(n) 优化
总结
拉格朗日插值以优雅的对称公式，将多项式插值问题转化为基函数的线性组合。对于连续整数点的特殊场景，通过阶乘预计算和前后缀积优化，可以将求值复杂度降至 O(n)。掌握拉格朗日插值，将帮助你在数列求通项、前缀和优化等问题中快速找到高效解法。
