莫比乌斯反演（洛谷P5221）

最新推荐文章于 2024-02-11 14:28:55 发布

原创最新推荐文章于 2024-02-11 14:28:55 发布 · 348 阅读

0 ·

CC 4.0 BY-SA版权

ACM 专栏收录该内容

30 篇文章

订阅专栏

本文深入解析了一道数学算法竞赛题目，针对大规模数论问题，通过巧妙的数学公式推导，将复杂度极高的原问题转化为易于处理的形式。详细介绍了如何利用莫比乌斯反演和欧拉函数解决∏i=1n∏j=1nlcm(i,j)gcd(i,j)(mod104857601)的问题，并提供了高效的C++代码实现。

题意：
给一个n（1<=n<=1000000）;
求
$∏i=1n∏j=1nlcm(i,j)gcd(i,j)(mod104857601)\prod_{i=1}^n\prod_{j=1}^n\frac{lcm(i,j)}{gcd(i,j)}(mod 104857601)$
题解：
公式推演：
$∏i=1n∏j=1nlcm(i,j)gcd(i,j)\prod_{i=1}^n\prod_{j=1}^n\frac{lcm(i,j)}{gcd(i,j)}$

$=∏i=1n∏j=1ni∗jgcd(i,j)2=\prod_{i=1}^n\prod_{j=1}^n\frac{i*j}{gcd(i,j)^2}$

$=(∏i=1n∏j=1ni∗j)(∏i=1n∏j=1ngcd(i,j)−2)=(\prod_{i=1}^n\prod_{j=1}^ni*j)(\prod_{i=1}^n\prod_{j=1}^ngcd(i,j)^{-2})$

$=(∏i=1nin∗n!)(∏i=1n∏j=1ngcd(i,j)−2)=(\prod_{i=1}^ni^n*n!)(\prod_{i=1}^n\prod_{j=1}^ngcd(i,j)^{-2})$

$=(n!)2n∗(∏i=1n∏j=1ngcd(i,j)−2)=(n!)^{2n}*(\prod_{i=1}^n\prod_{j=1}^ngcd(i,j)^{-2})$

把 $∏i=1n∏j=1ngcd(i,j)\prod_{i=1}^n\prod_{j=1}^ngcd(i,j)$ 拿出来单独看

$∏i=1n∏j=1ngcd(i,j)\prod_{i=1}^n\prod_{j=1}^ngcd(i,j)$

$=∏d=1n∏i=1n∏j=1n[gcd(i,j)==d]=\prod_{d=1}^n\prod_{i=1}^n\prod_{j=1}^n[gcd(i,j)==d]$

$=∏d=1nd∑i=1n∑j=1n[gcd(i,j)==d]=\prod_{d=1}^nd^{\sum_{i=1}^n\sum_{j=1}^n[gcd(i,j)==d]}$

$=∏d=1nd∑i=1⌊nd⌋∑j=1⌊nd⌋[gcd(i,j)==1]=\prod_{d=1}^nd^{\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor}[gcd(i,j)==1]}$

令 $sum[x]=∑i=1xϕ(i)sum[x]=\sum_{i=1}^x\phi(i)$

由莫比乌斯反演可以得到：
$(n!)2n∗(∏d=1nd2sum⌊nd⌋−1)−2(n!)^{2n}*(\prod_{d=1}^nd^{2sum\lfloor\frac{n}{d}\rfloor-1})^{-2}$

边上那个-2用欧拉降幂整掉，这题只给7MB空间，很容易被卡，用原地做法，分块也可以，分块确实可以更快点，但是得多开点空间，我这里就不用分块了
code：

#include<bits/stdc++.h>
#define ll long long
using namespace std;
const int N=1e6+5;
const int mod=104857601;
int prime[80000],phi[N],tot=0;
ll fac=1;
inline int read() {
    int x=0,w=0;char ch=getchar();
    while(!isdigit(ch))w|=ch=='-',ch=getchar();
    while(isdigit(ch))x=(x<<1)+(x<<3)+(ch^48),ch=getchar();
    return w?-x:x;
}
void init(int maxn){
    phi[1]=1;
    for(int i=2;i<maxn;i++){
        fac=fac*i%mod;
        if(!phi[i]){
            prime[tot++]=i;
            phi[i]=i-1;
        }
        for(int j=0;j<tot&&prime[j]*i<maxn;j++){
            if(i%prime[j]==0){
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }
            phi[i*prime[j]]=phi[i]*(prime[j]-1);
        }
    }
    for(int i=1;i<maxn;i++)phi[i]=(phi[i-1]+phi[i])%(mod-1);
}
ll fast(ll x,ll y=mod-1){
    ll ans=1;
    while(y){
        if(y&1)ans=ans*x%mod;
        x=x*x%mod;
        y>>=1;
    }
    return ans;
}
int main()
{
    int n=read();
    init(n+1);
    ll ans1,ans2=1;
    ans1=fast(fac,2*n);
    for(int i=1;i<=n;i++){
        ll power=((2-4*phi[n/i])%(mod-1)+mod-1)%(mod-1);
        ans2=ans2*fast(i,power)%mod;
    }
    printf("%lld",ans1*ans2%mod);
    return 0;
}