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板题，所以我们有：
但是一般的题都是让你差一个值，那样的复杂度是(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;
}
但实际上可以(O(n^2))，
首先考虑到这里都是一个很长的多项式乘一次多项式(i.e.(x-i))所以多项式相乘其实是(O(n))的。
所以预处理(prod_{j=0}^{p-1} (x-j)),花费是(O(n^2))。
然后在(Sigma_1^n)的时候用(O(n))的多项式除法，这样也是(o(n^2))的。细节在下面。
下面一共有三种解法：
插值
1优化 O(P^2)
令(S(i,x) =prod_{j eq i}(x-j))
当(x e i)时,$ S(i,x)= frac{prod_{j=0}^{p-1} (x-j)}{x-i}$.
根据lagrange插值定理有
我们(O(n^{2}))预处理出(prod_{j=0}^{p-1} (x-j)) .
然后用多项式除法(O(n))算出每个(S(i,x)),阶乘公式(O(1))算出(S(i,i)).再一个$Sigma $ 求和。总的复杂度为(O(n^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优化 dp O(p^2)
由langrange得到
设(dp[i][j])为(prodlimits_{j=0}^{i}{(x-j)})的(x_j)的系数。这个dp 可逆，于是我们可以删掉一个元素得到(prodlimits_{j ot= 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})，
#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;
}
3优化 O(p^2)
注意到:
(x=1) 时
(x!=1) 时,根据费马小定理
所以
于是(S(i,x)=frac{prod_{j=0}^{p-1} (x-j)}{x-i}=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;
}
差分
5 O(p^2)
如果对某个等式(f_i equiv a_i pmod p) 计算第(i)阶差分，对所有的(j<i，f')中的(b_jx^j)不存在了。
可以考虑下面的公式。
我们计算(p-1)阶差分，(f^{(p-1)})只有(b_{p-1}x^{p-1})项留了下来， 于是我们得到(b_{p-1})。
将(b_{p-1}x^{p-1})减掉后继续进行p-2阶差分...and so on...
对于(Delta^if(x))，我们可以用差分公式求（(f_{t_0}, f_{t_1}, dots, f_{t_i})时x连续的等式）
同时我们也可以用这个公式(而非直接求导)求出(x^{n})的n阶差分，令(f(x)=x^n)即可
令(y=Delta^if(x),x=Delta^ix^{i})于是我们得到
由于我们已经减掉了(j>i)的(b_jx^j)这个差分只有(b_ix^i)的项。
我们可以证明(sum_{j=0}^{i} (-1)^jC_i^j(i-j)^i=i!) ((ecause (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 ? '
' : ' ');
}
return 0;
}
构造
5 O(p^2 log{p})
(O(p^2 log{p}))
注意到上面的题目都没用到(a_i={0,1})的性质。我们考虑构造。
若将(a_i=0)的 (i) 放到数组A中，(a_i=1)的(i)放到数组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)(G(x) + 1))符合了(f(A_i)=0), 但是(f(B_i)=F(x))。
如何将F(x)变为1？费马小定理即可
所以构造如下：
细节：
这样构造的问题是会出现大于(p-1)次项，这样就要再用一次费马小定理把次数(n>p-1)项的系数加到(n mod (p-1))的系数上.
这样就会出现第二个问题，首先p-1次项的系数始终为0，其次(n mod (p-1)=0) 时会把系数加到常数项上，这样显然是错的，比如(p=2)时(f(x)=x^2)不等于(f_0(x)=1)。
考虑到常数项必为(a_0)，所以我们一开始就把常数项(a_0)拿出来，构造出除了常数项外的(f^*(x)=f(x)-a_0),将这个多项式的常数项作为原多项式p-1次的系数。最后再把(a_0)放回去。
这样出现第三个问题，(a_0=1)时(f'(A_i)=-1,f'(B_i)=0).于是分类讨论，a=0时一样构造，a=1时f(x)于是(A,B)就要反过来，最后载将所有系数取个反。
#include <bits/stdc++.h>
using namespace std;
const int P = 3005;
int p, u;
bool neg;
vector<int> fi = {1}, se = {1};
vector<int> operator*(vector<int> &a, vector<int> &b)
{
int sz = min((int)a.size() + (int)b.size() - 1, p - 1);
vector<int> ans(sz, 0);
for (int i = 0; i < (int)a.size(); i++)
for (int j = 0; j < (int)b.size(); j++)
(ans[(i + j) % (p - 1)] += a[i] * b[j]) %= p;
return ans;
}
vector<int> operator^(vector<int> cur, int p)
{
if (p == 1)
return cur;
vector<int> tmp = cur ^ (p / 2);
tmp = tmp * tmp;
if (p & 1)
tmp = tmp * cur;
return tmp;
}
int main()
{
ios_base::sync_with_stdio(false);
cin.tie(nullptr);
cin >> p >> u;
neg = (u == 1);
for (int i = 1; i < p; i++)
{
cin >> u;
vector<int> cur = {p - i, 1};
if (u == neg)
fi = fi * cur;
else
se = se * cur;
}
(++se[0]) %= p;
fi = fi ^ (p - 1);
fi = fi * se;
if (neg)
for (int &v : fi)
v = (p - v) % p;
fi.push_back(fi[0]);
fi[0] = neg;
for (int &v : fi)
cout << v << " ";
}