P5106 dkw的lcm

传送门

终于A了……细节真多……

首先我们发现这是个连乘,而且(phi)是个积性函数,所以我们可以考虑不同的质因子以及它的不同次数的贡献。简单来说就是把每一次的(phi(lcm(i_1,i_2,...)))拆成一堆(phi(p^c))的乘积

如果枚举到的(k)个数里质因子(p)的最大次数为(c),那么最终的答案就要乘上一个(phi(p^c))。所以我们的目的就是要求出这(k)个数里(p)的最大次数为(c)时的方案数(t),然后这部分对答案的贡献就是(phi(p^c)^t)

计算方案的部分可以用容斥解决。设(d_p(i))表示质因子(p)的最高次数为(i)的总方案数,(s_p(i))表示质因子(p)的最高次数不超过(i)的总方案数,那么(d_p(i)=s_p(i)-sum_{j=1}^{i-1}d_p(j)),边界条件为(d_p(0)=s_p(0)),前缀和优化一下就能快速计算了

最后是(s_p(i))怎么计算。我们可以单独考虑每一位的方案数最后(k)位乘起来就好了。先考虑质因子(p)的最高次数刚好为(i)的情况,对于某一位来说,能选的数(j)要满足(j imes p^ileq n)(j)不包含(p)这个质因子,那么满足条件的(j)的次数就是(leftlfloorfrac{n}{p^i} ight floor-leftlfloorfrac{n}{p^{i+1}} ight floor),于是前缀和一下就能得到对一个数来说的方案数为(n-leftlfloorfrac{n}{p^i} ight floor),那么(s_p(i)=(n-leftlfloorfrac{n}{p^i} ight floor)^k)

然后注意一个细节,因为上面(s_p)(d_p)的计算基本都是要炸精度所以要取模,然而因为这两个东西是作为次数出现的,所以根据欧拉定理取模的时候要模(phi(P))而不是(P),因为这个东西调死掉……

然后就没有然后了

//minamoto
#include<bits/stdc++.h>
#define R register
#define ll long long
#define inf 0x3f3f3f3f
#define fp(i,a,b) for(R int i=a,I=b+1;i<I;++i)
#define fd(i,a,b) for(R int i=a,I=b-1;i>I;--i)
#define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v)
using namespace std;
const int N=2e6+5,P=1e9+7,Phi=P-1;
int p[N],vis[N],phi[N],kkk[N];bool flag=0;
int n,k,m,ans;ll sum,t;
inline int add(R int x,R int y,R int P){return x+y>=P?x+y-P:x+y;}
inline int dec(R int x,R int y,R int P){return x-y<0?x-y+P:x-y;}
inline int ksm(R int x,R int y,R int P){
	int res=1;
	for(;y;y>>=1,x=1ll*x*x%P)if(y&1)res=1ll*res*x%P;
	return res;
}
void init(int n){
	phi[1]=kkk[1]=1;
	fp(i,2,n){
		if(!vis[i])p[++m]=i,phi[i]=i-1,kkk[i]=ksm(i,k,Phi);
        //预处理一下phi和所有数的k次幂
		for(R int j=1;j<=m&&1ll*i*p[j]<=n;++j){
			vis[i*p[j]]=1,kkk[i*p[j]]=1ll*kkk[i]*kkk[p[j]]%Phi;
			if(i%p[j]==0){phi[i*p[j]]=phi[i]*p[j];break;}
			phi[i*p[j]]=phi[i]*(p[j]-1);
		}
	}
}
void solve(int p){
	sum=kkk[n-n/p];
	for(R int i=p;i<=n;(1ll*i*p<=n?i*=p:i=inf)){
		t=dec(kkk[n-n/(1ll*i*p)],sum,Phi);
		ans=1ll*ans*ksm(phi[i],t+Phi,P)%P,
		sum=add(sum,t,Phi);
	}
}
int main(){
//	freopen("testdata.in","r",stdin);
	scanf("%d%d",&n,&k),ans=1,init(n);if(n==1)return puts("1"),0;
	fp(i,1,m)solve(p[i]);return printf("%d
",ans),0;
}
原文地址:https://www.cnblogs.com/bztMinamoto/p/10136447.html