我是靠谱客的博主 危机猎豹,最近开发中收集的这篇文章主要介绍SPOJ 20173 DIVCNT2 - Counting Divisors (square)(数论+莫比乌斯反演),觉得挺不错的,现在分享给大家,希望可以做个参考。

概述

Description

给出一整数 n n ,求S(n)=i=1nσ0(i2),其中 σ0(x) σ 0 ( x ) x x 的因子数

Input

第一行一整数T表示用例组数,每组用例输入一整数 n(1T104,1n1012) n ( 1 ≤ T ≤ 10 4 , 1 ≤ n ≤ 10 12 )

Output

输出 S(n) S ( n )

Sample Input

5
1
2
3
10
100

Sample Output

1
4
7
48
1194

Solution

x x 质因子分解有x=p1a1...pkak,进而

σ0(x2)=i=1k(2ai+1)=S{1,...,k}2|S|iSai=d|x2ω(d) σ 0 ( x 2 ) = ∏ i = 1 k ( 2 a i + 1 ) = ∑ S ⊂ { 1 , . . . , k } 2 | S | ∏ i ∈ S a i = ∑ d | x 2 ω ( d )

其中 ω(d) ω ( d ) 表示 d d 的素因子个数

(累乘iSai即为素因子集合为 S S 且每个素因子的指数不超过ai的数字个数,也即 x x 的满足素因子集合为S的因子个数,而这些数字每个对答案贡献均为 2|S| 2 | S | ,也即为这些数字的素因子个数)

注意到 2ω(d) 2 ω ( d ) 即为 d d 的为无平方因子数的因子个数,故有

2ω(d)=j|dμ2(j)

进而

S(n)=i=1nd|ij|dμ2(j)=d=1nndj|dμ2(j)=j=1μ2(j)i=1njnij S ( n ) = ∑ i = 1 n ∑ d | i ∑ j | d μ 2 ( j ) = ∑ d = 1 n ⌊ n d ⌋ ∑ j | d μ 2 ( j ) = ∑ j = 1 μ 2 ( j ) ∑ i = 1 ⌊ n j ⌋ ⌊ n i j ⌋

由于 nj ⌊ n j ⌋ 的取值只有 O(n) O ( n ) 种,对每种取值可以快速计算 μ2(j) μ 2 ( j ) 的区间和以及 i=1nni ∑ i = 1 n ⌊ n i ⌋ 即可

F(n)=i=1nμ2(i) F ( n ) = ∑ i = 1 n μ 2 ( i ) ,则有

F(n)=i=1nμ2(i)=i=1nd2|iμ(d)=d=1nμ(d)nd2 F ( n ) = ∑ i = 1 n μ 2 ( i ) = ∑ i = 1 n ∑ d 2 | i μ ( d ) = ∑ d = 1 n μ ( d ) ⌊ n d 2 ⌋

故可以 O(n) O ( n ) 求出 F(n) F ( n )

G(n)=i=1nni=i=1nσ0(i) G ( n ) = ∑ i = 1 n ⌊ n i ⌋ = ∑ i = 1 n σ 0 ( i ) ,则 G(n) G ( n ) 为积性函数,满足 G(1)=1,G(p)=2 G ( 1 ) = 1 , G ( p ) = 2 ,且分块后可以 O(n) O ( n ) 求出 G(n) G ( n )

预处理 n23 n 2 3 范围内的 F,G F , G ,对于较大的 n n 可以O(n)求出 F(n),G(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("%lldn",ans);
}
return 0;
}

最后

以上就是危机猎豹为你收集整理的SPOJ 20173 DIVCNT2 - Counting Divisors (square)(数论+莫比乌斯反演)的全部内容,希望文章能够帮你解决SPOJ 20173 DIVCNT2 - Counting Divisors (square)(数论+莫比乌斯反演)所遇到的程序开发问题。

如果觉得靠谱客网站的内容还不错,欢迎将靠谱客网站推荐给程序员好友。

本图文内容来源于网友提供,作为学习参考使用,或来自网络收集整理,版权属于原作者所有。
点赞(31)

评论列表共有 0 条评论

立即
投稿
返回
顶部