BZOJ 2693 jzptab

Crash的数字表格升级版。
其实就是如果答案式子有多重嵌套的话，
可以通过把最深的提到最上面，
其他的系数尝试线性筛（杜教筛）来完成。

#include<bits/stdc++.h>
#define mod 100000009
#define maxn 10000005
using namespace std;

/*

sigma(i,j,i*j/gcd(i,j))
=sigma(d,i,j,[gcd(i,j) == 1] * d * i * j)
=sigma(d,d * sigma(i,j,[gcd(i,j)==1] * i * j))

f(k) = sigma(i,j,[gcd(i,j)==k] * i * j)
g(k) = sigma(i,j,[k|gcd(i,j)] * i * j) = sigma(f(d) , k|d)
= k * k * (n/k+1)(n/k)/2 * (m/k+1)(m/k)/2
f(1) = sigma(d,mu[d]*g(d))
ans = sigma(d * sigma(p , mu[p] * p * p * (n/D+1)(n/D)/2 * (m/D+1)(m/D)/2))
= sigma( (n/D+1)(n/D)/2 * (m/D+1)(m/D)/2 * sigma(p|D , mu[p] * p * D))

*/

int mu[maxn],sum[maxn],pr[maxn],cnt_pr;
bool vis[maxn];

void sieve()
{
	mu[1] = 1 , sum[1]=1;
	for(int i=2;i<maxn;i++)
	{
		if(!vis[i]) pr[cnt_pr++]=i,mu[i]=-1,sum[i]=1-i;
		for(int j=0;j<cnt_pr&&pr[j]*i<maxn;j++)
		{
			vis[i * pr[j]] = 1;
			if(i % pr[j] == 0)
			{
				sum[i * pr[j]] = sum[i];
 				break;
			}
			mu[i * pr[j]] = -mu[i];
			sum[i * pr[j]] = (1ll * sum[i] * sum[pr[j]]) % mod;
		}
	}
	for(int i=2;i<maxn;i++) sum[i] = (sum[i-1] + 1ll * sum[i] * i) % mod;
}

inline int sm(int n)
{
	return 1ll * n * (n+1) / 2 % mod;
}

int main()
{
	sieve();
	int T,n,m;
	for(scanf("%d",&T);T--;)
	{
		scanf("%d%d",&n,&m);
		if(n > m) swap(n,m);
		int ans = 0;
		for(int i=1,nxt;i<=n;i=nxt+1)
		{
			nxt = min(n/(n/i),m/(m/i));
			ans = (ans + 1ll * (sum[nxt] - sum[i-1]) * sm(n/i) % mod * sm(m/i) %mod)%mod;
		}
		printf("%d\n",(ans+mod)%mod);
	}
}