模板汇总3.3_数学3

1.快速傅里叶变换(FFT)

①FFT用来干啥

加速向量卷积和多项式乘法(指在OI中)

②FFT原理

因为写的时候时间有点紧，这里默认前置知识都掌握了

考虑多项式的点值表达，这时两个多项式可以做到$O(n)$相乘

但是平时多项式都是用系数给出的，所以我们实际上要这样操作

只要①②(其中①叫做多项式求值，②叫做多项式插值)中较差的复杂度不超过$O(n^2)$即可相对快速完成计算

这里我们把①叫做DFT(离散傅里叶变换)，②叫做IDFT(逆DFT)，因为两者实现类似这里只说DFT。加速的方法是用单位根计算对应的点值

把单位根一个个带进去DFT显然还是$O(n^2)$的，但是傅里叶fa♂现了一个很强的东西(别问我怎么发现的，要不然为什么叫FFT啊=。=)：我们将单位根带入一个DFT F后按$k$($k<\frac{n}{2}$)奇偶分类得到两坨DFT F1和F2，假如F1F2都求出来了，那么这时我们可以$O(n)$合并出F，于是我们一路分治最后合起来$O(n\log n)$即可解决问题

那么如何$O(n)$合并出F？这就要依靠单位根的性质：$ω_{dn}^{dk}=ω_n^k(d>0)$，这样我们奇偶分完之后从奇数那坨里提出来一个$ω_n^k$就会得到我们想要的东西

$DFT(ω_n^k)=(even)DFT(ω_{\frac{n}{2}}^{k})+(odd)ω_n^kDFT(ω_{\frac{n}{2}}^{k})$

$DFT(ω_n^{k+\frac{n}{2}})=(even)DFT(ω_{\frac{n}{2}}^{k})-(odd)ω_n^kDFT(ω_{\frac{n}{2}}^{k})$

接下来是一些具体实现的问题

Q：万一中间分治的时候两边分不均怎么办

A：先用零补成2的整次幂再做

Q：直接做(指递归)怎么过不去啊，你是不是讲的假的啊

A：那是你写萎了，我都能过其实不一定要递归，我们观察最后递归结束时原多项式被划分成了什么样子，然后发现每个位置都跑到了它二进制翻转之后的那个位置，于是预处理一发然后就可以手动分位置了，省掉了递归的常数和数组空间，可以快不少

Q：说说非算法部分的一些东西的实现吧(什么鬼

A：首先你要写一个复数struct；然后三角函数很慢，完全可以预处理，单位根们**乘的时候注意下顺序，好像没啥了......

③FFT复杂度

时间复杂度：$O(n\log n)$

④FFT具体实现

递归版(总共2000+ms，最慢的点700+ms，比较慢，但是强烈建议先写递归版把算法搞懂)

 1 #include<cmath>
 2 #include<cstdio>
 3 #include<cctype>
 4 #include<cstring>
 5 #include<algorithm>
 6 using namespace std;
 7 const int N=4e6+6,M=30;
 8 const double pai=acos(-1);
 9 struct cpx
10 {
11     double x,y;
12 }a[N],b[N];
13 cpx operator + (cpx a,cpx b)
14 {
15     return (cpx){a.x+b.x,a.y+b.y};
16 } 
17 cpx operator - (cpx a,cpx b)
18 {
19     return (cpx){a.x-b.x,a.y-b.y};
20 }
21 cpx operator * (cpx a,cpx b)
22 {
23     double x1=a.x,x2=b.x,y1=a.y,y2=b.y;
24     return (cpx){x1*x2-y1*y2,x1*y2+x2*y1};
25 }
26 int n,m;
27 double Sin[M],Cos[M];
28 void Read(double &x)
29 {
30     int ret=0; x=0;
31     char ch=getchar();
32     while(!isdigit(ch))
33         ch=getchar();
34     while(isdigit(ch))
35         ret=(ret<<1)+(ret<<3)+(ch^48),ch=getchar();
36     x=ret; return ;
37 }
38 void Write(int x)
39 {
40     if(x>9) Write(x/10);
41     putchar(x%10^48);
42 }
43 int Round(double x)
44 { 
45     if(fabs(x)<0.4) return 0;
46     return x>0?(int)(x+0.5):(int)(x-0.5);
47 }
48 void Prework()
49 {
50     register int i;
51     scanf("%d%d",&n,&m);
52     for(i=0;i<=n;i++) Read(a[i].x);
53     for(i=0;i<=m;i++) Read(b[i].x);
54     m+=n,n=1; while(n<=m) n<<=1;
55     for(i=1;i<=24;i++)
56         Sin[i]=sin(2*pai/(1<<i)),Cos[i]=cos(2*pai/(1<<i));
57 }
58 void Trans(cpx *cop,int len,int typ)
59 {
60     if(len==1) return;
61     register int i;
62     cpx t1[len>>1],t2[len>>1];
63     for(i=0;i<len;i++)
64         (i&1)?t2[i>>1]=cop[i]:t1[i>>1]=cop[i];
65     Trans(t1,len>>1,typ),Trans(t2,len>>1,typ);
66     int newl=len>>1,lgg=log2(len);
67     cpx omg={Cos[lgg],typ*Sin[lgg]},ori={1,0},tmp;
68     for(i=0;i<newl;i++,ori=ori*omg)
69         tmp=ori*t2[i],cop[i]=t1[i]+tmp,cop[i+newl]=t1[i]-tmp;
70 }
71 int main()
72 {
73     register int i;
74     Prework();
75     Trans(a,n,1),Trans(b,n,1);
76     for(i=0;i<n;i++) a[i]=a[i]*b[i];
77     Trans(a,n,-1);
78     for(i=0;i<=m;i++) Write(Round(a[i].x/n)),putchar(' ');
79     return 0;
80 }
View Code

非递归版(总共1500+ms，最慢的点不到600ms，没有加某“三步变两步”优化因为我不会)

 1 //Really Fast FFT
 2 #include<cmath>
 3 #include<cstdio>
 4 #include<cctype>
 5 #include<cstring>
 6 #include<algorithm>
 7 using namespace std;
 8 const int N=4e6+6,M=30;
 9 const double pai=acos(-1);
10 struct cpx
11 {
12     double x,y;
13 }a[N],b[N];
14 cpx operator + (cpx a,cpx b)
15 {
16     return (cpx){a.x+b.x,a.y+b.y};
17 } 
18 cpx operator - (cpx a,cpx b)
19 {
20     return (cpx){a.x-b.x,a.y-b.y};
21 }
22 cpx operator * (cpx a,cpx b)
23 {
24     double x1=a.x,x2=b.x,y1=a.y,y2=b.y;
25     return (cpx){x1*x2-y1*y2,x1*y2+x2*y1};
26 }
27 int n,m,rev[N];
28 double Sin[M],Cos[M];
29 void Read(double &x)
30 {
31     int ret=0; x=0;
32     char ch=getchar();
33     while(!isdigit(ch))
34         ch=getchar();
35     while(isdigit(ch))
36         ret=(ret<<1)+(ret<<3)+(ch^48),ch=getchar();
37     x=ret; return ;
38 }
39 void Write(int x)
40 {
41     if(x>9) Write(x/10);
42     putchar(x%10^48);
43 }
44 int Round(double x)
45 { 
46     if(fabs(x)<0.4) return 0;
47     return x>0?(int)(x+0.5):(int)(x-0.5);
48 }
49 void Prework()
50 {
51     register int i;
52     scanf("%d%d",&n,&m);
53     for(i=0;i<=n;i++) Read(a[i].x);
54     for(i=0;i<=m;i++) Read(b[i].x);
55     m+=n,n=1; while(n<=m) n<<=1;
56     for(i=0;i<n;i++)
57         rev[i]=(rev[i>>1]>>1)+(i&1)*(n>>1);
58     for(i=1;i<=24;i++)
59         Sin[i]=sin(2*pai/(1<<i)),Cos[i]=cos(2*pai/(1<<i));
60 }
61 void Trans(cpx *cop,int len,int typ)
62 {
63     register int i,j,k;
64     for(i=0;i<len;i++)
65         if(rev[i]>i) swap(cop[i],cop[rev[i]]);
66     for(i=2;i<=len;i<<=1)
67     {
68         int lth=i>>1,lgg=log2(i);
69         cpx omg={Cos[lgg],Sin[lgg]*typ};
70         for(j=0;j<len;j+=i)
71         {
72             cpx ori={1,0},tmp;
73             for(k=j;k<j+lth;k++,ori=ori*omg)
74                 tmp=ori*cop[k+lth],cop[k+lth]=cop[k]-tmp,cop[k]=cop[k]+tmp;
75         }
76     }
77     if(typ==-1)
78         for(int i=0;i<=len;i++) cop[i].x/=len;
79 }
80 int main()
81 {
82     register int i;
83     Prework();
84     Trans(a,n,1),Trans(b,n,1);
85     for(i=0;i<n;i++) a[i]=a[i]*b[i];
86     Trans(a,n,-1);
87     for(i=0;i<=m;i++) Write(Round(a[i].x)),putchar(' ');
88     return 0;
89 }
View Code

2.快速数论变换(NTT)

①NTT用来干啥

“FFT太慢啦” “FN^3DP，我的FFT比NTT还快”

“FFT丢精度” “我写拆系数FFT”

算了这个引入太沙雕了，直接讲吧

NTT能在模意义下快速完成多项式卷积

②NTT原理

类似FFT的，现在我们要找一个能在模意义下代替单位根的东西，它就是— —

原根(有关阶和原根的更多姿势请参考“知识汇总”)

然后别的跟FFT完全一致=。=

所以当模数是$a2^b+1$时不妨想想NTT

③NTT复杂度

时间复杂度$O(n\log n)$

④NTT具体实现

 1 //Simple NTT
 2 #include<cmath>
 3 #include<cstdio>
 4 #include<cctype>
 5 #include<cstring>
 6 #include<algorithm>
 7 using namespace std;
 8 const int N=4000006,mod=998244353;
 9 int n,m,G,Gi,Ni,rev[N],a[N],b[N],pw[30][2];
10 inline void Read(int &x)
11 {
12     x=0;
13     char ch=getchar();
14     while(!isdigit(ch))
15         ch=getchar();
16     while(isdigit(ch))
17         x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
18     return ;
19 }
20 void Write(int x)
21 {
22     if(x>9) Write(x/10);
23     putchar(x%10^48);
24 }
25 int Qpow(int x,int k)
26 {
27     if(k==1) return x;
28     int tmp=Qpow(x,k/2);
29     return k%2?1ll*tmp*tmp%mod*x%mod:1ll*tmp*tmp%mod;
30 }
31 void Prework()
32 {
33     register int i; 
34     Read(n),Read(m);
35     for(i=0;i<=n;i++) Read(a[i]);
36     for(i=0;i<=m;i++) Read(b[i]);
37     m+=n,n=1; while(n<=m) n<<=1;
38     for(i=1;i<n;i++)
39         rev[i]=(rev[i>>1]>>1)+(i&1)*(n>>1);
40     G=3,Gi=Qpow(G,mod-2),Ni=Qpow(n,mod-2);
41     for(int i=1;i<=24;i++)
42     {
43         pw[i][0]=Qpow(G,(mod-1)/(1<<i));
44         pw[i][1]=Qpow(Gi,(mod-1)/(1<<i));
45     }
46 }
47 void Trans(int *arr,int len,int typ)
48 {
49     register int i,j,k;
50     for(i=0;i<len;i++)
51         if(rev[i]>i) swap(arr[rev[i]],arr[i]);
52     for(i=2;i<=len;i<<=1)
53     {
54         int lth=i>>1,ort=pw[(int)log2(i)][typ==-1];
55         for(j=0;j<len;j+=i)
56         {
57             int ori=1,tmp;
58             for(k=j;k<j+lth;k++,ori=1ll*ori*ort%mod)
59             {
60                 tmp=1ll*ori*arr[k+lth]%mod;
61                 arr[k+lth]=(arr[k]-tmp+mod)%mod;
62                 arr[k]=(arr[k]+tmp)%mod;
63             }
64         }
65     }
66     if(typ==-1)
67         for(i=0;i<=len;i++)
68             arr[i]=1ll*arr[i]*Ni%mod;
69 }
70 int main()
71 {
72     register int i;
73     Prework();
74     Trans(a,n,1),Trans(b,n,1);
75     for(i=0;i<n;i++) a[i]=1ll*a[i]*b[i]%mod;
76     Trans(a,n,-1);
77     for(i=0;i<=m;i++) Write(a[i]),putchar(' ');
78     return 0;
79 }
View Code

3.快速沃尔什-哈德玛变换(FWT)

①FWT用来干啥

快速进行集合的卷积

②FWT原理

类似FFT的(怎么又是类似FFT=。=)，我们构造一种变换把原来的集合展开成可以$O(n)$(即$O(2^n)$)乘完的形式，然后分治快速搞出来乘完再还原回去

怎么构造出来的不提，直接上结论(这里称原集合为$A$，分治的前半部分为$A_1$，后半部分为$A_2$)

或：$FWT_{or}(A)->(FWT_{or}(A_1),FWT_{or}(A_1+A_2))$

异或：$FWT_{xor}(A)->(FWT_{xor}(A_1)+FWT_{xor}(A_2)),FWT_{xor}(A_1)-FWT_{xor}(A_2))$

与：$FWT_{and}(A)->(FWT_{and}(A_1+A_2),FWT_{and}(A_2))$

IFWT显然就是反着做

或：$IFWT_{or}(A)->(IFWT_{or}(A_1),IFWT_{or}(A_2)-IFWT_{or}(A_1))$

异或：$IFWT_{xor}(A)->(\frac{IFWT_{xor}(A_1)+IFWT_{xor}(A_2)}{2}),\frac{IFWT_{xor}(A_1)-IFWT_{xor}(A_2)}{2})$

与：$IFWT_{and}(A)->(IFWT_{and}(A_1)-IFWT_{and}(A_2),IFWT_{and}(A_2))$

③FWT复杂度

时间复杂度$O(n2^n)$

④FWT具体实现

 1 #include<cstdio>
 2 #include<cctype>
 3 #include<cstring>
 4 #include<algorithm>
 5 using namespace std;
 6 const int N=(1<<18)+10,mod=998244353;
 7 int n,tot,inv,m1[N],m2[N],a[N],b[N];
 8 inline void read(int &x)
 9 {
10     x=0;
11     char ch=getchar();
12     while(!isdigit(ch))
13         ch=getchar();
14     while(isdigit(ch))
15         x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
16 }
17 void write(int x)
18 {
19     if(x>9) write(x/10);
20     putchar(x%10^48);
21 }
22 int qpow(int x,int k)
23 {
24     if(k==1) return x;
25     int tmp=qpow(x,k/2);
26     return k%2?1ll*tmp*tmp%mod*x%mod:1ll*tmp*tmp%mod;
27 }
28 void transform(int *s,int x,int t)
29 {
30     register int i,j,k;
31     if(x==1)
32     {
33         for(i=2;i<=tot;i<<=1)
34         {
35             int len=i>>1;
36             for(j=0;j<tot;j+=i)
37                 for(k=j;k<j+len;k++)
38                     s[k+len]=(1ll*s[k+len]+1ll*s[k]*t+mod)%mod;
39         }
40     }
41     else if(x==2)
42     {
43         for(i=2;i<=tot;i<<=1)
44         {
45             int len=i>>1;
46             for(j=0;j<tot;j+=i)
47                 for(k=j;k<j+len;k++)
48                     s[k]=(1ll*s[k]+1ll*s[k+len]*t+mod)%mod;
49         }
50     }
51     else
52     {
53         for(i=2;i<=tot;i<<=1)
54         {
55             int len=i>>1,tmp;
56             for(j=0;j<tot;j+=i) 
57                 for(k=j;k<j+len;k++)
58                     tmp=s[k+len],s[k+len]=(s[k]-tmp+mod)%mod,s[k]=(s[k]+tmp)%mod;
59             if(t<0) for(j=0;j<tot;j++) s[j]=1ll*s[j]*inv%mod;
60         }
61     }
62 }
63 void Solve(int typ)
64 {
65     register int i;
66     for(i=0;i<tot;i++)
67         a[i]=m1[i],b[i]=m2[i];
68     transform(a,typ,1),transform(b,typ,1);
69     for(i=0;i<tot;i++) a[i]=1ll*a[i]*b[i]%mod;
70     transform(a,typ,-1);
71     for(i=0;i<tot;i++) write(a[i]),putchar(' '); 
72     puts("");
73 }
74 int main()
75 {
76     register int i;
77     read(n),tot=1<<n,inv=qpow(2,mod-2);
78     for(i=0;i<tot;i++) read(m1[i]);
79     for(i=0;i<tot;i++) read(m2[i]);
80     for(i=1;i<=3;i++) Solve(i);
81     return 0;
82 }
View Code

4.拉格朗日插值

①拉格朗日插值用来干啥

$O(n^2)$插值

②拉格朗日插值的原理

依次对每个点值构造，使得对于一个点来说这一项只有它有值，其他的点都是零

就是这样的一个东西：

$f(x)=\sum\limits_{i=0}^ny_i\prod_{i!=j}\frac{x-x_j}{x_i-x_j}$

③拉格朗日插值复杂度

上面说过了，时间复杂度$O(n^2)$

④拉格朗日插值具体实现

直接带要求的值进去，把值插出来↓

 1 #include<cstdio>
 2 #include<cstring>
 3 #include<algorithm>
 4 using namespace std;
 5 const int N=2005,mod=998244353;
 6 int n,k,x[N],y[N];
 7 void exGCD(int a,int b,int &x,int &y)
 8 {
 9     if(!b) x=1,y=0;
10     else exGCD(b,a%b,y,x),y-=a/b*x;
11 }
12 int Inv(int x)
13 {
14     int xx,yy;
15     exGCD(x,mod,xx,yy);
16     return (xx%mod+mod)%mod;
17 }
18 int Lagrange(int p,int q)
19 {
20     int ret=0;
21     for(int i=1;i<=p;i++)
22     {
23         int l1=y[i],l2=1;
24         for(int j=1;j<=p;j++)
25             if(i!=j)
26             {
27                 l1=1ll*l1*(q-x[j])%mod;
28                 l2=1ll*l2*(x[i]-x[j]+mod)%mod;
29             }
30         ret=(ret+1ll*l1*Inv(l2)%mod)%mod;
31     }
32     return (ret%mod+mod)%mod;
33 }
34 int main()
35 {
36     scanf("%d%d",&n,&k);
37     for(int i=1;i<=n;i++)
38         scanf("%d%d",&x[i],&y[i]);
39     printf("%d",Lagrange(n,k));
40     return 0;
41 }
View Code

把多项式插出来↓

 1 #include<cstdio>
 2 #include<cstring>
 3 #include<algorithm>
 4 #define O 1ll
 5 using namespace std;
 6 const int N=2005,mod=998244353;
 7 int n,k,x[N],y[N],num[N],tmp[N],res[N],inv[N];
 8 void Add(int &x,int y)
 9 {
10     x+=y;
11     if(x>=mod) x-=mod;
12 }
13 void exGCD(int a,int b,int &x,int &y)
14 {
15     if(!b) x=1,y=0;
16     else exGCD(b,a%b,y,x),y-=a/b*x;
17 }
18 int Inv(int x)
19 {
20     int xx,yy;
21     exGCD(x,mod,xx,yy);
22     Add(xx,mod); return xx;
23 }
24 void Lagrange()
25 {
26     for(int i=1;i<=n;i++)
27     {
28         int den=1,lst=0;
29         for(int j=1;j<=n;j++)
30             if(i!=j) den=O*den*(x[i]-x[j]+mod)%mod;
31         den=O*y[i]*Inv(den)%mod;
32         for(int j=0;j<n;j++)
33         {
34             tmp[j]=O*(num[j]-lst+mod)*inv[i]%mod;
35             Add(res[j],O*den*tmp[j]%mod),lst=tmp[j];
36         }
37     }
38 }
39 void Pre()
40 {
41     num[0]=1;
42     for(int i=1;i<=n;swap(num,tmp),i++)
43     {
44         tmp[0]=0; inv[i]=Inv(mod-x[i]);
45         for(int j=1;j<=i;j++) tmp[j]=num[j-1];
46         for(int j=0;j<=i;j++) Add(tmp[j],O*num[j]*(mod-x[i])%mod);
47     }
48 }
49 int Calc(int x)
50 {
51     int ret=0,var=1;
52     for(int i=0;i<n;var=O*var*x%mod,i++)
53         Add(ret,O*var*res[i]%mod);
54     return ret;
55 }
56 int main()
57 {
58     scanf("%d%d",&n,&k);
59     for(int i=1;i<=n;i++)
60         scanf("%d%d",&x[i],&y[i]);
61     Pre(),Lagrange(),printf("%d",Calc(k));
62     return 0;
63 }
View Code

5.特征多项式

①用来干啥

优化常系数齐次线性递推，我们如果直接推是$O(nk)$的，矩乘是$O(k^3\log n)$的，那如果$k$达到$1e3$的级别怎么做呢？

于是有了这个东西— —等等，什么叫常系数齐次线性递推？

常系数— —就是字面意思(转移时系数是常数)，齐次— —这里你可以理解为系数下标+转移下标恒定，线性递推— —emmm......

也就是说是长这样的递推式：$f_n=\sum\limits_{i=1}^k a_i f_{n-i}$

②原理

这里需要Cayley-Hamilton定理：一个矩阵$A$的特征多项式$P(x)$满足$P(A)=0$

博主太菜，没有证明

假设我们现在已经求出来了一个常系数齐次线性递推的特征多项式，下一步怎么做？

发现因为$P(A)=0$，所以我们可以减去若干倍的$P(A)$答案不变，这就相当于在模$P(x)$意义下做，用$O(n^2)$的暴力取模+快速幂即可$O(k^2\log n)$求解

说了这些之后，常系数齐次线性递推的特征多项式怎么求？

不用求，它们都长这样：$P(x)=x^k-\sum\limits_{i=1}^k a_ix^{k-i}$

③复杂度

至少能做到$O(k^2\log n)$

如果你真的非常**的话，也可以用多项式卷积优化多项式取模做到$O(k\log k\log n)$

④具体实现

例题：NOI2017泳池

(但是这题本身的DP也不好推=。=）

 1 #include<cstdio>
 2 #include<cstring>
 3 #include<algorithm>
 4 using namespace std;
 5 const int N=1005,mod=998244353;
 6 int m,k,x,y,p,pp,dp[N],pw[N],mo[N],mu[N],cal[2*N],ret[N],f[N][N];
 7 void Add(int &x,int y)
 8 {
 9     x+=y;
10     if(x>=mod) x-=mod;
11 }
12 int Qpow(int x,int k)
13 {
14     if(k==1) return x;
15     int tmp=Qpow(x,k/2);
16     return k%2?1ll*tmp*tmp%mod*x%mod:1ll*tmp*tmp%mod;
17 }
18 void Polymul(int *a,int *b,int *m,int *r,int len)
19 {
20     for(int i=0;i<=len;i++)
21         for(int j=0;j<=len;j++)
22             Add(cal[i+j],1ll*a[i]*b[j]%mod);
23     for(int i=2*len;i>=len;i--)
24         if(cal[i]) 
25             for(int j=0;j<=len;j++)
26                 Add(cal[i+j-len],mod-1ll*cal[i]*m[j]%mod);
27     for(int i=0;i<=len;i++) r[i]=cal[i],cal[i]=0;
28 }
29 int Solve(int n)
30 {
31     for(int i=1;i<=n+2;i++) f[0][i]=1;
32     for(int i=1;i<=n;i++)
33     {
34         for(int j=2;i*(j-1)<=n;j++)
35         {
36             int tmp=0;
37             for(int h=1;h<=i;h++)
38                 Add(tmp,1ll*f[h-1][j+1]*f[i-h][j]%mod);
39             f[i][j]=1ll*tmp*pw[j-1]%mod*pp%mod;
40         }
41         f[i][n/i+2]=0;
42         for(int j=n/i+1;j>=2;j--)
43             Add(f[i][j],f[i][j+1]);
44     }
45 //    for(int i=0;i<=n;i++,puts(""))
46 //        for(int j=0;j<=n;j++)
47 //            printf("%d ",f[i][j]);
48     memset(dp,0,sizeof dp),dp[0]=1;
49     for(int i=1;i<=n;i++)
50     {
51         for(int j=1;j<=i;j++)
52             Add(dp[i],1ll*dp[i-j]*f[j-1][2]%mod*pp%mod);
53         Add(dp[i],f[i][2]);
54     }
55     for(int i=0;i<=n;i++) 
56         mo[i]=1ll*(mod-f[n-i][2])*pp%mod; mo[++n]=1;
57 //    for(int i=0;i<=n;i++) printf("%d ",mo[i]);puts("");
58     memset(mu,0,sizeof mu),memset(ret,0,sizeof ret);
59     mu[1]=1,ret[0]=1; int ept=m,ans=0; 
60     while(ept)
61     {
62         if(ept&1) Polymul(ret,mu,mo,ret,n);
63         Polymul(mu,mu,mo,mu,n),ept>>=1;
64     }//for(int i=0;i<=n;i++) printf("%d ",ret[i]);puts("");
65     for(int i=0;i<n;i++)
66         Add(ans,1ll*dp[i]*ret[i]%mod);
67     return ans;
68 }
69 int main()
70 {
71     scanf("%d%d%d%d",&m,&k,&x,&y);
72     p=1ll*x*Qpow(y,mod-2)%mod;
73     pp=(1-p+mod)%mod,pw[0]=1;
74     for(int i=1;i<=k;i++) 
75         pw[i]=1ll*pw[i-1]*p%mod;
76     printf("%d",(Solve(k)-Solve(k-1)+mod)%mod);
77     return 0;
78 }
View Code

 

 

posted @ 2018-12-11 10:33  Speranza_Leaf  阅读(202)  评论(0)    收藏  举报