Function [容斥原理]

链接

2020 年百度之星·程序设计大赛 - 初赛一
HDOJ6750
http://acm.hdu.edu.cn/showproblem.php?pid=6750

题意

令$f(n)=\sum _{i|n,gcd(i,\frac{n}{i})=1}i$
求$S(n)=f(1)+f(2)+...+f(n)$
$n\le 10^{12}，1\le test\le 10$

分析

首先对 $S(n)$ 进行化简：
$S(n)=-1+\sum _{(i,j)}i+j$
其中$i,j$满足$1\le i\le j\le n$ 且 $i\ast j\le n$ 且 $gcd(i,j)=1$

为什么会有一个$-1$呢？因为在 $i=1,j=1$ 的时候 $1$ 被计算了两遍
接下来我们可以发现，上面的式子中的 $i$ 满足 $i\le \sqrt{n}$

若$i>1$，那么 $j$ 满足 $i<j\le \frac{n}{i}$ ，因此对于 $\sqrt{n}$ 以内的 $i$
要求出 $[i+1,\frac{n}{i}]$ 内与 $i$ 互质的 $j$ 的个数以及总和

由于 $gcd(i,j)=1$ 等价于 $gcd(i,j-i)=1$ ，那么只需要算出
$k\in [1,\frac{n}{i}-i]$ 内与 $i$ 互质的 $k$ 的个数以及总和，就能算出 $j$ 的总和

对 $i$ 进行质因数分解，上述 $k$ 的个数以及总和都能用容斥原理算出。

代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef double db;
template<class T>inline void MAX(T &x,T y){if(y>x)x=y;}
template<class T>inline void MIN(T &x,T y){if(y<x)x=y;}
template<class T>inline void rd(T &x){
    x=0;char o,f=1;
    while(o=getchar(),o<48)if(o==45)f=-f;
    do x=(x<<3)+(x<<1)+(o^48);
    while(o=getchar(),o>47);
    x*=f;
}
const int M=1e6+5;
const int P=1e9+7;
int cas,prime[M/10],ptot,pre[M];
void calc(int *Q,int top,ll n,int &sum,int &cnt){
    for(int i=0;i<(1<<top);i++){
        int cur=0,num=1;
        for(int j=0;j<top;j++)if(i&1<<j)cur++,num*=Q[j];
        int tmp=n/num%P;
        if(cur%2==0)sum=(sum+1ll*tmp*(tmp+1)/2%P*num)%P,cnt=(cnt+tmp)%P;
        else sum=(sum-1ll*tmp*(tmp+1)/2%P*num+P)%P,cnt=(cnt-tmp+P)%P;
    }
}
static int Q[35];
int main(){
#ifndef ONLINE_JUDGE
    freopen("jiedai.in","r",stdin);
//    freopen("jiedai.out","w",stdout);
#endif
    for(int i=2;i<M;i++){
        if(pre[i]==0)pre[i]=prime[++ptot]=i;
        for(int j=1;j<=ptot;j++){
            int t=i*prime[j];
            if(t>=M)break;
            pre[t]=prime[j];
            if(i%prime[j]==0)break;
        }
    }
    rd(cas);
    while(cas--){
        ll n;
        rd(n);
        int ans=(n-1+(n%P)*(n%P+1)/2)%P;
        for(int i=2;1ll*i*i<=n;i++){
            int sum=0,cnt=0,top=0,num=i;
            while(num>1){
                int x=pre[num];
                Q[top++]=x;
                while(num%x==0)num/=x;
            }
            calc(Q,top,n/i-i,sum,cnt);
            ans=(ans+2ll*cnt*i%P+sum)%P;
        }
        printf("%d\n",ans);
    }
    return (0-0);
}

比赛的时候少敲了一个1ll*，交上去居然是MLE，赛后加上就过了，呜呜 。