有一种组合方向的考虑有没有dalao肯高抬啊？
题目大意
有一个初始为空的数组$a$，按照以下的流程进行操作：
- 在$1cdots m$中等概率选出一个数$x$并添加到$a$的末尾
- 如果$a$中所有元素的$gcd=1$则完成这个数组$a$的修改
- 重复这一过程
求数组$a$的期望长度，$m le 10^5,mod 10^9+7$
题目分析
质因数分解的期望dp题
以下介绍的两个做法中，第一个做法本人不会所以 求助会做的dalao麻烦高抬一手 ；
第二个做法是对推得的dp式子质因数分解求解————当然网上更多的是莫比乌斯反演的做法，此处就不介绍了。
未完的做法一：组合考虑
考虑计算$f(d)$为：数组$a$的前$len-1$个（也即结束前的那一序列）$gcd=d$时，对答案的期望长度贡献。
那么枚举结束前的长度$i$，记$h(d)$为$1cdots m$中与$d$互质的数的个数；$g(i,d)$为一个$gcd=d$的长度为$i$的序列的概率，有$f(d)=sumlimits_{i=1}^{∞}(i+1) imes g(i,d) imes h(d)$。
当然，这样的式子远远不够。$h(d)$的计算非常轻松但是问题在于$g(i,d)$有没有什么组合上的表达方式。我最初以为$g(i,d)=lfloorfrac{m}{d} floor^{i-1} imes frac{1}{m}$，意即将序列看做是无序取入的集合，确保$i-1$个元素都是${d,2d,3dcdots}$，并钦点一个$d$.但是这个处理有两个问题：1.序列不能够这样看成无序集合；2.并不是必须要有一个$d$才能使数组的$gcd=d$，例如${2d,3d,5d}$就有$gcd=d$。
这个方向大概不太能行吧，为求出式子，$g(i,d)$可能还有更严苛的表达限制。
做法二：质因数分解
由于这里的长度是无穷的，所以处理角度还当是从序列的$gcd$入手。
记$f_d$为当前数组已经$gcd=d$，到达$gcd=1$的状态的期望步数。（这个状态要比我上面那个状态要优秀得多）
容易得到转移 $f_d=1+frac{1}{m} sumlimits_{i=1}^mf_{gcd(i,d)}$，
稍作化简得到 $f_d=1+frac{1}{m} sumlimits_{imid d,i e d}(f_i imes c(i,d)+lfloorfrac{m}{d} floor f_d)$，其中$c(i,d)=sumlimits_{j=1}^m[gcd (d,j)=i]$.
通过移项，有 $f_d=frac{1+frac{1}{m} sumlimits_{imid d,i
e d}f_i imes c(i,d)}{1-frac{1}{m}lfloorfrac{m}{d}
floor}$
那么对于这里$c(i,d)$的处理，就有很多方法了。这里介绍质因数分解的做法：
由于 $c(i,d)=sumlimits_{j=1}^m[gcd (d,j)=i]=sumlimits_{j=1}^{lfloor frac{m}{i} floor}[gcd (lfloorfrac{d}{i} floor,j)=1]$，考虑子问题：$1cdots c$中与$q$互质的数的个数。
那么将$q$唯一分解，有$q=p_1^{k_1}p_2^{k_2}cdots$。注意到实际上我们只关心$p_i$而不关心$k_i$，所求的也就是$1cdots c$中不存在任何$p_i$因数的数的个数。至此转为一个经典容斥问题，可以预处理出$r$个素因数再$2^r$容斥解决。
接下来是一些题外话，
上面提到的分解素因数的过程，理论上由于$10^5$内的数最多只有$6$个素因数，那么枚举每个素数检查是否为其因数的方法看上去也不算太差，
然而获得了1700ms+的“好”成绩：
一开始我还以为是$2^r$的dfs飞天了，然而后来意识到上述的这个素数分解的方法常数的确是巨大的……
可以这么看：对于一个要处理的数$x$，共约要遍历$frac{x}{ln x}$个素数。那么在这一部分上的总期望时间花费就大致是$frac{m^2sqrt m}{ln m}$（好像这个上界太松没什么价值？）。
总之判断约数是否为质数的方法就要快很多：
1 #include<bits/stdc++.h> 2 #define MO 1000000007 3 typedef long long ll; 4 typedef std::vector<int> vec; 5 const int maxn = 100035; 6 7 int n,num,cnt,pr[maxn],phi[maxn]; 8 bool vis[maxn]; 9 ll ans,f[maxn]; 10 vec fac[maxn],opt; 11 12 void makePrime() 13 { 14 phi[1] = 1; 15 for (int i=2; i<maxn; i++) 16 { 17 if (!vis[i]) pr[++pr[0]] = i, phi[i] = i-1; 18 for (int j=1,chk=1; j<=pr[0]&&i*pr[j]<maxn&&chk; j++) 19 if (pr[j]%i) phi[i*pr[j]] = phi[i]*(pr[j]-1); 20 else chk = 0, phi[i*pr[j]] = phi[i]*pr[j]; 21 } 22 for (int i=1; i<maxn; i++) 23 for (int j=i; j<maxn; j+=i) fac[j].push_back(i); 24 } 25 int qmi(int a, int b) 26 { 27 int ret = 1; 28 for (; b; b>>=1,a=1ll*a*a%MO) 29 if (b&1) ret = 1ll*ret*a%MO; 30 return ret; 31 } 32 void dfs(int x, int c, int t) 33 { 34 if (x==opt.size()) cnt += ((t&1)?-num/c:num/c); 35 else{ 36 dfs(x+1, c, t); 37 dfs(x+1, c*opt[x], t+1); 38 } 39 }41 int calc(int p, int d) 42 { 43 vec().swap(opt), p /= d, num = n/d; 44 for (int pri; p>1; ) 45 { 46 pri = fac[p][1]; 47 if (vis[pri]) break; 48 opt.push_back(pri); 49 while (p%pri==0) p /= pri; 50 } 51 cnt = 0, dfs(0, 1, 0); 52 return cnt; 53 } 54 int main() 55 { 56 scanf("%d",&n), ans = 1, f[1] = 0; 57 makePrime(); 58 for (int i=2; i<=n; i++) 59 { 60 int lwr = (1-1ll*qmi(n, MO-2)*(n/i)%MO+MO)%MO; 61 for (int j=1; j<fac[i].size()-1; j++) 62 f[i] = (f[i]+1ll*f[fac[i][j]]*calc(i, fac[i][j])%MO)%MO; 63 f[i] = 1ll*(1ll*f[i]*qmi(n, MO-2)%MO+1)*qmi(lwr, MO-2)%MO; 64 ans = (ans+f[i]+1)%MO; 65 } 66 printf("%lld ",1ll*ans*qmi(n, MO-2)%MO); 67 return 0; 68 }
END