2154的加強版

可以發現 由莫比烏斯函式組成的積性函式 在 prime[j]|i 轉移時 多出來的u(p*p*a)=0部分不記的 只需考慮F(x) 中的x 乘上了p對之前的影響

#include<cstdio>
#include<cstdlib>
#include<algorithm>
using namespace std;
typedef long long ll;

inline char nc()
{
	static char buf[100000],*p1=buf,*p2=buf;
	if (p1==p2) { p2=(p1=buf)+fread(buf,1,100000,stdin); if (p1==p2) return EOF; }
	return *p1++;
}

inline void read(int &x)
{
	char c=nc(),b=1;
	for (;!(c>='0' && c<='9');c=nc()) if (c=='-') b=-1;
	for (x=0;c>='0' && c<='9';x=x*10+c-'0',c=nc()); x*=b;
}

const int maxn=10000000,P=100000009;
int prime[maxn],vst[maxn],tot;

ll sum[maxn];

inline void Init(){
	sum[1]=1;
	for (int i=2;i<=maxn;i++)
	{
		if (!vst[i]){
			prime[++tot]=i; sum[i]=(ll)i*(i-1)%P; sum[i]=(P-sum[i])%P;
		}
		for (int j=1;(ll)prime[j]*i<=maxn && j<=tot;j++)
		{
			vst[i*prime[j]]=1;
			if (i%prime[j]==0){
				sum[i*prime[j]]=sum[i]*prime[j]%P;
				break;
			}
			sum[i*prime[j]]=sum[i]*sum[prime[j]]%P;
		}
	}
	for (int i=1;i<=maxn;i++)
		(sum[i]+=sum[i-1])%=P;
}

inline ll Sum(ll x,ll y){
	return (x*(x+1)/2%P)*(y*(y+1)/2%P)%P;
}

inline void Solve(int n,int m){
	if (n>m) swap(n,m);
	ll ans=0; int last;
	for (int i=1;i<=n;i=last+1)
	{
		last=min(n/(n/i),m/(m/i));
		ans+=Sum(n/i,m/i)*(sum[last]-sum[i-1]);
		ans=(ans%P+P)%P;
	}
	printf("%lld\n",ans);
}

int main()
{
	int T,n,m;
	Init();
	read(T);
	while (T--){
		read(n); read(m);
		Solve(n,m);
	}
	return 0;
}