「余姚中学 2019 联测 Day 4」随机除法

好题，难就难在转移的高位前缀和

首先是一个浅显的$\text{dp}$状态，令$n=\Pi prime_i^{c_i}$

则状态只跟$\{c_i\}$有关，这是一个可重集合，强制定义$c_i\ge c_{i-1}$最小表示出所有不同状态

搜索一下$\text{dp}$状态，发现只有$170000$左右的状态数

直接枚举因数转移复杂度显然是升天的，直接枚举子集状态转移复杂度也很高，并且不好确定系数

所以用一个高位前缀和处理优化枚举因数的转移

再说一遍，是高位前缀和枚举因数的转移，不同的因数可能对应同一个状态

常见的高位前缀和是$dp_{i,S}=dp_{i-1,S}+dp_{i,S-\{S_i\} }$

转移具有单调性，对于状态排序之后，定义辅助数组$dp_{i,j}$表示对于$i$这个状态

它的子集(注意这个子集是未排序的)中和它不同的最低位置$\ge j$的总和

计算高位前缀和时，每次转移只会对于一个位置改变

枚举状态时，取得位置是$j$，调用时需要排序

而排完序之后$j$可能会后移，所以需要定义成$\ge j$的，否则会算多

比如转移$(1,1,1)\leftarrow (0,1,1),(1,0,1),(1,1,0)$

如果定义成$\le j$的状态，三种状态转移之后都变成$(1,1,0)$

原先在这个状态里的三个位置编号是$(0,1,2)$

如果都去$(1,1,0)$这个状态里转移过来，原先$(0,1,2)$对应的下标位置改变，变成

$(1,2,0)$

$(0,2,1)$

$(0,1,2)$

我们访问的时候访问的应该是子状态中不同位置$\le j$的总和

而从这个下标改变的状态里转移过来时，原先$>j$的下标被移移动进$\le j$的范围

再转移就错了

所以正确定义状态之后就可以高位前缀和了

存储和访问状态可以用$\text{Hash,Trie,int128}$三种方法存储，$\text{int128}$真香啊

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
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef __int128 Node;
#define rep(i,a,b) for(int i=a,i##end=b;i<=i##end;++i)
#define drep(i,a,b) for(int i=a,i##end=b;i>=i##end;--i)
#define Mod1(x) ((x>=P)&&(x-=P))

const int N=180000,P=1e9+7;

ll qpow(ll x,ll k=P-2) {
	ll res=1;
	for(;k;k>>=1,x=x*x%P) if(k&1) res=res*x%P;
	return res; 
}
Node Max,st[N];
int m,a[100],cnt,dp[N][20],F[N],pri[]={0,2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79};
char str[30];

int Find(Node s){ return lower_bound(st+1,st+cnt+1,s)-st; }
void Load(Node s){
	rep(i,1,20) a[i]=0;
	for(int i=1;i<=20 && s>1; ++i) while(s%pri[i]==0) s=s/pri[i],a[i]++;
}
void dfs(int p,int lst,Node s) { 
	st[++cnt]=s;
	rep(i,1,lst) {
		if((s*=pri[p])>Max) return;
		dfs(p+1,i,s);
	}
}


int main(){
	freopen("div.in","r",stdin),freopen("div.out","w",stdout);
	Max=1e12,Max=Max*Max;
	dfs(1,100,1);
	fprintf(stderr,"States Number = %d\n",cnt);
	sort(st+1,st+cnt+1);
	rep(i,2,cnt) {
		Node now; Load(now=st[i]);
		int c=1;
		rep(i,1,18) c=1ll*c*(a[i]+1)%P;
		drep(j,18,0) {
			if(!a[j]) dp[i][j]=dp[i][j+1];
			else {
				int k=j;
				while(k>1 && a[k-1]==a[k]) --k;
				int p=Find(now/pri[j]);
				drep(d,j,k) {
					dp[i][d]=dp[i][d+1]+dp[p][d];
					Mod1(dp[i][d]);
				}
				j=k;
			}
		}
		F[i]=(dp[i][1]+c)*qpow(c-1)%P;
		rep(j,1,18) dp[i][j]+=F[i],Mod1(dp[i][j]);
	}
	while(~scanf("%s%d",str,&m)) {
		Node n=0;
		for(int i=0;str[i];++i) n=n*10+(str[i]-'0');
		rep(i,1,m) {
			int x;scanf("%d",&x); a[i]=0;
			while(n%x==0) n=n/x,a[i]++;
		}
		sort(a+1,a+m+1,greater <int> ());
		n=1;
		rep(i,1,m) rep(j,1,a[i]) n=n*pri[i];
		printf("%d\n",F[Find(n)]);
	}
}