LOJ6183 看无可看

Problem

https://loj.ac/problem/6183

Solution

考虑这个递推式的特征方程$x^2=2x+3$,解出来$x1=3,x2=-1$

设递推式为$f_n=\alpha x_1^n+\beta x_2^n$
$$
\alpha +\beta=f_0\\
\alpha x_1+\beta x_2=f1\\
\alpha=\frac{f1-f0x_2}{x_1-x_2},\beta=\frac{f0x1-f1}{x1-x2}\\
\alpha=\frac{f1+f0}{4},\beta=\frac{3f0-f1}{4}\\
$$
最后要求的答案是
$$
\sum_{S’\in S,|S’|=k}\alpha x_1^{\sum_{x\in S’}x}+\beta x_2^{\sum_{x\in S’}x}\
=\alpha \sum_{S’\in S,|S’|=k}\prod_{x\in S’}3^x+\beta \sum_{x\in S’}\prod_{x\in S’}(-1)^x
$$

然后发现要算$\sum_{S’\in S,|S’|=k}\prod_{x\in S’}3^x$这么个东西。

考虑组合意义,设$f[i][j]$表示前$i$个数,选了$j$个的答案。

有$f[i][j]=f[i-1][j]+f[i-1][j-1]*(3^{a_i})$,有点像一个背包,表示选/不选。

可以对f[i]构造生成函数$F_i$为$1+3x$,这个转移就成了一个卷积的形式,$F_i=F_{i-1}*F_{i-1}$

然后把所有东西卷一下就好了,注意要分治fft,不然是$n^2\log$的。

由于毒瘤模数就写了mtt,但是要跑3s。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
#include <bits/stdc++.h>
using namespace std;
#define IL inline
IL int read() { int x; int _w = scanf("%d", &x); return x; }
IL int read(int &x) { int _w=scanf("%d", &x); }
IL void write(int x) { printf("%d\n", x); }
#define rep(i, j, k) for (int i = j; i <= k; ++i)
#define repd(i, j, k) for (int i = j; i >= k; --i)
#define pb push_back
#define mp make_pair
#define fr first
#define se second
#define ls (p<<1)
#define rs (p<<1|1)
#define ef else if
#define re real
#define im imag
#define fre(x) freopen(x".in","r",stdin),freopen(x".out", "w", stdout)
#define popcnt(x) __builtin_popcount(x)
#define ls (p<<1)
#define rs (p<<1|1)
#define cpx my_complex

typedef long long ll;
const int mod=99991;
const long double pi=acos(-1);
const int phip=99990;
const int inv4=24998;
IL int Mod(ll x){x=(x%mod+mod)%mod; return x;}
const int maxn=2e5+10;
int n, k;
ll a[maxn], F0, F1;
int fpw(int x,ll y,int r=1){y%=phip;for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)r=1ll*r*x%mod;assert(r>=0);return r;}

int fun(ll x){
if(!(x&1))return ((inv4*(F1+F0)%mod*fpw(3,x)%mod)+(inv4*(((3*F0-F1)%mod+mod)%mod)%mod))%mod;
else return Mod((inv4*(F1+F0)%mod*fpw(3,x)%mod)-(inv4*(((3*F0-F1)%mod+mod)%mod)%mod));
}

class my_complex {
public:
long double real,imag;
my_complex(long double a=0, long double b=0){real=a,imag=b;}
my_complex operator * (my_complex &a){return my_complex(real*a.real-imag*a.imag,real*a.imag+imag*a.real);}
my_complex operator - (my_complex &a){return my_complex(real-a.real,imag-a.imag);}
my_complex operator + (my_complex &a){return my_complex(real+a.real,imag+a.imag);}
};
int rev[maxn];
vector<cpx>f1,f2,g1,g2,h1,h2,h3,h4;
vector<int>tr[maxn<<2], gr[maxn<<2];
IL void dft(vector<cpx>&x, int o) {
int sz=x.size();
rep(i,0,sz-1)if(i<rev[i])swap(x[i],x[rev[i]]);
for(int i=1;i<sz;i<<=1){
cpx e(cos(pi/i),o*sin(pi/i));
for(int j=0;j<sz;j+=(i<<1)){
cpx xx(1,0);
for(int k=j;k<j+i;++k,xx=xx*e){
cpx a1=x[k];
cpx a2=xx*x[k+i];
x[k]=a1+a2;x[k+i]=a1-a2;
}
}
}
}
IL void conv(vector<int>&x, vector<int>&y, vector<int>&z){
int n=x.size(), m = y.size(), s = n + m - 1;
while(s!=(s&-s))s+=s&-s;
int b=1; while(1<<b<s)b++;
rep(i,0,s-1)rev[i]=rev[i>>1]>>1|((i&1)<<(b-1));
f1.clear(), f2.clear(), g2.clear(), g1.clear();

rep(i,0,n-1)f1.pb(cpx(x[i]&32767,0));rep(i,n,s-1)f1.pb(cpx(0,0));
rep(i,0,n-1)f2.pb(cpx(x[i]>>15,0));rep(i,n,s-1)f2.pb(cpx(0,0));
rep(i,0,m-1)g1.pb(cpx(y[i]&32767,0));rep(i,m,s-1)g1.pb(cpx(0,0));
rep(i,0,m-1)g2.pb(cpx(y[i]>>15,0));rep(i,m,s-1)g2.pb(cpx(0,0));

dft(f1,1);dft(g1,1);dft(f2,1);dft(g2,1);

h1.clear(); h1.resize(s); h2.clear(); h2.resize(s); h3.clear(); h3.resize(s); h4.clear(); h4.resize(s);

rep(i,0,s-1)h1[i]=f1[i]*g1[i];rep(i,0,s-1)h2[i]=f1[i]*g2[i];
rep(i,0,s-1)h3[i]=f2[i]*g1[i];rep(i,0,s-1)h4[i]=f2[i]*g2[i];

dft(h1,-1);dft(h2,-1);dft(h3,-1);dft(h4,-1);
z.clear(); z.resize(s);
rep(i,0,s-1) {
ll a1=ll(h1[i].re/s+0.5)%mod;ll a2=ll(h2[i].re/s+0.5)%mod;
ll a3=ll(h3[i].re/s+0.5)%mod;ll a4=ll(h4[i].re/s+0.5)%mod;
z[i]=(a1+((a2+a3)<<15)+(a4<<30))%mod;
}
}
IL void bd(int p,int x,int y,int base,int base2){
if(x==y){tr[p].pb(1);tr[p].pb(fpw(base,a[x]));
gr[p].pb(1);gr[p].pb(fpw(base2,a[x]));return;}
int mid=(x+y)>>1;
bd(ls,x,mid,base,base2);bd(rs,mid+1,y,base,base2);
conv(tr[ls],tr[rs],tr[p]);
conv(gr[ls],gr[rs],gr[p]);
}
int main(){
//fre("see");
//freopen("0.in","r",stdin);
read(n),read(k);
rep(i,1,n)a[i]=read();
F0=read(),F1=read();
int ans=0;
if(n<=0){
rep(i,0,(1<<n)-1){
int all=popcnt(i);
if(all!=k)continue;
ll tot=0;
rep(j,1,n){
if(i&(1<<j-1)){tot+=a[j];}
else continue;
}
ans=ans+fun(tot)%mod;
ans%=mod;
}
printf("%d",ans);
} else {
bd(1,1,n,3,mod-1);
int a1=tr[1][k];int a2=gr[1][k];
cerr<<a1<<' '<<a2<<endl;
int ans=((inv4*(F1+F0)%mod*a1%mod)+((inv4*(((3*F0-F1)%mod+mod)%mod)%mod)*a2%mod))%mod;
printf("%d",ans);
}
return 0;
}

可以对f[i]构造生成函数$F_i$为$1+3x$,这个转移就成了一个卷积的形式,$F_i=F_{i-1}*F_{i-1}$

然后把所有东西卷一下就好了,注意要分治fft,不然是$n^2\log$的。

由于毒瘤模数就写了mtt,但是要跑3s。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
#include <bits/stdc++.h>
using namespace std;
#define IL inline
IL int read() { int x; int _w = scanf("%d", &x); return x; }
IL int read(int &x) { int _w=scanf("%d", &x); }
IL void write(int x) { printf("%d\n", x); }
#define rep(i, j, k) for (int i = j; i <= k; ++i)
#define repd(i, j, k) for (int i = j; i >= k; --i)
#define pb push_back
#define mp make_pair
#define fr first
#define se second
#define ls (p<<1)
#define rs (p<<1|1)
#define ef else if
#define re real
#define im imag
#define fre(x) freopen(x".in","r",stdin),freopen(x".out", "w", stdout)
#define popcnt(x) __builtin_popcount(x)
#define ls (p<<1)
#define rs (p<<1|1)
#define cpx my_complex

typedef long long ll;
const int mod=99991;
const long double pi=acos(-1);
const int phip=99990;
const int inv4=24998;
IL int Mod(ll x){x=(x%mod+mod)%mod; return x;}
const int maxn=2e5+10;
int n, k;
ll a[maxn], F0, F1;
int fpw(int x,ll y,int r=1){y%=phip;for(;y;y>>=1,x=1ll*x*x%mod)if(y&1)r=1ll*r*x%mod;assert(r>=0);return r;}

int fun(ll x){
if(!(x&1))return ((inv4*(F1+F0)%mod*fpw(3,x)%mod)+(inv4*(((3*F0-F1)%mod+mod)%mod)%mod))%mod;
else return Mod((inv4*(F1+F0)%mod*fpw(3,x)%mod)-(inv4*(((3*F0-F1)%mod+mod)%mod)%mod));
}

class my_complex {
public:
long double real,imag;
my_complex(long double a=0, long double b=0){real=a,imag=b;}
my_complex operator * (my_complex &a){return my_complex(real*a.real-imag*a.imag,real*a.imag+imag*a.real);}
my_complex operator - (my_complex &a){return my_complex(real-a.real,imag-a.imag);}
my_complex operator + (my_complex &a){return my_complex(real+a.real,imag+a.imag);}
};
int rev[maxn];
vector<cpx>f1,f2,g1,g2,h1,h2,h3,h4;
vector<int>tr[maxn<<2], gr[maxn<<2];
IL void dft(vector<cpx>&x, int o) {
int sz=x.size();
rep(i,0,sz-1)if(i<rev[i])swap(x[i],x[rev[i]]);
for(int i=1;i<sz;i<<=1){
cpx e(cos(pi/i),o*sin(pi/i));
for(int j=0;j<sz;j+=(i<<1)){
cpx xx(1,0);
for(int k=j;k<j+i;++k,xx=xx*e){
cpx a1=x[k];
cpx a2=xx*x[k+i];
x[k]=a1+a2;x[k+i]=a1-a2;
}
}
}
}
IL void conv(vector<int>&x, vector<int>&y, vector<int>&z){
int n=x.size(), m = y.size(), s = n + m - 1;
while(s!=(s&-s))s+=s&-s;
int b=1; while(1<<b<s)b++;
rep(i,0,s-1)rev[i]=rev[i>>1]>>1|((i&1)<<(b-1));
f1.clear(), f2.clear(), g2.clear(), g1.clear();

rep(i,0,n-1)f1.pb(cpx(x[i]&32767,0));rep(i,n,s-1)f1.pb(cpx(0,0));
rep(i,0,n-1)f2.pb(cpx(x[i]>>15,0));rep(i,n,s-1)f2.pb(cpx(0,0));
rep(i,0,m-1)g1.pb(cpx(y[i]&32767,0));rep(i,m,s-1)g1.pb(cpx(0,0));
rep(i,0,m-1)g2.pb(cpx(y[i]>>15,0));rep(i,m,s-1)g2.pb(cpx(0,0));

dft(f1,1);dft(g1,1);dft(f2,1);dft(g2,1);

h1.clear(); h1.resize(s); h2.clear(); h2.resize(s); h3.clear(); h3.resize(s); h4.clear(); h4.resize(s);

rep(i,0,s-1)h1[i]=f1[i]*g1[i];rep(i,0,s-1)h2[i]=f1[i]*g2[i];
rep(i,0,s-1)h3[i]=f2[i]*g1[i];rep(i,0,s-1)h4[i]=f2[i]*g2[i];

dft(h1,-1);dft(h2,-1);dft(h3,-1);dft(h4,-1);
z.clear(); z.resize(s);
rep(i,0,s-1) {
ll a1=ll(h1[i].re/s+0.5)%mod;ll a2=ll(h2[i].re/s+0.5)%mod;
ll a3=ll(h3[i].re/s+0.5)%mod;ll a4=ll(h4[i].re/s+0.5)%mod;
z[i]=(a1+((a2+a3)<<15)+(a4<<30))%mod;
}
}
IL void bd(int p,int x,int y,int base,int base2){
if(x==y){tr[p].pb(1);tr[p].pb(fpw(base,a[x]));
gr[p].pb(1);gr[p].pb(fpw(base2,a[x]));return;}
int mid=(x+y)>>1;
bd(ls,x,mid,base,base2);bd(rs,mid+1,y,base,base2);
conv(tr[ls],tr[rs],tr[p]);
conv(gr[ls],gr[rs],gr[p]);
}
int main(){
//fre("see");
//freopen("0.in","r",stdin);
read(n),read(k);
rep(i,1,n)a[i]=read();
F0=read(),F1=read();
int ans=0;
if(n<=0){
rep(i,0,(1<<n)-1){
int all=popcnt(i);
if(all!=k)continue;
ll tot=0;
rep(j,1,n){
if(i&(1<<j-1)){tot+=a[j];}
else continue;
}
ans=ans+fun(tot)%mod;
ans%=mod;
}
printf("%d",ans);
} else {
bd(1,1,n,3,mod-1);
int a1=tr[1][k];int a2=gr[1][k];
cerr<<a1<<' '<<a2<<endl;
int ans=((inv4*(F1+F0)%mod*a1%mod)+((inv4*(((3*F0-F1)%mod+mod)%mod)%mod)*a2%mod))%mod;
printf("%d",ans);
}
return 0;
}