首页 > 代码库 > bzoj 4961: 除除除
bzoj 4961: 除除除
Description
我们定义一种操作是将一个正整数n(n>1)用某个整数d替换，这个数必须是n的约数(1≤d≤n)。给你一个正整数n，
你需要确定操作进行的期望次数，如果我们希望不断地使用这种操作来将n变成1，假设每次操作选择每个可能的d
的概率均等。为了便于计算，输入将给出n和它的所有不同质因数p_1,p_2,?p_m，保证n恰好有m个不同的质因数。
为了便于输出，设答案是有理数a/b，并且有bc≡1(mod10^9+7)，你只需要输出ac对10^9+7取模的值。例如当n=351
384000时，期望运算的次数为
1384855049944986283970053414177036273994739277918823/282971529872677632598150446595770345000925504317000≈4.893973081
但你只需要输出321468106即可。
Input
输入包含多组测试数据，以EOF结束。
对于每组测试数据：
第一行包含两个正整数n和m，其中m表示n的不同质因数个数，满足2≤n≤10^24。
第二行包含m个质数p_1,p_2,...,p_m，对于i=1,2,...,m满足2≤p_i≤10^6。
约200000组数据。
Output
对于每组测试数据，输出一行一个整数，表示题目要求输出的值。
对一个询问，答案只和每个质因子的幂次构成的可重集有关，在数据范围内只有约170000个本质不同的询问，因此可以预处理递推出答案
递推式为（f为答案，d为约数个数）：
$f(n)=\frac{d(n)+\sum_{i|n and i<n}f(d)}{d(n)-1}$
为了优化这个递推式，将n的质因子的幂次降序排列为p(n,1),p(n,2),p(n,3),...，求和部分考虑记录一个前缀和g(n,k)辅助计算，表示满足 当x<=k时p(a,x)=p(n,x)，否则p(a,x)<=p(n,x) 的f(a)之和，按p(n)的字典序升序处理，适当用hash存储g可以使转移复杂度与g的状态数（约1300000）同阶
#include<cstdio>typedef unsigned int u32;typedef unsigned long long u64;u64 pp[107];const u64 _eq=1844677;const u32 P=1e9+7;struct num{ u64 v0,v1,v2; void push(u32 x){ v0=v0*10+x,v1*=10,v2*=10; v1+=v0>>32,v0=u32(v0); v2+=v1>>32,v1=u32(v1); } bool div(u32 x){ u64 a0=v0,a1=v1,a2=v2; a1+=a2%x<<32,a2/=x; a0+=a1%x<<32,a1/=x; if(a0%x)return 0; v0=a0/x,v1=a1,v2=a2; return 1; } bool read(){ v0=v1=v2=0; int c=getchar(); while(c<48){ if(c<0)return 0; c=getchar(); } while(c>47)push(c-48),c=getchar(); return 1; }};int _(){ int x=0,c=getchar(); while(c<48)c=getchar(); while(c>47)x=x*10+c-48,c=getchar(); return x;}const u32 M=1<<22;struct hmp{ u64 hx[M]; int hy[M],y0; int&operator[](u64 x){ if(!x)return y0; int w=x&(M-1); while(hx[w]){ if(hx[w]==x)return hy[w]; w=(w+1237)&(M-1); } hx[w]=x; return hy[w]; }}H;int ps[]={2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73};int ss[33];u32 iv[100007],va=0;u32 inv(u32 a){ if(a<=100000)return iv[a]; u32 v=1; for(u32 n=P-2;n;n>>=1,a=u64(a)*a%P)if(n&1)v=u64(v)*a%P; return v;}void dfs(int w,int d,double s,u64 h){ if(w){ u64 h1=h,s0=0; u32 c=1; for(int i=0;i<w;++i){ int x=ss[i]; c*=(x+1); s0+=H[h1+pp[x-1]-pp[x]]; h1+=(_eq-1)*pp[x]; } va+=H[h1]=(s0+c)%P*inv(c-1)%P; for(int i=w-1;i>=0;--i){ int x=ss[i]; u64 h2=h1+(1-_eq)*pp[x]; H[h2]=(H[h1]+H[h2+pp[x-1]-pp[x]])%P; h1=h2; } } for(int i=1;i<=d;++i){ s*=ps[w]; if(s>1.01e24)return; ss[w]=i; dfs(w+1,i,s,h+pp[i]); }}int main(){ num x; u32 m; pp[1]=1; iv[1]=1; for(int i=2;i<=100000;++i)iv[i]=u64(P-P/i)*iv[P%i]%P; for(u32 i=2;i<=100;++i)pp[i]=pp[i-1]*149; dfs(0,100,1,0); while(x.read()){ u64 h=0; for(m=_();m;--m){ u32 y=_(),t=0; while(x.div(y))++t; h+=pp[t]; } printf("%d\n",H[h*_eq]); } return 0;}
bzoj 4961: 除除除
声明：以上内容来自用户投稿及互联网公开渠道收集整理发布，本网站不拥有所有权，未作人工编辑处理，也不承担相关法律责任，若内容有误或涉及侵权可进行投诉： 投诉/举报 工作人员会在5个工作日内联系你，一经查实，本站将立刻删除涉嫌侵权内容。