[海军国际项目办公室]生之花
生之花
S
Y
D
e
v
i
l
\rm\color{black}{S}\color{red}{YDevil}
SYDevil出的恶心题，题面也太臭了。
题目描述
这里略去一张奈亚子。
题解
不得不说，手*魔鬼写的题面真的太恶臭了。
首先，对于
f
(
A
,
n
)
f(A,n)
f(A,n)的转移，我们很容易想到背包上。
我们定义
d
p
i
,
j
dp_{i,j}
dpi,j前
i
i
i个中选了
j
j
j个加入我们的函数转移的所有方案的函数值之和，显然，转移中还会用到选择的总方案，所以我们还需记录
g
i
,
j
g_{i,j}
gi,j表示选择函数转移的方案数。
容易得到转移方程式，
d
p
i
,
j
=
a
d
p
i
−
1
,
j
−
1
+
d
p
i
−
1
,
j
+
b
g
i
−
1
,
j
dp_{i,j}=adp_{i-1,j-1}+dp_{i-1,j}+bg_{i-1,j}
dpi,j=adpi−1,j−1+dpi−1,j+bgi−1,j
g
i
,
j
=
g
i
−
1
,
j
−
1
+
g
i
−
1
,
j
g_{i,j}=g_{i-1,j-1}+g_{i-1,j}
gi,j=gi−1,j−1+gi−1,j
这样的话单次询问是
O
(
n
2
)
O\left(n^2\right)
O(n2)，明显是可以优化的。
我们可以将我们的带
x
x
x的部分与不带
x
x
x的部分分别计算，分别记作
f
a
i
,
j
fa_{i,j}
fai,j与
f
b
i
,
j
fb_{i,j}
fbi,j。
f
a
fa
fa中只记录
x
x
x的系数，
f
b
fb
fb记录不带
x
x
x的部分的大小。
显然，最后的答案等于
x
f
a
+
f
b
xfa+fb
xfa+fb，转移也比较好想，
f
a
i
,
j
=
a
f
a
i
−
1
,
j
−
1
+
f
a
i
−
1
,
j
,
f
b
i
,
j
=
a
f
b
i
−
1
,
j
−
1
+
f
b
i
−
1
,
j
+
b
g
i
−
1
,
j
−
1
fa_{i,j}=afa_{i-1,j-1}+fa_{i-1,j},fb_{i,j}=afb_{i-1,j-1}+fb_{i-1,j}+bg_{i-1,j-1}
fai,j=afai−1,j−1+fai−1,j,fbi,j=afbi−1,j−1+fbi−1,j+bgi−1,j−1
最后可以用差分求出我们询问的区间。
显然，上面的转移过程可以通过矩阵进行优化。
我们可以将第一维去掉，记
f
i
f_{i}
fi表示
(
f
a
i
,
f
b
i
,
g
i
)
\left(fa_{i},fb_{i},g_{i}\right)
(fai,fbi,gi)。
当我们加入
(
a
j
,
b
j
)
(a_{j},b_{j})
(aj,bj)时有转移矩阵，
A
j
=
(
a
j
,
0
,
0
0
,
a
j
,
0
0
,
b
j
,
1
)
A_{j}=\left(\begin{array}{cc}a_{j},0,0\\0,a_{j},0\\0,b_{j},1\end{array}\right)
Aj=⎝⎛aj,0,00,aj,00,bj,1⎠⎞，表示我们我们增加一个时改变的方案。
显然，转移是
f
i
=
A
j
f
i
−
1
+
I
f
i
f_{i}=A_{j}f_{i-1}+If_{i}
fi=Ajfi−1+Ifi
由于我们必须维护选择的点的个数这一维，所以我们必须将
A
A
A与
I
I
I的转移分开，维护所有的
F
i
F_{i}
Fi，但我们有没有更简洁的维护方法呢？
其实有了上面的式子，是很容易向生成函数上靠的。
定义
F
=
∑
f
i
x
i
F=\sum f_{i}x^i
F=∑fixi，每次转移相当于给我们的
F
F
F乘上
(
I
+
A
j
x
)
(I+A_{j}x)
(I+Ajx)。
初始矩阵
f
0
=
(
1
,
0
,
1
)
f_0=\left(1,0,1\right)
f0=(1,0,1)，答案是
∑
i
=
l
r
(
f
0
∏
(
I
+
A
j
x
)
)
[
x
i
]
\sum_{i=l}^{r} (f_0\prod(I+A_{jx}))[x^i]
∑i=lr(f0∏(I+Ajx))[xi]
显然，我们目标的多项式
f
0
∏
(
I
+
A
j
x
)
f_0\prod(I+A_{jx})
f0∏(I+Ajx)是可以算出来的，很容易想到分治的方法去求。
但我们的矩阵是不能直接卷积的，我的模数又是一个输入的数，相当于我们只能把矩阵中所有的数都 MTT 一次，再按矩阵乘法的形式合起来的话也可以做，但这样就常数太大了，据
R
a
i
n
y
b
u
n
n
y
\rm\color{black}{R}\red{ainybunny}
Rainybunny计算，每层乘法平均要进行
7.5
7.5
7.5次 FFT ，常数未必也太大了。
我们考虑还有没有其他方法可以求着东西。
学过小学乘法的人都知道，我们一般进行乘法是像竖式乘法一样做的，如果我们要进行
A
×
B
A\times B
A×B，我们就要将
A
A
A的每一位与
B
B
B的每一位相乘，再加起来，这样的话，如果
A
A
A与
B
B
B长度都是
n
n
n的话，我们整个乘法的时间复杂度是
O
(
n
2
)
O\left(n^2\right)
O(n2)的。
但其实我们可以对这个乘法进行优化。
我们可以
A
A
A
B
B
B都分成两段，
A
=
1
0
m
A
′
+
A
′
′
,
B
=
1
0
m
B
′
+
B
′
′
A=10^mA'+A'',B=10^mB'+B''
A=10mA′+A′′,B=10mB′+B′′，两段长度相等。
显然，根据乘法分配律，
A
B
=
1
0
2
m
A
′
B
′
+
1
0
m
(
A
′
B
′
′
+
A
′
′
B
′
)
+
A
′
′
B
′
′
AB=10^{2m}A'B'+10^m(A'B''+A''B')+A''B''
AB=102mA′B′+10m(A′B′′+A′′B′)+A′′B′′
如果我们要分别处理
A
′
B
′
,
A
′
B
′
′
,
A
′
′
B
′
,
A
′
′
B
′
′
A'B',A'B'',A''B',A''B''
A′B′,A′B′′,A′′B′,A′′B′′的话，显然并没有优化时间复杂度，但我们是否可以将
A
′
B
′
′
A'B''
A′B′′和
A
′
′
B
′
A''B'
A′′B′合并起来处理。
由于
(
A
′
+
A
′
′
)
(
B
′
+
B
′
′
)
=
A
′
B
′
+
A
′
′
B
′
′
+
(
A
′
B
′
′
+
A
′
′
B
′
)
(A'+A'')(B'+B'')=A'B'+A''B''+(A'B''+A''B')
(A′+A′′)(B′+B′′)=A′B′+A′′B′′+(A′B′′+A′′B′)，所以我们可以用
A
′
+
A
′
′
A'+A''
A′+A′′乘上
B
′
+
B
′
′
B'+B''
B′+B′′减去
A
′
B
′
A'B'
A′B′和
A
′
′
B
′
′
A''B''
A′′B′′，得到我们中间这一项的和。
这样我们就从原先的四次乘法变成了上次乘法加两次加法，如果只这样做一次的话最多只优化了常数，但我们完全可以再多做几次嘛，将我们的乘法不断递归下去，那么就有，
T
(
n
)
=
3
T
(
n
2
)
+
2
n
T(n)=3T(\frac{n}{2})+2n
T(n)=3T(2n)+2n，跟据主定理可以求出，
T
(
n
)
=
n
log
2
3
≈
n
n
T(n)=n^{\log_{2}3}\approx n\sqrt{n}
T(n)=nlog23≈nn。
这也就是
k
a
r
a
t
s
u
b
a
karatsuba
karatsuba乘法，一种常用于高精乘法的操作。
事实上该乘法还可以用于优化矩阵乘法，将矩阵乘法同样分成几部分去乘，可以将其从
O
(
n
3
)
O\left(n^3\right)
O(n3)优化到
O
(
n
log
2
7
)
O\left(n^{\log_{2}7}\right)
O(nlog27)。看起来没啥*用
这中做法显然可以迁移到我们这种矩阵的背包上面。
我们也可以像这样分治去进行矩阵乘法背包合并，算出
∏
(
I
+
A
x
)
\prod(I+Ax)
∏(I+Ax)，记作
G
G
G。
将
g
i
g_{i}
gi做个前缀和，就可以回答询问了。
事实上上面的做法在
n
n
n比较小的时候跑得比较慢，但
n
n
n较小的部分我们可以暴力合并，较大的部分我们再这样分治下去就可以了。
时间复杂度
O
(
n
log
2
3
log
n
)
O\left(n^{\log_{2}3}\log\,n\right)
O(nlog23logn)。
虽然常数有点大，但还是吊打暴力。
源码
#include<bits/stdc++.h>
using namespace std;
#define MAXN 20005
#define lowbit(x) (x&-x)
#define pb push_back
#define mkpr make_pair
#define fir first
#define sec second
#define lson (rt<<1)
#define rson (rt<<1|1)
#define debug(x) cerr<<#x<<"="<<x<<'\n'
typedef long long LL;
typedef unsigned long long uLL;
const LL INF=0x3f3f3f3f3f3f3f3f;
const int mo=998244353;
const int inv2=499122177;
const int jzm=2333;
const int n1=50;
const int zero=10000;
const int orG=3,invG=332748118;
const double Pi=acos(-1.0);
const double eps=1e-5;
typedef pair<int,int> pii;
template<typename _T>
_T Fabs(_T x){return x<0?-x:x;}
template<typename _T>
void read(_T &x){
_T f=1;x=0;char s=getchar();
while(s>'9'||s<'0'){if(s=='-')f=-1;s=getchar();}
while('0'<=s&&s<='9'){x=(x<<3)+(x<<1)+(s^48);s=getchar();}
x*=f;
}
template<typename _T>
void print(_T x){if(x>9)print(x/10);putchar(x%10+'0');}
LL gcd(LL a,LL b){return !b?a:gcd(b,a%b);}
int add(int x,int y,int p){return x+y<p?x+y:x+y-p;}
void Add(int &x,int y,int p){x=add(x,y,p);}
int qkpow(int a,int s,int p){int t=1;while(s){if(s&1LL)t=1ll*a*t%p;a=1ll*a*a%p;s>>=1LL;}return t;}
int n,m,P,stak,a[MAXN],b[MAXN];
struct matrix{
int c[5][5];matrix(){for(int i=1;i<4;i++)for(int j=1;j<4;j++)c[i][j]=0;}
void clear(){for(int i=1;i<4;i++)for(int j=1;j<4;j++)c[i][j]=0;}
matrix operator + (const matrix &rhs){
matrix res;res.clear();
for(int i=1;i<4;i++)for(int j=1;j<4;j++)
res.c[i][j]=add(c[i][j],rhs.c[i][j],P);
return res;
}
matrix operator - (const matrix &rhs){
matrix res;res.clear();
for(int i=1;i<4;i++)for(int j=1;j<4;j++)
res.c[i][j]=add(c[i][j],P-rhs.c[i][j],P);
return res;
}
matrix operator * (const matrix &rhs){
matrix res;res.clear();
for(int i=1;i<4;i++)
for(int k=1;k<4;k++)if(c[i][k])
for(int j=1;j<4;j++)
Add(res.c[i][j],1ll*c[i][k]*rhs.c[k][j]%P,P);
return res;
}
void print(){
for(int i=1;i<4;i++,puts(""))
for(int j=1;j<4;j++)printf("%d ",c[i][j]);
}
}sum[MAXN],I,A;
vector<matrix>tr[MAXN<<2],sta[MAXN],tmp;
void work(int x,int y){
int sizx=sta[x].size(),sizy=sta[y].size();
if(sizx<8||sizy<8){
tmp.clear();tmp.resize(sizx+sizy-1);
for(int i=0;i<sizx;i++)
for(int j=0;j<sizy;j++)
tmp[i+j]=(tmp[i+j]+(sta[x][i]*sta[y][j]));
sta[x]=tmp;return ;
}
int len=min(sizx/2,sizy/2);
int li1=++stak;sta[li1].resize(max(len,sizx-len));
int li2=++stak;sta[li2].resize(len);
int li3=++stak;sta[li3].resize(sizx-len);
int ri1=++stak;sta[ri1].resize(max(len,sizy-len));
int ri2=++stak;sta[ri2].resize(len);
int ri3=++stak;sta[ri3].resize(sizy-len);
for(int i=0;i<(int)sta[li1].size();i++)
if(i<len&&i+len<sizx)
sta[li1][i]=sta[x][i]+sta[x][i+len],
sta[li2][i]=sta[x][i],sta[li3][i]=sta[x][i+len];
else sta[li1][i]=sta[x][i+len],sta[li3][i]=sta[x][i+len];
for(int i=0;i<(int)sta[ri1].size();i++)
if(i<len&&i+len<sizy)
sta[ri1][i]=sta[y][i]+sta[y][i+len],
sta[ri2][i]=sta[y][i],sta[ri3][i]=sta[y][i+len];
else sta[ri1][i]=sta[y][i+len],sta[ri3][i]=sta[y][i+len];
work(li3,ri3);sta[ri3].clear();stak--;
work(li2,ri2);sta[ri2].clear();stak--;
work(li1,ri1);sta[ri1].clear();stak--;
sta[x].clear();sta[x].resize(sizx+sizy-1);
for(int i=0;i<(int)sta[li2].size();i++)
sta[x][i]=sta[x][i]+sta[li2][i],
sta[li1][i]=sta[li1][i]-sta[li2][i];
for(int i=0;i<(int)sta[li3].size();i++)
sta[x][i+len+len]=sta[x][i+len+len]+sta[li3][i],
sta[li1][i]=sta[li1][i]-sta[li3][i];
for(int i=0;i<(int)sta[li1].size();i++)
sta[x][i+len]=sta[x][i+len]+sta[li1][i];
sta[li3].clear();stak--;
sta[li2].clear();stak--;
sta[li1].clear();stak--;
}
void sakura(int rt,int l,int r){
int mid=l+r>>1;tr[rt].resize(r-l+2);
if(l==r){
tr[rt][0]=I;tr[rt][1].c[3][3]=1;
tr[rt][1].c[1][1]=tr[rt][1].c[2][2]=a[l];
tr[rt][1].c[3][2]=b[l];return ;
}
sakura(lson,l,mid);sakura(rson,mid+1,r);
int lid=++stak;sta[lid]=tr[lson];
int rid=++stak;sta[rid]=tr[rson];
work(lid,rid);sta[rid].clear();stak--;
for(int i=0;i<=r-l+1;i++)tr[rt][i]=sta[lid][i];
sta[lid].clear();stak--;
}
signed main(){
freopen("cthulhu.in","r",stdin);
freopen("cthulhu.out","w",stdout);
read(n);read(m);read(P);
for(int i=1;i<=3;i++)I.c[i][i]=1;
for(int i=1;i<=n;i++)read(a[i]),read(b[i]),a[i]%=P,b[i]%=P;
sakura(1,1,n);
for(int i=0;i<=n;i++)sum[i]=tr[1][i];
for(int i=1;i<=n;i++)sum[i]=sum[i-1]+sum[i];
for(int i=1;i<=m;i++){
int x,l,r;read(x);read(l);read(r);
A.clear();A.c[1][1]=x;A.c[1][2]=0;A.c[1][3]=1;
A=A*(sum[r]-sum[l-1]);
print(add(A.c[1][1],A.c[1][2],P));putchar('\n');
}
return 0;
}
浙公网安备 33010602011771号