拉格朗日插值和求多项式系数

拉格朗日介绍

先说说拉格朗日是啥吧
首先 拉格朗日插值是给你 n+1 个点 $(x,y)$ 然后根据这n个点可以$O(n^2)$的求出多项式的系数。也就是解出这个多项式的答案。

假设给你一个多项式
$y=a_0+a_1x+a_2x^2$
然后给你3个解$(x1,y1)(x2,y2)(x3,y3)$你第一个想法是怎么解？解方程啊是不是
代进去是不是这样

$$\begin{cases} y_1=a_0+a_1x_1+a_2x_1^2\\ y_2=a_0+a_1x_2+a_2x_2^2\\ y_3=a_0+a_1x_3+a_2x_3^2 \end{cases}$$

然年后解这个方程？
解这个方程复杂度多少,高斯消元O(n^3)很显然复杂度高了。
拉格朗日就比较厉害了他能O(n^2)解决
首先 假设一个多项式$f_1(x)= b_0 + b_1x+b_2x^2$
当他$x=x_1$解是1,$x_2$ $x_3$ 解是 0
同理再假设 $f_2(x)$ $f_3(x)$
然后$L(x)=y_1f_1(x)+y_2f_2(x)+y_3f_3(x)$,这个就是最开始那个方程，不信?你分别把$x_1,x_2,x_3$ 带进去解绝对是 $y_1,y_2,y_3$。
那么问题来了后面$f_1(x)$这个多项式怎么求出来？？？？
这就是拉格朗日基本公式

再来一遍，这就是拉格朗日基本公式

$$f_j(x)=l_j(x)=\prod_{j=0,j\neq i}^{n-1}\frac{x-x_j}{x_j-x_i}$$

把这个多项式展开会发现非常神奇的事，当$x=x_j$的时候刚好等于1否则等于0,刚好满足了原来所需要的方程式。就是下面这样：

$$f_j(x)=l_j(x)=\prod_{j=0,j\neq i}^{k}\frac{x-x_j}{x_i-x_j}=\frac{x-x_0}{x_0-x_j}\cdots\frac{x-x_0}{x_j-x_{j-1}}\frac{x-x_{j+1}}{x_j-x_{j+1}}\cdots\frac{x-x_k}{x_j-x_k}$$

最终表达式

$$L(x)=\sum_{i=0}^{n-1}y_i\prod_{j=0,j\neq i}^{n-1}\frac{x-x_j}{x_j-x_i}$$

怕你还是看不懂，举个例子给你看

$f(4)=10\ f(5)=5.25\ f(6)=1$ 求 $f(18)$

首先写出拉格朗日基本多项式

$$l_0=\frac{(x-5)(x-6)}{(4-5)(5-6)}$$ $$l_1=\frac{(x-4)(x-6)}{(5-4)(5-6)}$$ $$l_2=\frac{(x-4)(x-5)}{(6-4)(6-5)}$$ $$\begin{cases} p(x)&=f(4)l_0(x)+f(5)l_1(x)+f(6)l_2(x)\\ &=10\frac{(x-5)(x-6)}{(4-5)(5-6)}+5.25\frac{(x-4)(x-6)}{(5-4)(5-6)}+1*\frac{(x-4)(x-5)}{(6-4)(6-5)}\\ &=\frac{1}{4}(x^2-28x+136) \end{cases}$$

此时代入18:$f(18)=p(18)=-11$
还有一个问题你们肯定很想问。。。知道公式之后怎么解。。。
对于这个问题
分母 是不是每次算一下就行了，答案是固定的
分子是不是一个大的多项式里面少了一个，就预处理出总的多项式然后，模拟除一下$(x+c)$的多项式

理论知识全部搞定，下面就给你贴模板了

板子

//这个是杜教的板子 我打了点注释

/// 注意mod，使用前须调用一次 polysum::init(int M);
/// 注意mod，使用前须调用一次 polysum::init(int M);
namespace polysum {
#define rep(i,a,n) for (int i=a;i<n;i++)
#define per(i,a,n) for (int i=n-1;i>=a;i--)
typedef long long ll;
const ll mod = 1e9 + 7; /// 取模值
ll powmod(ll a, ll b) {
    ll res = 1;
    a %= mod;
    assert(b >= 0);
    for (; b; b >>= 1) {
        if (b & 1)res = res*a%mod;
        a = a*a%mod;
    }
    return res;
}

const int D = 1010000; /// 最高次限制
ll a[D], f[D], g[D], p[D], p1[D], p2[D], b[D], h[D][2], C[D];
ll calcn(int d, ll *a, ll n) {          //根据前 d 项 求 第n项
    if (n <= d) return a[n];
    p1[0] = p2[0] = 1;
    rep(i, 0, d + 1) {
        ll t = (n - i + mod) % mod;
        p1[i + 1] = p1[i] * t%mod;
    }
    rep(i, 0, d + 1) {
        ll t = (n - d + i + mod) % mod;
        p2[i + 1] = p2[i] * t%mod;
    }
    ll ans = 0;
    rep(i, 0, d + 1) {
        ll t = g[i] * g[d - i] % mod*p1[i] % mod*p2[d - i] % mod*a[i] % mod;
        if ((d - i) & 1) ans = (ans - t + mod) % mod;
        else ans = (ans + t) % mod;
    }
    return ans;
}
void init(int M) { /// M：最高次
    f[0] = f[1] = g[0] = g[1] = 1;
    rep(i, 2, M + 5) f[i] = f[i - 1] * i%mod;
    g[M + 4] = powmod(f[M + 4], mod - 2);
    per(i, 1, M + 4) g[i] = g[i + 1] * (i + 1) % mod;   //逆元
}

ll polysum(ll n, ll *arr, ll m) { /// a[0].. a[m] \sum_{i=0}^{n-1} a[i]
    for (int i = 0; i <= m; i++)
        a[i] = arr[i];
    a[m + 1] = calcn(m, a, m + 1);
    rep(i, 1, m + 2) a[i] = (a[i - 1] + a[i]) % mod;
    return calcn(m + 1, a, n - 1);
}
ll qpolysum(ll R, ll n, ll *a, ll m) { /// a[0].. a[m] \sum_{i=0}^{n-1} a[i]*R^i
    if (R == 1) return polysum(n, a, m);
    a[m + 1] = calcn(m, a, m + 1);
    ll r = powmod(R, mod - 2), p3 = 0, p4 = 0, c, ans;
    h[0][0] = 0;
    h[0][1] = 1;
    rep(i, 1, m + 2) {
        h[i][0] = (h[i - 1][0] + a[i - 1])*r%mod;
        h[i][1] = h[i - 1][1] * r%mod;
    }
    rep(i, 0, m + 2) {
        ll t = g[i] * g[m + 1 - i] % mod;
        if (i & 1) p3 = ((p3 - h[i][0] * t) % mod + mod) % mod, p4 = ((p4 - h[i][1] * t) % mod + mod) % mod;
        else p3 = (p3 + h[i][0] * t) % mod, p4 = (p4 + h[i][1] * t) % mod;
    }
    c = powmod(p4, mod - 2)*(mod - p3) % mod;
    rep(i, 0, m + 2) h[i][0] = (h[i][0] + h[i][1] * c) % mod;
    rep(i, 0, m + 2) C[i] = h[i][0];
    ans = (calcn(m, C, n)*powmod(R, n) - c) % mod;
    if (ans<0) ans += mod;
    return ans;
}
}

然后下面这个是求多项式的

LL temp[maxn];

void mul(LL *f, int len, LL t) { //len为多项式的次数+1，函数让多项式f变成f*(x+t)
    for (int i = len; i > 0; --i) {
        temp[i] = f[i];
        f[i] = f[i - 1];
    }
    temp[0] = f[0], f[0] = 0;
    for (int i = 0; i <= len; ++i) {
        f[i] = (f[i] + t * temp[i]) % mod;
    }
}

void dev(LL *f, LL *r, LL t, int len) { //f是被除多项式的系数，r保存f除以x+t的结果 len是最高次项
    for (int i = 0; i <= len; ++i) {
        temp[i] = f[i];
    }
    for (int i = len; i > 0; --i) {
        r[i - 1] = temp[i];
        temp[i - 1] = (temp[i - 1] - t * temp[i]) % mod;

    }
    return;
}

LL a[maxn], b[maxn], c[maxn];
LL x[maxn], y[maxn]; //x,y输入从 1开始到n
int n;

void lglr() {
    memset(a, 0, sizeof a);
    b[1] = 1, b[0] = -x[1];
    for (int i = 2; i <= n; ++i) {
        mul(b, i, -x[i]);
    }//预处理(x-x1)*(x-x2)...*(x-xn)
    for (int i = 1; i <= n; ++i) {
        LL fz = 1;
        for (int j = 1; j <= n; ++j) {
            if (j == i) continue;
            fz = fz * (x[i] - x[j]) % mod;
        }
        fz = qm(fz, mod - 2);
        fz = fz * y[i] % mod;//得到多项式系数
        dev(b, c, -x[i], n);//得到多项式，保存在b数组
        for (int j = 0; j < n; ++j) a[j] = (a[j] + fz * c[j]) % mod;
    }
}

LL cal(LL k) {  //计算第x=k值
    LL ans = 0;
    LL res = 1;
    for (int i = 0; i < n; ++i) {
        ans = (ans + res * a[i]) % mod;
        res = res * k % mod;
    }
    ans = (ans + mod) % mod;
    return ans;
}