bzoj 2693 jzptab(线性筛预处理+反演)

$$\sum_{i=1}^a\sum_{j=1}^blcm(i,j)$$
$$\sum_{i=1}^a\sum_{j=1}^b\frac{ij}{gcd(i,j)}$$
$$\sum_{d}d\sum_{i=1}^{\frac{a}{d}}\sum_{j=1}^{\frac{b}{d}}ij[gcd(i,j)=1]$$
$$\sum_{d}\sum_{i=1}^{\frac{a}{d}}\sum_{j=1}^{\frac{b}{d}}ij\sum_{d^{'}|i且d^{'}|j}u(d^{'})$$
$$\sum_{d}d\sum_{d^{'}}u(d^{'})d^{'2}\sum_{i=1}^{\frac{a}{dd^{'}}}\sum_{j=1}^{\frac{b}{dd^{'}}}ij$$
$$\frac{1}{4}\sum_{d}d\sum_{d^{'}}u(d^{'})d^{'2}\frac{a}{dd^{'}}(\frac{a}{dd^{'}}+1)\frac{b}{dd^{'}}(\frac{a}{dd^{'}}+1)$$
$$dd^{'}=T$$
$$\frac{1}{4}\sum_{T=1}T\sum_{d|T}u(d)d\frac{a}{dd^{'}}(\frac{a}{dd^{'}}+1)\frac{b}{dd^{'}}(\frac{a}{dd^{'}}+1)$$
$可以用老方法nlogn的枚举因子去筛，不过这题数据量1000W，可能会很慢（nlogn筛跑8s，线性筛1s以内）$
$线性筛的话就是筛F(T)=\sum_{d|T}Tu(d)d$
$首先F(T)肯定是积性函数，所以没有出现过的素因子F(i*p)=F(i)*F(p)$
$当i\%p==0时候，怎么解呢$
$因为u(d)如果里面d的某个素因子超过1次的话，u(d)=0，所以p乘上之后u(d*p)=0$
$所以多出来的几项都是0，只有T变成了T*p$
$所以F(i*p)=F(i)*p$

---

代码：

#include <map>
#include <set>
#include <stack>
#include <queue>
#include <cmath>
#include <string>
#include <vector>
#include <cstdio>
#include <cctype>
#include <cstring>
#include <sstream>
#include <cstdlib>
#include <iostream>
#include <algorithm>
#pragma comment(linker,"/STACK:102400000,102400000")

using namespace std;
#define   MAX           10000005
#define   MAXN          1000005
#define   maxnode       205
#define   sigma_size    2
#define   lson          l,m,rt<<1
#define   rson          m+1,r,rt<<1|1
#define   lrt           rt<<1
#define   rrt           rt<<1|1
#define   middle        int m=(r+l)>>1
#define   LL            long long
#define   ull           unsigned long long
#define   mem(x,v)      memset(x,v,sizeof(x))
#define   lowbit(x)     (x&-x)
#define   pii           pair<int,int>
#define   bits(a)       __builtin_popcount(a)
#define   mk            make_pair
#define   limit         10000

//const int    prime = 999983;
const int    INF   = 0x3f3f3f3f;
const LL     INFF  = 0x3f3f;
//const double pi    = acos(-1.0);
const double inf   = 1e18;
const double eps   = 1e-9;
const LL     mod   = 1e8+9;
const ull    mx    = 133333331;

/*****************************************************/
inline void RI(int &x) {
      char c;
      while((c=getchar())<'0' || c>'9');
      x=c-'0';
      while((c=getchar())>='0' && c<='9') x=(x<<3)+(x<<1)+c-'0';
 }
/*****************************************************/

bool prime[MAX];
int pr[MAX];
int tot;
LL T[MAX];
void init(){
    mem(prime,0);
    mem(T,0);tot=0;
    T[1]=1;
    for(int i=2;i<MAX;i++){
        if(!prime[i]){
            pr[tot++]=i;
            T[i]=(i-(LL)i*i)%mod;
        }
        for(int j=0;j<tot&&pr[j]*i<MAX;j++){
            prime[pr[j]*i]=1;
            if(i%pr[j]==0){
                T[i*pr[j]]=T[i]*pr[j]%mod;    
                break;
            }
            T[i*pr[j]]=T[i]*T[pr[j]]%mod;
        }
    }
    for(int i=1;i<MAX;i++){
        T[i]=(T[i]+T[i-1])%mod;
        if(T[i]<0) T[i]+=mod;
    }
}

LL qpow(LL a,LL n){
    LL ans=1;
    while(n){
        if(n&1) ans=ans*a%mod;
        a=a*a%mod;
        n>>=1;
    }
    return ans;
}
int main() {
    freopen("bzoj_2693.in","r",stdin);
    freopen("bzoj_2693.out","w",stdout);
    int t;
    init();
    cin>>t;
    while(t--){
        int n,m;
        scanf("%d%d",&n,&m);
        LL ans=0;
        for(int i=1;i<=min(n,m);){
            int x=n/(n/i);
            int y=m/(m/i);
            int z=min(x,y);
            LL aa=(LL)(n/i)*(n/i+1)/2%mod;
            LL bb=(LL)(m/i)*(m/i+1)/2%mod;
            ans+=((T[z]-T[i-1])*aa%mod)*bb%mod;
            if(ans>=mod) ans-=mod;
            if(ans<0) ans+=mod;
            i=z+1;
        }
        cout<<(ans%mod+mod)%mod<<endl;
    }
    return 0;
}