约数个数和「SDOI2015」

题意

已知(n,m)，求(sum^n_asum^m_b d(ab))，其中(d(x))表示(x)的约数个数。

---

思路

结论：(d(nm)=sum_{i|n}sum_{j|m}[gcd(i,j)==1])。（证明见后）

由此可得原式等价于(sum^n_asum^m_bsum_{i|n}sum_{j|m}[gcd(i,j)==1])

即(sum^n_isum^m_jfrac{n}{i}frac{m}{j}[gcd(i,j)==1])

设(f(d)=sum^n_isum^m_jfrac{n}{i}frac{m}{j}[gcd(i,j)==d])

(F(d)=sum^n_isum^m_jfrac{n}{i}frac{m}{j}[gcd(i,j)|d]=sum^{frac{n}{d}}_isum^{frac{m}{d}}_j frac{n}{id}frac{m}{jd})

那么有(F(d)=sum_{d|i} f(i))

所以(f(i)=sum_{i|d}F(d)mu{frac{d}{i}})

代码

#include <bits/stdc++.h>

using namespace std;

namespace StandardIO {

	template<typename T> inline void read (T &x) {
		x=0;T f=1;char c=getchar();
		for (; c<'0'||c>'9'; c=getchar()) if (c=='-') f=-1;
		for (; c>='0'&&c<='9'; c=getchar()) x=x*10+c-'0';
		x*=f;
	}
	template<typename T> inline void write (T x) {
		if (x<0) putchar('-'),x=-x;
		if (x>=10) write(x/10);
		putchar(x%10+'0');
	}

}

using namespace StandardIO;

namespace Solve {
	
	const int N=50001;
	
	int n,top;
	int a,b;
	int isprime[N],prime[N];
	long long mu[N],pre[N];
	
	inline void init () {
		mu[1]=1;
		for (register int i=2; i<=N; ++i) {
			if (!isprime[i]) {
				prime[++top]=i,mu[i]=-1;
			}
			for (register int j=1; j<=top&&i*prime[j]<=N; ++j) {
				isprime[i*prime[j]]=1;
				if (i%prime[j]==0) {
					mu[i*prime[j]]=0;
					break;
				}
				mu[i*prime[j]]=-mu[i];
			}
		}
		for (register int i=1; i<=N; ++i) mu[i]+=mu[i-1];
		for (register int x=1; x<=N; ++x) {
			for (register int i=1,r; i<=x; i=r+1) {
				r=x/(x/i);
				pre[x]+=(long long)(x/i)*(r-i+1);
			}
		}
	}
	inline long long calc (int m,int n) {
		long long sum=0;
		for (register int i=1,r; i<=min(m,n); i=r+1) {
			r=min(n/(n/i),m/(m/i));
			sum+=pre[n/i]*pre[m/i]*(mu[r]-mu[i-1]);
		}
		return sum;
	}

	inline void MAIN () {
		init();
		read(n);
		for (register int i=1; i<=n; ++i) {
			read(a),read(b);
			write(calc(a,b)),putchar('
');
		}
	}
	
}

int main () {
	Solve::MAIN();
}

---

证明

orz Siyuan