【莫比乌斯反演】[HYSBZ/BZOJ2693]jzptab

题目大意就是求∑ni=1∑mj=1lcm(i,j)，但是有多组输入数据.
如果之前做过【莫比乌斯反演】[HYSBZ\BZOJ2154]Crash的数字表格,就会发现，对于每一个询问，有O(n)的做法，但显然不够快。
在上一道题中(sum的定义见【莫比乌斯反演】[HYSBZ\BZOJ2154]Crash的数字表格)。

ans=∑d=1min(n,m)d∑k=1min(m,n)μ(k)∗k2∗sum(⌊nkd⌋,⌊mkd⌋)

令D=kd,稍加变形可得

ans=∑D=1min(n,m)sum(⌊nD⌋,⌊mD⌋)∑i|dμ(i)∗i2∗Di

当我们对r∑i|dμ(i)∗i2∗Di在做筛法的同时求一个前缀和，很显然就可以用分块优化以O(n√)的时间回答每次询问。

#include<cstdio>
#include<algorithm>
using namespace std;
#define MAXN 10000000
#define MOD 100000009
#define Sum(x,y) ((1ll*x*(x+1)/2%MOD)*(1ll*y*(y+1)/2%MOD)%MOD)
bool f[MAXN+10];
int p[MAXN+10],mu[MAXN+10],sum[MAXN+10],pcnt,T,n,m,ans;
void prepare(){
    int i,j;
    mu[1]=sum[1]=1;
    for(i=2;i<=MAXN;i++){
        if(!f[i])
            p[++pcnt]=i,mu[i]=-1,sum[i]=(i*1*1+1ll*i*i*-1)%MOD;
        for(j=1;i*p[j]<=MAXN;j++){
            f[i*p[j]]=1;
            if(i%p[j]==0){
                mu[i*p[j]]=0;
                sum[i*p[j]]=(1ll*sum[i]*p[j])%MOD;
                break;
            }
            sum[i*p[j]]=((1ll*sum[i]*p[j]-1ll*sum[i]*p[j]%MOD*p[j])%MOD+MOD)%MOD;
        }
        sum[i]=(sum[i-1]+sum[i])%MOD;
    }
}
void Read(int &x){
    char c;
    while(c=getchar(),c!=EOF)
        if(c>='0'&&c<='9'){
            x=c-'0';
            while(c=getchar(),c>='0'&&c<='9')
                x=x*10+c-'0';
            ungetc(c,stdin);
            return;
        }
}
void solve(){
    int i,last,t=min(n,m);
    ans=0;
    for(i=1;i<=t;i=last+1){
        last=min(n/(n/i),m /(m /i));
        ans=((ans+(sum[last]-sum[i-1])*Sum(n/i,m /i)%MOD)%MOD+MOD)%MOD;
    }
}
int main()
{
    Read(T);
    prepare();
    while(T--){
        Read(n),Read(m);
        solve();
        printf("%d\n",ans);
    }
}