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

题解里给出了一种解偏微分方程得出生成函数的方法
最后应该是
H(z,x)=sinz+xcoszcoszxsinz
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)=n1k=0b(k)(nk)2nk1,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[N];
struct NTT{
  int P,G;
  int num,w[2][N];
  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[N],b2[N],b3[N],_b[3][N],c[N];

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[N],inv[N];
ll pw[N];

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[n]*inv[m]%P*inv[n-m]%P;
}

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

inline int Solve(ll n,int m){
  static ll pw[N];
  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;
}
  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值