AtCoder 137 F 插值 想法题

题意：AtCoder 137 F

已知：\(f(x)\)在p个点的值：\(f(i) \equiv a_i \pmod p\) \((0 \leq i \leq p-1)\)

求：\(f(x) = b_{p-1} x^{p-1} + b_{p-2} x^{p-2} + \ldots + b_0\)的所有系数\(b_i\)

保证：\(2 \leq p \leq 2999\) ，\(0 \leq a_i \leq 1\)

要求： \(0 \leq b_i \leq p-1\)

想法

来源 下面的n代指p，即多项式的最高次数+1

首先注意到这个题是lagrange板题，所以我们有：
\[ f(x)=\sum_{i=0}^{p-1}\left(\prod_{j\neq i}{\frac{x-x_j}{x_i-x_j}}\right)a_i \]
但是一般的题都是让你差一个值，那样的复杂度是\(O(n^2)\)(考虑到这题x连续可以优化到\(O(n)\))

而这个题让你把系数都求出来，于是你要用形如这样的\(O(n^4)\)板子。

Poly Lagrange(vector<type> x,vector<type> y){
    int n=x.size()-1;
    Poly ans;
    rep(k,0,n){
        Poly t,p;
        t=y[k];
        rep(j,0,n)if(j!=k){
            p=(vector<type>){frac(0,1)-x[j]/(x[k]-x[j]),frac(1,1)/(x[k]-x[j])};
            t=t*p;
        }
        ans=ans+t;
    }
    return ans;
}

插值

优化1 O(P^2)

令\(S(i,x) = \frac{\prod_{j=0}^{p-1} (x-j)}{x-i}\)得到
\[ f(x)=\sum_{i=0}^{p-1} \frac{\prod_{j=0}^{p-1} (x-j)}{x-i}*(\frac{a[i]}{S(j,j)})=\sum_{i=0}^{p-1}S(i,x)*\frac{a[i]}{S(i,i)} \]
注意到我们算\(f(t)\)时，若\(i \ne t\)，则\(S(i,t)*\frac{a[i]}{S(i,i)}\)总是为0.

所以得到\(S(t,t)*\frac{a[t]}{S(t,t)} = a[t]\)。 复杂度为O(P^2)

代码：

#include<bits/stdc++.h>
using namespace std;
const int maxn=3010;
int mod;
void plusle(int &a,int b){
    a+=b;if(a>=mod)a-=mod; return;
}
void minun(int &a,int b){
    a-=b; if(a<0)a+=mod; return;
}
int add(int a,int b) {
  a+=b;
  return a>=mod?a-mod:a;
}
int sub(int a,int b) {
  a-=b;
  return a<0?a+mod:a;
}
int mul(int a,int b) {
  return (int)((long long)a*b%mod);
}
int power(int a,int b) {
  int res=1;
  while (b>0) {
    if (b&1) {
      res=mul(res,a);
    }
    a=mul(a,a);
    b>>=1;
  }
  return res;
}
int inv(int x){
    return power(x,mod-2);
}
vector<int> mul(const vector<int> &p,int j){
    ///c[0]+(c[1]*x)+(c[2]*(x^2)+...)c[n-1]*(x^(n-1))
    ///(x-j)
    int n=(int)p.size();
    vector<int> s(n+1,0);
    for(int i=0;i<n;i++){
        plusle(s[i+1],p[i]);
    }
    for(int i=0;i<n;i++){
        plusle(s[i],mul(p[i],j));
    }
    return s;
}
vector<int> d( vector<int> p,int j){
    ///c[0]+c[1]*x+...+c[n-1]*(x^(n-1))
    ///x-j
    int n=(int)p.size();
    vector<int> s(n-1,0);
    for(int i=n-1;i>0;i--){
        plusle(p[i-1],mul(j,p[i]));
        s[i-1]=p[i];
    }
    return s;
}
int get(const vector<int> &g,int x){
    int res=0;
    int now=1;
    for(int i=0;i<(int)g.size();i++){
        plusle(res,mul(g[i],now));
        now=mul(now,x);
    }
    return res;
}
int a[maxn];

int main(){
    int p;
    scanf("%d",&p);
    mod=p;
    for(int i=0;i<p;i++){
        scanf("%d",&a[i]);
    }
    vector<int> t={0,1};
    for(int i=1;i<p;i++){
        t=mul(t,i);
    }
    vector<int> aa(p,0);
    for(int i=0;i<p;i++){
        if(a[i]==0)continue;
        vector<int> s=d(t,i);
        int cur=get(s,i);
        cur=inv(cur);
        for(int i=0;i<(int)s.size();i++){
            plusle(aa[i],mul(s[i],cur));
        }
    }
    for(int i:aa){
        printf("%d ",i);
    }
}
/*
    Good Luck
        -Lucina
*/

优化2 O(p^2)

注意到
\[ \sum_{i=1}^{p-1} x^i \equiv \begin{cases} -1, &x=1\\ 0. &x\ne 1 \end{cases} \]
于是\(S(i,x)=\sum_{j=1}^{p-1}-i^{p-j-1}x^j\)

#include <bits/stdc++.h>
using namespace std;

int r[3333];
int main() {
    int p;
    cin >> p;
    for (int i = 0; i < p; i++) {
        int a;
        cin >> a;
        if (!i) (r[0] += a) %= p;
        a = p - a;
        for (int d = 1; d < p; d++) {
            (r[p - d] += a) %= p;
            (a *= i) %= p;
        }
    }
    for (int i = 0; i < p; i++)
        cout << r[i] << ' ';
    cout << endl;
    return 0;
}

优化3 dp O(p^2)

由langrange得到
\[ y=\sum\limits_{i=0}^{p-1}a_i\prod\limits_{j \not= i}\frac{x-j}{i-j}\=\sum\limits_{i=0}^{p-1}a_i\prod\limits_{j \not= i}{(x-j)}\prod\limits_{j \not= i}\frac{1}{i-j} \]
设\(dp[i][j]\)为\(\prod\limits_{j=0}^{i}{(x-j)}\)的\(x_j\)的系数。这个dp 可逆，于是我们可以删掉一个元素得到\(\prod\limits_{j \not= i}{(x-j)}\)

\(dp[i][0]=dp[i-1][0]*i\)

\(dp[i][j]=dp[i-1][j]*i+dp[i-1][j-1]\)

实际上

\(dp[i][0]=\frac{dp[i+1][0]}{i+1}\)

\(dp[i][j]=\frac{dp[i+1][j]-dp[i][j-1]}{i+1}\)，所以我们可以删掉一个元素得到\(\prod\limits_{j=0}^{p-1}{(x-j)}\)

#include<bits/stdc++.h>
#define ld double
#define ull unsigned long long
#define ll long long
#define pii pair<int,int >
#define iiii pair<int,pii >
#define mp make_pair
#define INF 1000000000
#define MOD 1000000007
#define rep(i,x) for(int (i)=0;(i)<(x);(i)++)
inline int getint(){
    int x=0,p=1;char c=getchar();
    while (c<=32)c=getchar();
    if(c==45)p=-p,c=getchar();
    while (c>32)x=x*10+c-48,c=getchar();
    return x*p;
}
using namespace std;
//ruogu
const int N=3500;
int n,mod,a[N],inv[N],inv2[N],res[N];
int dp[N][N];
//
inline int add(int x,int y){x+=y;if(x>=mod)x-=mod;return x;}
inline int sub(int x,int y){x-=y;return (x<0)?x+mod:x;}
inline int mul(int x,int y){int ans=x*y;return ans%mod;}
inline int modpow(int x,int y){
    int ans=1;
    while(y){
        if(y&1)ans=mul(ans,x);
        x=mul(x,x);
        y>>=1;
    }
    return ans;
}
inline int modinv(int x){
    return modpow(x,mod-2);
}
int main(){
    n=getint();mod=n;
    rep(i,n)a[i]=getint();
    rep(i,N)inv[i]=modinv(i);
    rep(i,N)inv2[i]=modinv(sub(0,i));
    dp[n][0]=1;
    for(int i=n-1;i>=0;i--)for(int j=0;j<=n-i;j++){
        dp[i][j]=mul(dp[i+1][j],sub(0,i));
        if(j)dp[i][j]=add(dp[i][j],dp[i+1][j-1]);
        //cout<<i<<" "<<j<<" "<<dp[i][j]<<endl;
    }
    rep(i,n)if(a[i]){
        int ans=1;
        rep(j,i)ans=mul(ans,inv[i-j]);
        for(int j=i+1;j<n;j++)ans=mul(ans,inv2[j-i]);
        int tmp=0;
        if(!i)tmp=dp[1][0];
        res[0]=add(res[0],mul(tmp,ans));
        for(int j=1;j<n;j++){
            tmp=mul(sub(dp[0][j],tmp),inv2[i]);
            if(!i)tmp=dp[1][j];
            res[j]=add(res[j],mul(tmp,ans));
        }
    }
    rep(i,n)printf("%d ",res[i]);
    return 0;
}

构造

4 O(p^2 log{p})

常数较小的\(O(p^2 \log{p})\)

我们将\(a_i=0\)的项放到数组A中，\(a_i=1\)的项放到数组B中，构造两个多项式：

\(F(x) = (x - A_1)(x - A_2)...(x - A_k)\)

\(G(x) = (x - B_1)(x - B_2)...(x - B_{p - k})\)

则
\[ f(x)=F(x)^{p - 1}(G(x) + 1) \]

5 O(p^2)

对任意\(a_i\)取值的方法

计算第\(i\)阶差分，对所有的\(j<i，b_jx^j\)不存在了。

我们计算\(p-1\)阶差分，只留下了\(b_{p-1}x^{p-1}\)项， 于是得到\(b_{p-1}\)，然后将\(b_{p-1}x^{p-1}\)减掉。

继续进行p-2阶差分。。。

对于连续的方程\(f_{t_0}, f_{t_1}, \dots, f_{t_i}\)我们可以用公式\(\sum_{j=0}^{i} (-1)^jC_i^jf_{t_j}(x)\)。由于我们已经减掉了\(j>i\)的\(b_jx^j\)这个差分只有\(b_ix^i\)的项。

我们可以证明\(\sum_{j=0}^{i} (-1)^jC_i^j(i-j)^i=i!\) (\(\because (i!, p)=1\)),于是只有一个解。

#include <iostream>

using namespace std;

int bexp(int a, int x, int p) {
    if (x == 0) {
        return 1;
    }
    if (x % 2 == 1) {
        return a * bexp(a, x - 1, p) % p;
    }
    int t = bexp(a, x / 2, p);
    return t * t % p;
}

int inv(int a, int p) {
    return bexp(a, p - 2, p);
}

int main() {
    int p;
    cin >> p;
    int a[p];
    for (int i = 0; i < p; i++) {
        cin >> a[i];
    }

    int c[p][p];
    for (int i = 0; i < p; i++) {
        c[i][0] = c[i][i] = 1;
        for (int j = 1; j < i; j++) {
            c[i][j] = (c[i-1][j-1] + c[i-1][j]) % p;
        }
    }

    int e[p][p];
    for (int i = 0; i < p; i++) {
        for (int j = 0; j < p; j++) {
            e[i][j] = j == 0 ? 1 : e[i][j-1] * i % p;
        }
    }

    int b[p];
    for (int i = p - 1; i >= 0; i--) {
        int x = 0, y = 0;
        int neg = 1;
        for (int j = i; j >= 0; j--) {
            y = ((y + a[j] * c[i][j] * neg) % p + p) % p;
            x = ((x + e[j][i] * c[i][j] * neg) % p + p) % p;
            neg = -neg;
        }
        b[i] = x == 0 ? 0 : y * inv(x, p) % p;

        for (int j = 0; j < p; j++) {
            a[j] = ((a[j] - b[i] * e[j][i]) % p + p) % p;
        }
    }

    for (int i = 0; i < p; i++) {
        cout << b[i] << (i == p - 1 ? '\n' : ' ');
    }

    return 0;
}

AtCoder 137 F 插值 想法题

原文：https://www.cnblogs.com/SuuT/p/11349024.html