[SDOI2015]约数个数和

Description
设d(x)为x的约数个数,给定N、M,求\(\sum\limits_{i=1}^n\sum\limits_{j=1}^md(ij)\)

Input
输入文件包含多组测试数据。
第一行,一个整数T,表示测试数据的组数。
接下来的T行,每行两个整数N、M。

Output
T行,每行一个整数,表示你所求的答案。

Sample Input
2
7 4
5 6

Sample Output
110
121

HINT
\(1\leqslant N, M\leqslant50000\)
\(1\leqslant T\leqslant50000\)


\(d\)太丑了,我之后都写成\(\sigma_0\)

首先你需要知道\(\sigma_0(nm)=\sum\limits_{i|n}\sum\limits_{j|m}[\gcd(i,j)=1]\)

(我就是不知道这个推了1h无果)

然后我们证明一下这个东西(我只能证明等式左右两边相等)

\(n=\prod\limits_{i=1}^kp_i^{a_i},m=\prod\limits_{i=1}^kp_i^{b_i}\)\(a_i,b_i\)可以为0

然后有\(\sigma_0(nm)=\prod\limits_{i=1}^k(a_i+b_i+1)\),我们考虑某个质因子\(p_i\)对答案的贡献, \(p_i\)要么在\(n\)中有\(a_i\)中取法,要么在\(m\)中有\(b_i\)种取法,要么不取有一种方法,所以共计\(a_i+b_i+1\)中取法,再乘法原理累乘一下即为\(\sigma_0(n,m)\)

然后我们来推式子

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m\sum\limits_{u|i}\sum\limits_{v|j}[\gcd(u,v)=1]\]

我们调整一下枚举顺序

\[\sum\limits_{u=1}^n\sum\limits_{v=1}^m\lfloor\dfrac{n}{u}\rfloor\lfloor\dfrac{m}{v}\rfloor[\gcd(u,v)=1]\]

\(u,v\)太丑啦,我们改一下变量名,并且对最后用面的式子反演一下

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^m\lfloor\dfrac{n}{i}\rfloor\lfloor\dfrac{m}{j}\rfloor\sum\limits_{x|i,x|j}\mu(x)\]

我们把枚举\(x\)提到前面,并且进行分组,得到

\[\sum\limits_{x=1}^n\mu(x)\sum\limits_{i=1}^{\lfloor\frac{n}{x}\rfloor}\lfloor\dfrac{n}{ix}\rfloor\sum\limits_{j=1}^{\lfloor\frac{m}{x}\rfloor}\lfloor\dfrac{m}{jx}\rfloor\]

如果我们能够求出\(\sum\limits_{i=1}^n\lfloor\dfrac{n}{i}\rfloor\),我们就可以对原式进行数论分块了

然后我们参考一下[AHOI2005]约数研究这题,可以发现\(\sum\limits_{i=1}^n\lfloor\dfrac{n}{i}\rfloor=\sum\limits_{i=1}^n\sigma_0(i)\),因为\(\lfloor\dfrac{n}{i}\rfloor\)相当于统计\(1\sim n\)中有多少个数的因数有\(i\)

\(\sigma_0(n)\)显然是个积性函数,因此我们可以线筛求出前缀和,然后对原式进行数论分块即可

/*program from Wolfycz*/
#include<cmath>
#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
#define inf 0x7f7f7f7f
using namespace std;
typedef long long ll;
typedef unsigned int ui;
typedef unsigned long long ull;
inline char gc(){
    static char buf[1000000],*p1=buf,*p2=buf;
    return p1==p2&&(p2=(p1=buf)+fread(buf,1,1000000,stdin),p1==p2)?EOF:*p1++;
}
inline int frd(){
    int x=0,f=1; char ch=gc();
    for (;ch<'0'||ch>'9';ch=gc())   if (ch=='-')    f=-1;
    for (;ch>='0'&&ch<='9';ch=gc()) x=(x<<3)+(x<<1)+ch-'0';
    return x*f;
}
inline int read(){
    int x=0,f=1; char ch=getchar();
    for (;ch<'0'||ch>'9';ch=getchar())  if (ch=='-')    f=-1;
    for (;ch>='0'&&ch<='9';ch=getchar())    x=(x<<3)+(x<<1)+ch-'0';
    return x*f;
}
inline void print(int x){
    if (x<0)    putchar('-'),x=-x;
    if (x>9)    print(x/10);
    putchar(x%10+'0');
}
const int N=5e4;
int prime[N+10],mu[N+10],smu[N+10],f[N+10],Frm[N+10];
ll sf[N+10];
bool inprime[N+10];
void prepare(){
    int tot=0;
    mu[1]=f[1]=smu[1]=sf[1]=1;
    for (int i=2;i<=N;i++){
        if (!inprime[i])    prime[++tot]=i,mu[i]=-1,f[i]=Frm[i]=2;
        for (int j=1;j<=tot&&i*prime[j]<=N;j++){
            inprime[i*prime[j]]=1;
            if (i%prime[j]==0){
                mu[i*prime[j]]=0;
                Frm[i*prime[j]]=Frm[i]+1;
                f[i*prime[j]]=f[i]/Frm[i]*Frm[i*prime[j]];
                break;
            }
            mu[i*prime[j]]=-mu[i];
            Frm[i*prime[j]]=2;
            f[i*prime[j]]=f[i]<<1;
        }
        sf[i]=sf[i-1]+f[i];
        smu[i]=smu[i-1]+mu[i];
    }
}
int main(){
    prepare();
    for (int T=read();T;T--){
        int n=read(),m=read(); ll Ans=0;
        if (n>m)    swap(n,m);
        for (int i=1,pos;i<=n;i=pos+1){
            pos=min(n/(n/i),m/(m/i));
            Ans+=(smu[pos]-smu[i-1])*sf[n/i]*sf[m/i];
        }
        printf("%lld\n",Ans);
    }
    return 0;
}

转载于:https://www.cnblogs.com/Wolfycz/p/10270454.html

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值