您的位置:首页 > 其它

[多项式求逆] 51Nod TalkingData数据科学精英夏令营挑战赛 F 驴蛋蛋与老孙与微分式

2017-06-18 14:09 281 查看
题解里给出了一种解偏微分方程得出生成函数的方法

最后应该是

H(z,x)=sinz+xcoszcosz−xsinz

H(z)=sinzcosz=tanz=x+13x3+215x5+o(x5)

实际上我们打一个表 把它喂给OEIS 然后就得到了A009006

如果令F(0)=1,他确实是
E.g.f. 1 + tan(x).


得到递推公式b(0)=1,b(n)=−∑n−1k=0b(k)∗(nk)∗2n−k−1,a(n)=|b(n)|

无论怎样最后都是多项式求逆

一个多项式题 模数是109+7就不好了

#include<cstdio>
#include<cstdlib>
#include<algorithm>
using namespace std;
typedef long long ll;

const int P=1e9+7;
const int INV2=(P+1)/2;
const int M[]={998244353,1004535809,469762049};
const int G[]={3,3,3};
const ll _M=(ll)M[0]*M[1];

inline ll Pow(ll a,int b,int p){
ll ret=1;
for (;b;b>>=1,a=a*a%p) if (b&1) ret=ret*a%p;
return ret;
}
inline ll mul(ll a,ll b,ll p){
a%=p; b%=p;
return ((a*b-(ll)((ll)((long double)a/p*b+1e-3)*p))%p+p)%p;
}

const int m1=M[0],m2=M[1],m3=M[2];
const int inv1=Pow(m1%m2,m2-2,m2),inv2=Pow(m2%m1,m1-2,m1),inv12=Pow(_M%m3,m3-2,m3);
inline int CRT(int a1,int a2,int a3){
ll A=(mul((ll)a1*m2%_M,inv2,_M)+mul((ll)a2*m1%_M,inv1,_M))%_M;
ll k=((ll)a3+m3-A%m3)*inv12%m3;
return (k*(_M%P)+A)%P;
}

const int N=550000;
int R
;
struct NTT{
int P,G;
int num,w[2]
;
void Pre(int _P,int _G,int m){
num=m; P=_P; G=_G;
int g=Pow(G,(P-1)/num,P);
w[1][0]=1; for (int i=1;i<num;i++) w[1][i]=(ll)w[1][i-1]*g%P;
w[0][0]=1; for (int i=1;i<num;i++) w[0][i]=w[1][num-i];
}
void FFT(int *a,int n,int r){
for (int i=0;i<n;i++) if (i<R[i]) swap(a[i],a[R[i]]);
for (int i=1;i<n;i<<=1)
for (int j=0;j<n;j+=(i<<1))
for (int k=0;k<i;k++){
int x=a[j+k],y=(ll)a[j+i+k]*w[r][num/(i<<1)*k]%P;
a[j+k]=(x+y)%P; a[j+i+k]=(x+P-y)%P;
}
if (!r) for (int i=0,inv=Pow(n,P-2,P);i<n;i++) a[i]=(ll)a[i]*inv%P;
}
}ntt[3];

int tmp
,b2
,b3
,_b[3]
,c
;

inline void GetInv(int *a,int *b,int n){
if (n==1) return void(b[0]=CRT(Pow(a[0],m1,m1-2),Pow(a[0],m2,m2-2),Pow(a[0],m3,m3-2)));
GetInv(a,b,n>>1);
int L=0; while (!(n>>L&1)) L++;
for (int i=1;i<(n<<1);i++) R[i]=(R[i>>1]>>1)|((i&1)<<L);

for (int i=0;i<n;i++) b2[i]=b3[i]=b[i],tmp[i]=a[i],b2[i+n]=b3[i+n]=0;
ntt[0].FFT(b,n<<1,1);
ntt[1].FFT(b2,n<<1,1);
ntt[2].FFT(b3,n<<1,1);

for (int I=0;I<3;I++){
int *d=I==0?b:(I==1?b2:b3);
for (int i=0;i<n;i++) tmp[n+i]=0,tmp[i]=a[i];
ntt[I].FFT(tmp,n<<1,1);
for (int i=0;i<(n<<1);i++)
tmp[i]=(ll)tmp[i]*d[i]%M[I];
ntt[I].FFT(tmp,n<<1,0);
for (int i=0;i<(n<<1);i++) _b[I][i]=tmp[i];
}

for (int i=0;i<(n<<1);i++) c[i]=CRT(_b[0][i],_b[1][i],_b[2][i]),c[i]=c[i]==0?0:P-c[i];
(c[0]+=2)%=P;

for (int I=0;I<3;I++){
int *d=I==0?b:(I==1?b2:b3);
for (int i=0;i<n;i++) tmp[n+i]=0,tmp[i]=c[i];
ntt[I].FFT(tmp,n<<1,1);
for (int i=0;i<(n<<1);i++)
tmp[i]=(ll)tmp[i]*d[i]%M[I];
ntt[I].FFT(tmp,n<<1,0);
for (int i=0;i<(n<<1);i++) _b[I][i]=tmp[i];
}

for (int i=0;i<n;i++) b[i]=CRT(_b[0][i],_b[1][i],_b[2][i]),b[n+i]=0;
}

ll fac
,inv
;
ll pw
;

inline void Init(int n){
fac[0]=1; for (int i=1;i<=n;i++) fac[i]=fac[i-1]*i%P;
inv[1]=1; for (int i=2;i<=n;i++) inv[i]=(ll)(P-P/i)*inv[P%i]%P;
inv[0]=1; for (int i=1;i<=n;i++) inv[i]=inv[i]*inv[i-1]%P;
pw[0]=1; for (int i=1;i<=n;i++) pw[i]=pw[i-1]*2%P;
}
inline ll C(int n,int m){
return fac
*inv[m]%P*inv[n-m]%P;
}

int n,m;
int a
,b
;
int _a[3]
,B
;

inline int Solve(ll n,int m){
static ll pw
;
pw[0]=1; for (int i=1;i<=m+1;i++) pw[i]=pw[i-1]*(n%P)%P;
ll ret=((ll)m+1)*((P+1)>>1)%P*pw[m]%P;
for (int k=0;k<=m;k+=2)
ret+=C(m+1,k)*B[k]%P*pw[m+1-k]%P;
return ret%P*Pow(m+1,P-2,P)%P;
}

int main(){
n=100000;
for (m=1;m<=n;m<<=1);
for (int i=0;i<3;i++) ntt[i].Pre(M[i],G[i],m<<1);
Init(m+1);
for (int i=1;i<m;i++) a[i]=inv[i]*INV2%P;
a[0]=1;
GetInv(a,b,m);
for (int i=1;i<=n;i++) B[i]=(ll)b[i]*fac[i]%P*pw[i]%P;
for (int i=1;i<=n;i+=4) B[i]=(P-B[i])%P;
for (int i=1;i<=n;i++) (B[i]+=B[i-1])%=P;
int Q,_n;
scanf("%d",&Q);
while (Q--){
scanf("%d",&_n);
printf("%d\n",B[_n]);
}
return 0;
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: