SPOJ 20173 DIVCNT2 - Counting Divisors (square)（数论+莫比乌斯反演）

Description

给出一整数$n$，求$S(n)=\sum\limits_{i=1}^n\sigma_0(i^2)$，其中$\sigma_0(x)$为$x$的因子数

Input

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

Output

输出$S(n)$

Sample Input

5
1
2
3
10
100

Sample Output

1
4
7
48
1194

Solution

对$x$质因子分解有$x=p_1^{a_1}...p_k^{a_k}$，进而

$$\sigma_0(x^2)=\prod\limits_{i=1}^k(2a_i+1)=\sum\limits_{S\subset \{1,...,k\}}2^{|S|}\prod\limits_{i\in S}a_i=\sum\limits_{d|x}2^{\omega(d)}$$
其中$\omega(d)$表示$d$的素因子个数

（累乘$\prod\limits_{i\in S}a_i$即为素因子集合为$S$且每个素因子的指数不超过$a_i$的数字个数，也即$x$的满足素因子集合为$S$的因子个数，而这些数字每个对答案贡献均为$2^{|S|}$，也即为这些数字的素因子个数）

注意到$2^{\omega(d)}$即为$d$的为无平方因子数的因子个数，故有

$$2^{\omega(d)}=\sum\limits_{j|d}\mu^2(j)$$
进而

$$S(n)=\sum\limits_{i=1}^n\sum\limits_{d|i}\sum\limits_{j|d}\mu^2(j)=\sum\limits_{d=1}^n\lfloor\frac{n}{d}\rfloor\sum\limits_{j|d}\mu^2(j)=\sum\limits_{j=1}\mu^2(j)\sum\limits_{i=1}^{\lfloor\frac{n}{j}\rfloor}\lfloor\frac{n}{ij}\rfloor$$
由于$\lfloor\frac{n}{j}\rfloor$的取值只有$O(\sqrt{n})$种，对每种取值可以快速计算$\mu^2(j)$的区间和以及$\sum\limits_{i=1}^n\lfloor\frac{n}{i}\rfloor$即可

设$F(n)=\sum\limits_{i=1}^n\mu^2(i)$，则有

$$F(n)=\sum\limits_{i=1}^n\mu^2(i)=\sum\limits_{i=1}^n\sum\limits_{d^2|i}\mu(d)=\sum\limits_{d=1}^{\sqrt{n}}\mu(d)\lfloor\frac{n}{d^2}\rfloor$$
故可以$O(\sqrt{n})$求出$F(n)$

设$G(n)=\sum\limits_{i=1}^n\lfloor\frac{n}{i}\rfloor=\sum\limits_{i=1}^n\sigma_0(i)$，则$G(n)$为积性函数，满足$G(1)=1,G(p)=2$，且分块后可以$O(\sqrt{n})$求出$G(n)$

预处理$n^{\frac{2}{3}}$范围内的$F,G$，对于较大的$n$可以$O(\sqrt{n})$求出$F(n),G(n)$

Code

#include<cstdio>
#include<algorithm>
#include<vector>
#include<cmath>
using namespace std;
typedef long long ll;
#define maxn 100000005
int m=100000000;
int mu[maxn],num[maxn],f[maxn],p[maxn/10],res=0;
ll g[maxn];
bool vis[maxn];
void get_prime(int n)
{
    mu[1]=1,f[1]=1,g[1]=1;
    for(int i=2;i<=n;i++)
    {
        if(!vis[i])p[res++]=i,mu[i]=-1,g[i]=2,num[i]=1;
        f[i]=f[i-1]+(mu[i]?1:0);
        for(int j=0;j<res&&i*p[j]<=n;j++)
        {
            vis[i*p[j]]=1;
            if(i%p[j])
            {
                mu[i*p[j]]=-mu[i];
                g[i*p[j]]=g[i]*2;
                num[i*p[j]]=1;
            }
            else
            {
                mu[i*p[j]]=0;
                num[i*p[j]]=num[i]+1;
                g[i*p[j]]=g[i]/(num[i]+1)*(num[i]+2);
                break;
            }
        }
    }
    for(int i=2;i<=n;i++)g[i]+=g[i-1];
}
ll F(ll n)
{
    if(n<=m)return f[n];
    int nn=(int)sqrt(n+0.5); 
    ll ans=0;
    for(int d=1;d<=nn;d++)
        if(mu[d])ans+=mu[d]*(n/d/d);
    return ans;
}
ll G(ll n)
{
    if(n<=m)return g[n];
    ll ans=0;
    for(ll i=1,pre;i<=n;i=pre+1)
    {
        pre=n/(n/i);
        ans+=(pre-i+1)*(n/i);
    }
    return ans;
}
ll a[10005],mx;
int main()
{
    int T;
    scanf("%d",&T);
    for(int i=1;i<=T;i++)scanf("%lld",&a[i]),mx=max(mx,a[i]);
    mx=min(mx,(ll)m);
    get_prime(mx);
    for(int j=1;j<=T;j++)
    {
        ll ans=0,last=0,now,n=a[j];
        for(ll i=1,pre;i<=n;i=pre+1)
        {
            pre=n/(n/i);
            now=F(pre);
            ans+=(now-last)*G(n/i);
            last=now;
        }
        printf("%lld\n",ans);
    }
    return 0;
}