题解「Luogu5221 Product」

题意

求这个东西:

[prod_{i=1}^Nprod_{j=1}^Nfrac{{ m{lcm}}(i,j)}{{ m{gcd}}(i,j)} ({ m{mod}} 104857601) ]

题解

根据

[{ m{lcm}}(i,j)=frac{i imes j}{{ m{gcd}}(i,j)} ]

化简式子:

[egin{aligned} ext{Ans}&=prod_{i=1}^Nprod_{j=1}^Nfrac{i imes j}{{ m{gcd}}^2(i,j)}\ &=frac{prod_{i=1}^Nprod_{j=1}^Nij}{prod_{i=1}^Nprod_{j=1}^N{ m{gcd}}^2(i,j)} end{aligned} ]

先看位于分子的式子:

[egin{aligned} prod_{i=1}^Nprod_{j=1}^Nij&=prod_{i=1}^Ni^N (N!)\ &=(N!)^N(N!)^N\ &=(N!)^{2N} end{aligned} ]

快速幂求即可。

然后是分母的式子:

先不考虑 ({ m{gcd}}(i,j)) 的次数。枚举 ({ m{gcd}}(i,j)) ,若有数对 (i,j) 满足 ({ m{gcd}}(i,j)=d) ,那么 (d) 对答案会产生贡献。有:

[prod_{i=1}^Nprod_{j=1}^N{ m{gcd}}(i,j)=prod_{d=1}^Nd^{f(d)} ]

[egin{aligned} f(d)&=sum_{i=1}^Nsum_{j=1}^N[{ m{gcd}}(i,j)=d]\ &=sum_{i=1}^{lfloorfrac{N}{d} floor}sum_{j=1}^{lfloorfrac{N}{d} floor}sum_{k|i,k|j}mu(k)\ &=sum_{k=1}^{lfloorfrac{N}{d} floor}mu(k)lfloorfrac{N}{kd} floor^2 end{aligned} ]

枚举 (d) ,整除分块求 (f(d)) 即可。


( ext{Code}:)

#include <iostream>
#include <cstring>
#include <cstdio>
#include <algorithm>
#include <cmath>
#define maxn 1000005
#define Rint register int
#define INF 0x3f3f3f3f
using namespace std;
typedef long long lxl;
const lxl mod=104857601;

template <typename T>
inline T read()
{
	T x=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9') {if(ch=='-') f=-1;ch=getchar();}
	while(ch>='0'&&ch<='9') {x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
	return x*f;
}

int prime[maxn],cnt;
bool flag[maxn];
int mu[maxn];

inline void sieve()
{
	mu[1]=1;
	for(int i=2;i<maxn;++i)
	{
		if(!flag[i]) prime[++cnt]=i,mu[i]=-1;
		for(int j=1;j<=cnt&&i*prime[j]<maxn;++j)
		{
			flag[i*prime[j]]=true;
			if(!(i%prime[j])) break;
			mu[i*prime[j]]=-mu[i];
		}
	}
	for(int i=1;i<maxn;++i)
		mu[i]+=mu[i-1];
}

inline lxl fmi(lxl a,lxl b)
{
	lxl ans=1;
	a%=mod;
	while(b>0)
	{
		if(b&1) ans=(ans*a)%mod;
		a=(a*a)%mod;
		b>>=1;
	}
	return ans;
}

inline lxl calcu(lxl N)
{
	lxl res=0;
	for(lxl l=1,r=0;l<=N;l=r+1)
	{
		r=N/(N/l);
		res+=(mu[r]-mu[l-1])*(N/l)*(N/l);
	}
	return res;
}

inline lxl solve(int N)
{
	lxl res=1;
	for(lxl i=1;i<=N;++i) res=(res*i)%mod;
	res=fmi(res,2*N);
	lxl flr=1;
	for(lxl d=1;d<=N;++d)
		flr=flr*fmi(d,calcu(N/d))%mod;
	flr=fmi(flr*flr%mod,mod-2);
	return res*flr%mod;
}

int main()
{
	// freopen("P5221.in","r",stdin);
	sieve();
	int N=read<int >();
	printf("%lld
",solve(N));
	return 0;
}

原文地址:https://www.cnblogs.com/syc233/p/13562004.html