设d(x)为x的约数个数,给定N、M,求∑i=1N∑j=1Md(ij)
链接🔗
「SDOI2015」约数和个数
题解
约数个数公式 (不会证鸭)
d(ij)=x∣i∑y∣j∑[gcd(x,y)=1]
原式等于
i=x∑nj=y∑mi∣x∑j∣y∑[gcd(i,j)==1]
将枚举因数的位置调换一下
i=1∑nj=1∑m⌊in⌋⌊jm⌋[gcd(i,j)=1]
然后我们设
f(x)=i=1∑nj=1∑m⌊in⌋⌊jm⌋[gcd(i,j)=x]
再设
g(x)=x∣d∑f(d)
等式展开
g(x)=i=1∑nj=1∑m⌊in⌋⌊jm⌋[x∣gcd(i,j)]
再将x提出来
g(x)=i=1∑⌊xn⌋j=1∑⌊xm⌋⌊ixn⌋⌊jxm⌋
拆开来
g(x)=i=1∑⌊xn⌋⌊ixn⌋j=1∑⌊xm⌋⌊jxm⌋
我们设 h(x)=∑i=1n⌊ix⌋
那么g(x)=h(⌊xn⌋)h(⌊xm⌋)
我们要求的东西是 f(1) ,然后通过莫比乌斯反演,可以得出我们要求得东西是这个玩意
f(x)=x∣d∑μ(xdg(d))
f(1)=1∣d∑nμ(1d)h(⌊dn⌋)h(⌊dm⌋)
μ 是可以预先筛出来的,然后我们可以预处理楚h(i) 的值,然后就可以数论分块,将每次查询降为 O(n)
代码
#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
const int maxd = 5e4+10;
int m,n,ans,k,t;
int mu[maxd], prime[maxd],summu[maxd],cnt;
long long f[maxd];
bool vis[maxd];
void init()
{
mu[1] = 1;
for(int i=2;i<maxd;i++)
{
if(!vis[i]) prime[++cnt] = i, mu[i] = -1;
for(int j=1;j<=cnt && i*prime[j] <= maxd;j++)
{
vis[i*prime[j]] = 1;
if(i%prime[j] == 0) break;
mu[i*prime[j]] = -mu[i];
}
}
for(int i=1;i<maxd;i++) summu[i] = summu[i-1] + mu[i];
for(int i=1;i<=5e4;i++) //预处理h(i)
for(int l=1,r;l<=i;l = r+1) r = i/(i/l),f[i]+=1LL * (r-l+1) * (i/l);
}
long long sum(int m,int n) //数论分块
{
long long t = min(n, m), ans = 0;
for (int i = 1, j = 1; i <= t; i = j + 1)
{
j = min(n / ( n / i), m / (m / i));
ans += (summu[j] - summu[i - 1]) * f[n / i] * f[m / j];
}
return ans;
}
int main()
{
// freopen("a.in","r",stdin);
// freopen("k.out","w",stdout);
init();scanf("%d",&t);
for(int i=1;i<=t;i++)
{
scanf("%d %d",&n,&m);
printf("%lld\n",sum(n,m));
}
return 0;
}