BZOJ4591: [Shoi2015]超能粒子炮·改

BZOJ4591: [Shoi2015]超能粒子炮·改

Description

曾经发明了脑洞治疗仪&超能粒子炮的发明家SHTSC又公开了他的新发明:超能粒子炮·改--一种可以发射威力更加强大的粒子流的神秘装置。
超能粒子炮·改相比超能粒子炮,在威力上有了本质的提升。它有三个参数n,k。它会向编号为0到k的位置发射威力为C(n,k) mod 2333的粒子流。
现在SHTSC给出了他的超能粒子炮·改的参数,让你求其发射的粒子流的威力之和模2333。

Input

第一行一个整数t。表示数据组数。
之后t行,每行二个整数n,k。含义如题面描述。
k<=n<=10^18,t<=10^5

Output

t行每行一个整数,表示其粒子流的威力之和模2333的值。

Sample Input

1
5 5

Sample Output

32

题解Here!

组合数还带质数取模,当然用$Lucas$!
前置技能——$Lucas$定理:$$C_{n}^{m}\%p=C_{n/p}^{m/p} imes C_{n\%p}^{m\%p}\%p$$
题目要求:$$sum_{i=0}^kC_{n}^i\%p$$
设$F(n,k)=sum_{i=0}^kC_{n}^i$。
则:$$F(n,k)=sum_{i=0}^kC_{n}^i=sum_{i=0}^kC_{n/p}^{i/p} imes C_{n\%p}^{i\%p}$$
这个东西好熟悉啊。。。
没错!这玩意包含一个式子:$n/p$!
这是啥?
知道莫比乌斯反演就会知道这玩意是数论分块的标志。
于是:
$$F(n,k)=C_{n/p}^0sum_{i=0}^{p-1}C_{n\%p}^i+C_{n/p}^1sum_{i=0}^{p-1}C_{n\%p}^i+cdots+C_{n/p}^{k/p}sum_{i=0}^{k\%p}C_{n\%p}^i$$
前面有$0$到$k/p-1$这些整块,先单独挑出来搞事。
我们发现它们有共同的系数$sum_{i=0}^{p-1}C_{n\%p}^i$。
于是我们可以将$sum_{i=0}^{p-1}C_{n\%p}^i$提出来:
$$sum_{i=0}^{p-1}C_{n\%p}^i imes (C_{n/p}^0+C_{n/p}^1+...C_{n/p}^{k/p-1})$$
再一波操作就成了:
$$F(n\%p,p-1) imes F(n/p,k/p-1)$$
再加上那个不完整的块:
$$sum_{i=0}^{k\%p}C_{n\%p}^i=F(n\%p,k\%p)$$
所以最终我们得到了一个式子:
$$F(n,k)=F(n\%p,p-1) imes F(n/p,k/p-1)+C_{n/p}^{k/p} imes F(n\%p,k\%p)$$
而$n\%p,k\%pleq 2333$。
所以$F(n\%p,p-1),F(n\%p,k\%p)$可以直接预处理出来。
至于那个$C_{n/p}^{k/p}$就丢给了$Lucas$。
时间复杂度$O(p^2+Tlog_{2333}^2n)$。
一开始把数组开小了还$RE imes3$。。。
附代码:
#include<iostream>
#include<algorithm>
#include<cstdio>
#define MAXN 2510
#define MOD 2333LL
using namespace std;
long long n,k;
long long C[MAXN][MAXN],F[MAXN][MAXN];
inline long long read(){
	long long date=0;char c=0;
	while(c<'0'||c>'9')c=getchar();
	while(c>='0'&&c<='9'){date=date*10+c-'0';c=getchar();}
	return date;
}
void make(){
	int m=MAXN-10;
	C[0][0]=1;
	for(int i=1;i<=m;i++){
		C[i][0]=C[i][i]=1;
		for(int j=1;j<i;j++)C[i][j]=(C[i-1][j-1]+C[i-1][j])%MOD;
	}
	for(int i=0;i<=m;i++){
		F[i][0]=1;
		for(int j=1;j<=m;j++)F[i][j]=(F[i][j-1]+C[i][j])%MOD;
	}
}
long long Lucas(long long n,long long m){
	if(n<m)return 0;
	if(m==0||m==n)return 1;
	if(m==1||m==n-1)return n;
	return C[n%MOD][m%MOD]*Lucas(n/MOD,m/MOD)%MOD;
}
long long solve(long long n,long long k){
	if(k<0)return 0;
	if(n==0||k==0)return 1;
	if(n<MOD&&k<MOD)return F[n][k];
	return (solve(n/MOD,k/MOD-1)*F[n%MOD][MOD-1]%MOD+Lucas(n/MOD,k/MOD)*F[n%MOD][k%MOD]%MOD)%MOD;
}
int main(){
	int t=read();
	make();
	while(t--){
		n=read();k=read();
		printf("%lld
",solve(n,k));
	}
	return 0;
}
原文地址:https://www.cnblogs.com/Yangrui-Blog/p/10545849.html