约数个数和「SDOI2015」

题意

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


思路

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

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

\(\sum^n_i\sum^m_j\frac{n}{i}\frac{m}{j}[gcd(i,j)==1]\)

\(f(d)=\sum^n_i\sum^m_j\frac{n}{i}\frac{m}{j}[gcd(i,j)==d]\)

\(F(d)=\sum^n_i\sum^m_j\frac{n}{i}\frac{m}{j}[gcd(i,j)|d]=\sum^{\frac{n}{d}}_i\sum^{\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('\n');
		}
	}
	
}

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

证明

orz Siyuan

搜狗截图20190814211034.png

posted @ 2019-08-14 20:38  Ilverene  阅读(278)  评论(0编辑  收藏  举报