jzoj6005. 【PKUWC2019模拟2019.1.17】数学（生成函数&分治FFT）

题目描述

Description

Input

Output

Sample Input
4 5
Sample Output
2
样例解释：

Data Constraint

20%

随便搞个dp
设s[i][j]表示当前已经乘了i个数，末尾的数为j
显然可以写一个O(n³)的dp，然后用前缀和优化到O(n²)

45%

考虑求答案的生成函数

生成函数是™什么东西？
就是一个奇♂妙的函数，这个函数的x^p完全没有卵用
有用的只是每一项的系数
一般用于与排列组合有关的问题，x的指数根据不同的问题有不同的含义，在本题中就是k（选择的数的个数），那么对应的系数就是f(n,k)
则
$G(x)={\prod }_{i=0}^n(x+i)$
因为本题要求的是不被p的整除f(n,k)，当i=0时相当于整体位移，所以可以忽略
$G(x)={\prod }_{i=1}^n(x+i)$
那么45%就直接分治FFT/多模数NTT/MTT
mdzz我只会FFT

---

然而这个分治FFT并不是网上f(i)=∑f(j)g(i-j)的高级做法，实际就™是个分治
考虑把${\prod }_{i=1}^n(x+i)$从中间断开，之后再两两合并，具体实现可以用栈来模拟
（要预处理$\omega$不然会被卡爆）
这样做的深度log n，每层的总长度为n，再加上FFT的log n时间，总复杂度为O(n logn²)
妙啊

100%

因为n非常大（有10¹⁰⁰⁰），但是p很小，考虑转化成45%的问题
要求
$G(x)={\prod }_{i=1}^n(x+i)$
则
$G(x)={\prod }_{i=1}^{ap+b}(x+i)$（a=n/p，b=a%p，因为i是在模p意义下，所以直接判断i是否为0）
$G(x)=({\prod }_{i=1}^p(x+i){)}^a{\prod }_{i=1}^b(x+i)$（a=n/p，b=a%p）
因为${\prod }_{i=1}^p(x+i)$在1~p-1的取值中模p意义下与$x(x^{p-1}-1)$相等，所以二者等价（p是质数）
（其实中间有推导过程但既然大佬也不会那我就以后再说吧）
然后替换一下
$G(x)=(x(x^{p-1}-1){)}^a{\prod }_{i=1}^b(x+i)$
$G(x)=x^a(x^{p-1}-1{)}^a{\prod }_{i=1}^b(x+i)$
因为$x^a$是整体位移，不会影响到系数，所以可以忽略
$G(x)=(x^{p-1}-1{)}^a{\prod }_{i=1}^b(x+i)$
二项式展开
$G(x)=({\sum }_{i=0}^a{x}^{(p-1)i}{C}_a^i(-1{)}^{a-i}){\prod }_{i=1}^b(x+i)$

---

先考虑前半部分
因为$x^{(p-1)i}$和$(-1{)}^{a-i}$都不会把系数变为0，所以就变成了判断$C_a^i$是否为0
根据Lucas定理，$C_a^i=C_{\lfloor a/p\rfloor }^{\lfloor i/p\rfloor }\ast C_{a\%p}^{i\%p}$，i在不停模p的每一位都要小于a对应的位，其实就相当于小于a在p进制下的对应位
而位与位之间又不会互相影响（因为i只要有一位大于a$C_a^i$就会为0），所以可以直接计算
$\prod (mod[i]+1)$，就是把a在p进制下的每一位可能的方案数相乘（要加上0）
显然可以直接搞，时间复杂度为O(n的位数*log n)

---

考虑后半部分
因为前半部分两两x的指数相差至少为(p-1)，再乘上${\prod }_{i=1}^b(x+i)$就相当于给每一个x都乘了一遍
于是可以根据b的取值来讨论
①b<p-1
那么显然乘上之后不会影响其他的x^(p-1)i，所以可以直接分治FFT求出后半部分，然后与前半部分的方案相乘
②b=p-1
因为x^(p-1)ix^p-1=x^(p-1)(i+1)，所以会影响到其它的x^(p-1)i，不能直接乘
由前面可得，${\prod }_{i=1}^p(x+i)=x(x^{p-1}-1)$，则
${\prod }_{i=1}^{p-1}(x+i)=(x^{p-1}-1)$（在模p意义下）
当b=p-1时，原式等于
$G(x)=x^a(x^{p-1}-1{)}^a{\prod }_{i=1}^{p-1}(x+i)$
$G(x)=x^a(x^{p-1}-1{)}^a(x^{p-1}-1)$
去掉整体位移后
$G(x)=(x^{p-1}-1{)}^{a+1}$
做法同上

关于分治FFT的高级用法和多模数NTT/MTT和真·生成函数就以后再填坑吧

code

#include <iostream>
#include <cstdlib>
#include <cstring>
#include <cstdio>
#include <cmath>
#define fo(a,b,c) for (a=b; a<=c; a++)
#define fd(a,b,c) for (a=b; a>=c; a--)
#define mod 998244353
using namespace std;

struct C{
	long double a,b;
	C (long double _a=0,long double _b=0) {a=_a,b=_b;}
};
C operator+(C x,C y){return C(x.a+y.a,x.b+y.b);}
C operator-(C x,C y){return C(x.a-y.a,x.b-y.b);}
C operator*(C x,C y){return C(x.a*y.a-x.b*y.b,x.a*y.b+x.b*y.a);}

int n[1002];
C A[131073];
C c[18][131073];
int d[18];
int D[18];
C w[18][131073];
int b[4001];
int N,S,len,p,i,j,k,l,Mod,Len;
char ch;
long long s,ans;

void dft(int N,int S,int a,int type)
{
	int i,j,k,l,s1=N/2,s2=1,s3=2;
	C W,v,u;
	
	fo(i,0,N-1)
	{
		j=i,k=0;
		
		fo(l,1,S)
		{
			k=k*2+(j&1);
			j>>=1;
		}
		A[k]=c[a][i];
	}
	fo(i,0,N-1) c[a][i]=A[i];
	
	fo(i,1,S)
	{
		fo(j,0,s1-1)
		{
			fo(k,0,s2-1)
			{
				int _s=j*s3;
				W=w[i][k];W.b*=type;
				
				u=c[a][_s+k];
				v=c[a][_s+k+s2]*W;
				
				c[a][_s+k]=u+v;
				c[a][_s+k+s2]=u-v;
			}
		}
		
		s1>>=1;
		s2<<=1,s3<<=1;
	}
}

void chen(int N,int S,int a,int b)
{
	int i;
	
	dft(N,S,a,1);
	dft(N,S,b,1);
	
	fo(i,0,N-1)
	c[a][i]=c[a][i]*c[b][i];
	
	dft(N,S,a,-1);
	
	fd(i,N-1,0)
	{
		long long _s;
		
		c[a][i].a/=N;
		_s=floor(c[a][i].a+0.1);
		c[a][i].a=_s%p;
		
		c[a][i].b=0;
		c[b][i].a=0;
		c[b][i].b=0;
	}
}

void swap(int &x,int &y)
{
	int z=x;
	x=y;
	y=z;
}

void work() //(x+1)(x+2)...(x+Mod)
{
	int i,j,k,l,len=1;
	
	d[1]=2;
	D[1]=1;
	c[1][0].a=1;
	c[1][1].a=1;
	
	fo(i,2,Mod)
	{
		while (len>1 && d[len]==d[len-1])
		{
			chen(d[len-1]*2,D[len-1]+1,len-1,len);
			--len;
			d[len]*=2;
			++D[len];
		}
		
		++len;
		d[len]=2;
		D[len]=1;
		c[len][0].a=i;
		c[len][1].a=1;
	}
	
	while (len>1)
	{
		chen(N,S,len-1,len);
		--len;
	}
}

void work2()
{
	int i;
	
	while (len)
	{
		n[0]=0;
		fd(i,len,1)
		{
			n[i-1]+=(n[i]%p)*10;
			n[i]/=p;
		}
		
		while (len && !n[len]) --len;
		b[++Len]=n[0]/10;
	}
	
	Mod=b[1];
}

int main()
{
	freopen("math.in","r",stdin);
	freopen("math.out","w",stdout);
	
	ch=getchar();
	while (ch>='0' && ch<='9')
	{
		n[++len]=ch-'0';
		ch=getchar();
	}
	fo(i,1,len/2)
	swap(n[i],n[len-i+1]);
	
	scanf("%d",&p);
	
	work2();
	N=pow(2,ceil(log(Mod+1)/log(2)));S=floor((log(N)+0.00001)/log(2));
	
	k=1;l=2;
	fo(i,1,S)
	{
		fo(j,0,k-1)
		{
			w[i][j].a=cos((2*M_PI)*j/l);
			w[i][j].b=sin((2*M_PI)*j/l);
		}
		k<<=1;l<<=1;
	}
	
	if (Mod<p-1)
	{
		work();
		
		fd(i,N-1,0)
		{
			j=floor(c[1][i].a+0.1);
			ans+=((j%p)>0);
		}
		
		fo(i,2,Len)
		ans=ans*(b[i]+1)%mod;
		
		printf("%d\n",ans);
	}
	else
	{
		++b[2];
		fo(i,2,Len)
		if (b[i]>=p)
		{
			++b[i+1];
			b[i]-=p;
		}
		if (b[Len+1]) ++Len;
		
		ans=1;
		fo(i,2,Len)
		ans=ans*(b[i]+1)%mod;
		
		printf("%d\n",ans);
	}
	
	fclose(stdin);
	fclose(stdout);
	
	return 0;
}