bzoj 2693 jzptab 莫比乌斯反演

#include<cmath>
#include<cstdio>
#include<iostream>
#define LL long long
#define maxn 10000005
using namespace std;
LL mod=100000009;
LL sum(LL x,LL y)
{
x=(x*(x+1)/2)%mod;
y=(y*(y+1)/2)%mod;
return x*y%mod;
}
LL h[maxn];
int mu[maxn];
bool use[maxn];
int prime[maxn],cnt;
void get_mu()
{
mu[1]=1;h[1]=1;
int N=10000000;
for(int i=2;i<=N;i++)
{
if(!use[i])
{
prime[++cnt]=i;mu[i]=-1;
h[i]=i-(LL)i*i;
h[i]%=mod;
}
for(int j=1;j<=cnt;j++)
{
if((LL)i*prime[j]>(LL)N) break;
use[i*prime[j]]=1;
if(i%prime[j]==0)
{
h[i*prime[j]]=(h[i]*prime[j])%mod;
mu[i*prime[j]]=0;
break;
}
h[i*prime[j]]=(h[i]*h[prime[j]])%mod;
mu[i*prime[j]]=-mu[i];
}
}
for(int i=1;i<=N;i++)
{
h[i]+=h[i-1];
h[i]%=mod;
}
}
int main()
{
LL n,m;
get_mu();
int T;scanf("%d",&T);
while(T--)
{
scanf("%lld%lld",&n,&m);
if(n>m) swap(n,m);
LL l=1,r,ans=0;
while(l<=n)
{
int p1=n/l,p2=m/l;
r=min(n/p1,m/p2);
ans+=(sum(n/l,m/l)*(h[r]-h[l-1]+mod))%mod;
ans%=mod;l=r+1;
}
if(ans<0) ans+=mod;
ans%=mod;
printf("%lld\n",ans);
}
return 0;
}