AtCoder 137 F 插值 差分 或构造
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\neq i}(x-j)\)
当\(x\ne 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]\)为\(\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}\)，
#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!\) (\(\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;
}
构造
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 << " ";
}
