[数论]HDU 6750 Function 百度之星2020初赛第一场H题

HDU 6750 Function

百度之星2020初赛第一场H题

题目大意

设 (f(n)=sum_{d|n}d[gcd(d,frac{n}{d})=1]),
求 (S(n)=sum_{i=1}^{n}f(i))
(1leq nleq 10^{12})

题解

不难发现 (f(n)) 是积性函数，证明如下：
令 (gcd(n,m)=1),
(f(n)=sum_{d|n}d[gcd(d,frac{n}{d})=1]),
(f(m)=sum_{d|m}d[gcd(d,frac{m}{d})=1]),
(f(nm)=sum_{d|nm}d[gcd(d,frac{nm}{d})=1]),
(f(n)f(m)=left(sum_{d|n}d[gcd(d,frac{n}{d})=1] ight)left(sum_{d|m}d[gcd(d,frac{m}{d})=1] ight))
(=sum_{d|n}sum_{p|m}dp[gcd(d,frac{n}{d})=1][gcd(p,frac{m}{p})=1])
(=sum_{d|n}sum_{p|m}dp[gcd(dp,frac{nm}{dp})=1])
(=sum_{dp|nm}dp[gcd(dp,frac{nm}{dp})=1])
(=sum_{d|nm}d[gcd(d,frac{nm}{d})=1])
(=f(nm))
所以 (f(n)) 是积性函数，我们要求的 (S(n)) 是一个积性函数的前缀和。

考虑构造两个积性函数 (g) 和 (h),使得 (f=g*h),
即 (f(n)=sum_{d|n}g(frac{n}{d})h(d))。
那么 (S(n)=sum_{i=1}^{n}f(i)=sum_{i=1}^{n}sum_{d|i}g(frac{i}{d})h(d)=sum_{d=1}^{n}sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}g(i)h(d))

因为 (n) 的范围非常大，式子中的 (h(d)) 不太好处理。我们期望降低 (d) 的枚举范围。

考虑当 (n=p,pin Prime) 的情况,
则 (f(p)=sum_{d|p}g(d)h(frac{p}{d})=g(1)h(p)+g(p)h(1)=h(p)+g(p)),
又由 (f(n)=sum_{d|n}d[gcd(d,frac{n}{d})=1]) 得, (f(p)=1+p)
我们希望(forall pin Prime,h(p)=0), 此时(f(p)=g(p)=1+p)。
发现 (sigma_1(n)=sum_{d|n}d, sigma_1(p)=1+p)，所以我们不妨令 (g(n)=sigma_1(n))。
现在的问题是 (h(n)) 是什么。

再考虑当 (n=p^2,pin Prime) 的情况,
(f(p^2)=sum_{d|p^2}sigma_1(d)h(frac{p^2}{d})=sigma_1(1)h(p^2)+sigma_1(p)h(p)+sigma_1(p^2)h(1))
(=h(p^2)+sigma_1(p^2)=h(p^2)+1+p+p^2)
又由于 (f(p^2)=1+p^2)， 所以 (h(p^2)=-p)

再考虑当 (n=p^e,pin Prime,e>2) 的情况。
(f(p^e)=sum_{d|p^e}sigma_1(d)h(frac{p^e}{d})=sigma_1(p^0)h(p^e)+sigma_1(p)h(p^{e-1})+dots+sigma_1(p^{e-2})h(p^2)+sigma_1(p^{e-1})h(p^1)+sigma_1(p^e)h(p^0))
(=sigma_1(p^0)h(p^e)+sigma_1(p)h(p^{e-1})+dots+sigma_1(p^{e-3})h(p^3)-psigma_1(p^{e-2})+sigma_1(p^e))
(=sigma_1(p^0)h(p^e)+sigma_1(p)h(p^{e-1})+dots+sigma_1(p^{e-3})h(p^3)-p(1+p+dots+p^{e-2})+(1+p+dots+p^e))
(=sigma_1(p^0)h(p^e)+sigma_1(p)h(p^{e-1})+dots+sigma_1(p^{e-3})h(p^3)+1+p^e)
又因为 (f(p^e)=1+p^e)
所以当 (pin Prime,e>2) 时，(h(p^e)=0)

综上 (h(1)=1,h(p)=0,h(p^2)=-p,h(p^e)=0,e>2)
不妨令 (t(n)=h(n^2))
发现 (h(n)) 对素数的平方比较特殊，考虑到到莫比乌斯函数 (mu(n)) 的性质，得 (h(p^2)=mu(p)p=-p)
考虑到 (h(n)) 也是积性函数，
则 (h(n)=h(p_1^{k_1}p_2^{k_2}dots p_m^{k_m})=h(p_1^{k_1})h(p_2^{k_2})dots h(p_m^{k_m}))
根据前面的性质，只有当n的所有质因子的质数 (k_i) 均为2时，(h(n)) 才不为0。所以对于不为0的 (h(n)),(n)的最大质因子一定小于等于 (sqrt n)。
所以 (S(n)=sum_{d=1}^{n}sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}sigma_1(i)h(d)=sum_{d=1}^{sqrt n} h(d^2)sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}sigma_1(i)=sum_{d=1}^{sqrt n} mu(d)dsum_{i=1}^{leftlfloorfrac{n}{d^2} ight floor}sigma_1(i))
显然对于后面的 (sum_{i=1}^{leftlfloorfrac{n}{d^2} ight floor}sigma_1(i)) 我们又可以整除分块来做。
约数个数和函数 (sigma_1(n)) 的前缀和同样也可以由 (sum_{i=1}^{n}sigma_1(n)=sum_{i=1}^{n}sum_{d|i}d=sum_{d=1}^{n}dsum_{i=1}^{leftlfloorfrac{n}{d} ight floor}1=sum_{d=1}^{n}dleftlfloorfrac{n}{d} ight floor) 使用整除分块求出。
官方题解指出时间复杂度为 (O(sqrt n log n))，但我并不会证。

Code

#include <iostream>
#include <algorithm>
#include <cstring>
#include <string>
#include <cstdio>
#include <vector>
#include <cmath>
using namespace std;

#define RG register int
#define LL long long
const LL MOD=1000000007LL;
const LL Inv2=500000004LL;

template<typename elemType>
inline void Read(elemType &T){
    elemType X=0,w=0; char ch=0;
    while(!isdigit(ch)) {w|=ch=='-';ch=getchar();}
    while(isdigit(ch)) X=(X<<3)+(X<<1)+(ch^48),ch=getchar();
    T=(w?-X:X);
}

const int maxn=1000000;
vector<int> Prime;
int Mu[maxn+5];
LL Sum[maxn+5];
bool not_Prime[maxn+5];

void Get_Mu(int Len){
    not_Prime[1]=1;
    Mu[1]=1;
    for(RG i=2;i<=Len;i++) {
        if(!not_Prime[i]){
            Prime.push_back(i);Mu[i]=-1;
        }
        for(RG j=0;j<(int)Prime.size();++j) {
            LL mid=(LL)Prime[j]*i;
            if(mid>Len) break;
            not_Prime[mid]=1;
            if(i%Prime[j]==0) {Mu[mid]=0;break;}
            Mu[mid]=-Mu[i];
        }
    }
    for(RG i=1;i<=maxn;++i)
        Sum[i]=((Sum[i-1]+(LL)Mu[i]*i%MOD)%MOD+MOD)%MOD;
    return;
}

LL Sigma1(LL n){
    LL Res=0;
    for(LL L=1,R;L<=n;L=R+1){
        R=n/(n/L);
        LL temp=(L+R)%MOD*((R-L+1)%MOD)%MOD*Inv2%MOD;
        Res=((Res+temp*(n/L)%MOD)%MOD+MOD)%MOD;
    }
    return Res;
}

LL Calc(LL n){
    LL Res=0;
    LL m=sqrt(n);
    for(LL L=1,R;L<=m;L=R+1){
        R=min((LL)sqrt(n/(n/L/L)),m);
        Res=((Res+(Sum[R]-Sum[L-1])%MOD*Sigma1(n/L/L)%MOD)%MOD+MOD)%MOD;
    }
    return Res;
}

int main(){
    Get_Mu(maxn);
    int T;LL n;
    Read(T);
    while(T--){
        Read(n);
        printf("%lld
",Calc(n));
    }
    return 0;
}