【[SDOI2014]数表】
求
\[\sum_{i=1}^N\sum_{j=1}^Mσ(gcd(i,j))[σ(gcd(i,j))<=a]
\]
\(σ\)表示约数和函数
感觉非常难求的样子
先把套路搞出来
\[f(n)=\sum_{i=1}^N\sum_{j=1}^M[(i,j)=n]
\]
\[F(n)=\sum_{i=1}^N\sum_{j=1}^M[n|(i,j)]=\sum_{n|d}f(d)=\left \lfloor \frac{N}{n} \right \rfloor\times\left \lfloor \frac{M}{n}\right \rfloor
\]
我们要求的是
\[Ans=\sum_{i=1}^Nf(i)σ(i)=\sum_{i=1}^Nσ(i)\sum_{i|d}\mu(\frac{d}{i})\left \lfloor \frac{N}{d} \right \rfloor\times\left \lfloor \frac{M}{d}\right \rfloor
\]
可以把\(d\)放到前面来枚举
\[Ans=\sum_{d=1}^N \left \lfloor \frac{N}{d} \right \rfloor\times\left \lfloor \frac{M}{d}\right \rfloor\sum_{i|d}\mu(\frac{d}{i})σ(i)
\]
我们发现后面就是\(σ\times \mu\)
非常显然
\[σ=I\times id
\]
两边卷上\(\mu\)
\[σ\times \mu=I\times \mu\times id=id
\]
所以如果只是求上面那个柿子甚至线筛都不用
但是有一个限制
\[σ(i)<=a
\]
我们得好好考虑一下了
首先由于是多组询问，我们可以直接离线掉，按照每组询问的\(a\)从小到大排序
之后把所有的数按照\(σ\)来排序
每次把所有的\(σ<=a\)的数加入进去，显然这个数加进去的时候会影响其所有倍数，同时我们整除分块的时候需要前缀和，于是这里用一个树状数组维护
之后就很好做了
代码
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
#define re register
#define maxn 100005
#define uint long long
#define LL long long
#define max(a,b) ((a)>(b)?(a):(b))
#define min(a,b) ((a)<(b)?(a):(b))
#define lowbit(x) ((x)&(-x))
const LL mod=1ll<<31;
inline int read()
{
char c=getchar();
int x=0;
while(c<'0'||c>'9') c=getchar();
while(c>='0'&&c<='9') x=(x<<3)+(x<<1)+c-48,c=getchar();
return x;
}
struct ASK
{
int n,m,rk,A;
uint ans;
}a[20005];
inline int cmp(ASK A,ASK B)
{
return A.A<B.A;
}
int mu[maxn],f[maxn],p[maxn>>1],pre[maxn],d[maxn];
int T,N;
struct Node
{
int num,p;
}b[maxn];
uint c[maxn];
uint Ans[maxn];
inline void add(int x,uint val){for(re int i=x;i<=N;i+=lowbit(i)) c[i]=(c[i]+val)%mod;}
inline uint ask(int x){uint cnt=0;for(re int i=x;i;i-=lowbit(i)) cnt=(cnt+c[i])%mod;return cnt;}
inline int cop(Node A,Node B)
{
return A.num<B.num;
}
void write(uint x)
{
if(x>9) write(x/10);
putchar(x%10+48);
}
int main()
{
T=read();
for(re int i=1;i<=T;i++) a[i].n=read(),a[i].m=read(),a[i].rk=i,scanf("%d",&a[i].A);
for(re int i=1;i<=T;i++) if(a[i].n>a[i].m) std::swap(a[i].n,a[i].m);
for(re int i=1;i<=T;i++) N=max(N,a[i].m);
std::sort(a+1,a+T+1,cmp);
f[1]=1,mu[1]=1;
for(re int i=2;i<=N;i++)
{
if(!f[i]) p[++p[0]]=i,mu[i]=-1;
for(re int j=1;j<=p[0]&&p[j]*i<=N;j++)
{
f[p[j]*i]=1;
if(i%p[j]==0) break;
mu[p[j]*i]=-1*mu[i];
}
}
for(re int i=1;i<=N;i++)
for(re int j=1;j*i<=N;j++) d[i*j]+=i;
for(re int i=1;i<=N;i++) b[i].p=i,b[i].num=d[i];
std::sort(b+1,b+N+1,cop);
int now=1;
for(re int i=1;i<=T;i++)
{
while(now<=N&&b[now].num<=a[i].A)
{
for(re int j=1;j*b[now].p<=N;j++) add(b[now].p*j,((b[now].num*mu[j])%mod+mod)%mod);
now++;
}
for(re uint l=1,r;l<=a[i].n;l=r+1)
{
r=min(a[i].n/(a[i].n/l),a[i].m/(a[i].m/l));
a[i].ans+=(uint)(a[i].n/l)%mod*(uint)(a[i].m/l)%mod*(ask(r)-ask(l-1)+mod)%mod;
a[i].ans%=mod;
}
Ans[a[i].rk]=a[i].ans;
}
for(re int i=1;i<=T;i++) write(Ans[i]),putchar(10);
return 0;
}
