HDU 5942 Just a Math Problem（莫比乌斯反演+容斥原理）

Description

已知$f(k)$为$k$的素因子个数，$g(k)=2^{f(k)}$，给一正整数$n$，求$\sum\limits_{i=1}^ng(i)$

Input

第一行一整数$T$表示用例组数，每组用例输入一整数$n(1\le T\le 50,1\le n\le 10^{12})$

Output

输出$\sum\limits_{i=1}^ng(i)$，结果模$10^9+7$

Sample Input

3
1
10
100

Sample Output

Case #1: 1
Case #2: 23
Case #3: 359

Solution1

令$k=p_1^{a_1}...p_m^{a_m}$，则$g(k)=2^{f(k)}=2^m$，即为集合$\{p_1^{a_1},...,p_m^{a_m}\}$的子集个数，该集合一个子集的乘积$p$与该子集的补集的乘积$q$显然互素，且该对应是一一的，故$g(k)=|\{(p,q)=1,pq=k\}|$

对于$h(n)=\sum\limits_{i=1}^ng(i)$，同样的有$h(n)=|\{(p,q)=1,pq\le n\}|=2|\{(p,q)=1,pq\le n,p<q\}|+1$

由于$pq\le n,p<q$，故$p\le \sqrt{n}$，进而有$|\{(p,q)=1,pq\le n,p<q\}|=\sum\limits_{p=1}^{\sqrt{n}}(Count(n/p,p)-\varphi(p))$，其中$Count(m,p)$为$1$~$m$中与$p$互素的数的个数，用容斥原理即可求出

Code1

#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<vector>
#include<queue>
#include<map>
#include<set>
#include<ctime>
using namespace std;
typedef long long ll;
#define maxn 1000001
#define mod 1000000007ll
int euler[maxn],prime[maxn],res;
void get_euler(int n)//求n以内所有数的欧拉函数值,顺便得到n以内所有素数共res个
{
    memset(euler,0,sizeof(euler));
    euler[1]=1;
    res=0;
    for(int i=2;i<=n;i++)
    {
        if(!euler[i])euler[i]=i-1,prime[res++]=i;
        for(int j=0;j<res&&prime[j]*i<=n;j++)
        {
            if(i%prime[j]) euler[prime[j]*i]=euler[i]*(prime[j]-1);
            else
            {
                euler[prime[j]*i]=euler[i]*prime[j];
                break;
            }
        }
    }
}
int factor[maxn][15],cnt[maxn];
void get_factor(int n)
{
    cnt[n]=0;
    int p=n;
    for(int i=2;i*i<=p;i++)
        if(p%i==0)
        {
            factor[n][cnt[n]++]=i;
            while(p%i==0)p/=i;    
        }
    if(p>1)factor[n][cnt[n]++]=p;
}
ll count(ll n,int p)//|{x|(x,p)=1,1<=x<=n}|
{
    ll ans=n;
    if(cnt[p]==-1)get_factor(p);
    int N=1<<cnt[p];
    for(int i=1;i<N;i++)
    {
        int num=0;
        ll temp=1;
        for(int j=0;j<cnt[p];j++)
            if(i&(1<<j))num++,temp*=factor[p][j];
        if(num&1)ans-=n/temp;
        else ans+=n/temp;
    }
    return ans%mod;
}
ll Solve(ll n)
{
    ll ans=0;
    for(int p=1;1ll*p*p<n;p++)ans=((ans+count(n/p,p)-euler[p])%mod+mod)%mod;
    return (ans*2+1)%mod;
}
int main()
{
    get_euler(1000000);
    memset(cnt,-1,sizeof(cnt));
    int T,Case=1;
    scanf("%d",&T);
    while(T--)
    {
        ll n;
        scanf("%lld",&n);
        printf("Case #%d: %lld\n",Case++,Solve(n));
    }    
    return 0;
}

Solution2

$g(i)=|\{(p,q)=1,pq=i\}|=\sum\limits_{d|i}[(d,\frac{i}{d})=1]=\sum\limits_{d|i}\sum\limits_{k|(d,\frac{i}{d})}\mu(k)$

记$d=sk,i=stk^2$，$h(n)=\sum\limits_{i=1}\sum\limits_{d|i}\sum\limits_{k|(d,\frac{i}{d})}\mu(k)=\sum\limits_{k=1}^{\sqrt{n}}\mu(k)\sum\limits_{s=1}^{\lfloor\frac{n}{k^2}\rfloor}\lfloor\frac{n}{k^2i}\rfloor$

记$S(n)=\sum\limits_{i=1}^n\lfloor\frac{n}{i}\rfloor$，则$h(n)=\sum\limits_{k=1}^{\sqrt{n}}\mu(k)S(\lfloor\frac{n}{k^2}\rfloor)$

对于较小的$n$预处理$S(n)$，对于较大的 n n <script type="math/tex" id="MathJax-Element-88">n</script>分块求

Code2

#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<vector>
#include<queue>
#include<map>
#include<set>
#include<ctime>
using namespace std;
typedef long long ll;
#define maxn 1000001
#define mod 1000000007ll
bool check[maxn];
int prime[maxn],mu[maxn];
ll S[maxn];
void Moblus(int n)
{
    memset(check,0,sizeof(check));
    mu[1]=1;
    int tot=0;
    for(int i=2;i<=n;i++)
    {
        if(!check[i])
        {
            prime[tot++]=i;
            mu[i]=-1;
        }
        for(int j=0;j<tot;j++)
        {
            if(i*prime[j]>n)break;
            check[i*prime[j]]=1;
            if(i%prime[j]==0)
            {
                mu[i*prime[j]]=0;
                break;
            }
            mu[i*prime[j]]=-mu[i];
        }
    }
}
ll Sum(ll n)
{
    if(n<maxn&&S[n]!=-1)return S[n];
    ll ans=0;
    for(ll i=1,pre;i<=n;i=pre+1)
    {
        pre=n/(n/i);
        ans=(ans+1ll*(pre-i+1)*(n/i)%mod)%mod;
    }
    if(n<maxn)S[n]=ans;
    return ans;
}
int main()
{
    Moblus(1000000);
    memset(S,-1,sizeof(S));
    int T,Case=1;
    scanf("%d",&T);
    while(T--)
    {
        ll n,ans=0;
        scanf("%lld",&n);
        for(int i=1;1ll*i*i<=n;i++)
            if(mu[i])ans+=mu[i]*Sum(n/i/i),ans=(ans+mod)%mod;
        printf("Case #%d: %lld\n",Case++,ans);
    }    
    return 0;
}