[bzoj2693] jzptab

Description

Input

一個正整數T表示數據組數

接下來T行 每行兩個正整數 表示N、M

Output

T行 每行一個整數 表示第i組數據的結果

Sample Input

1
4 5

Sample Output

122

HINT

T <= 10000
N, M<=10000000

Solution

前置知識：莫比烏斯反演

推一下題目中的式子，莫比烏斯反演下，可得：
\[ \begin{align} ans&=\sum_{i=1}^n\sum_{j=1}^m\frac{i*j}{gcd(i,j)}\\ &=\sum_{d=1}^{min(n,m)}d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}i*j[gcd(i,j)=1]\\ &=\sum_{d=1}^{min(n,m)}d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}i*j \sum_{d^\prime|i\&d^\prime|j}\mu(d^\prime) \\ &=\sum_{d=1}^{min(n,m)}d\sum_{d^\prime}\mu(d^\prime)f(\lfloor\frac{n}{dd^\prime}\rfloor,\lfloor\frac{m}{dd^\prime}\rfloor){d^\prime}^2 \\ &=\sum_{T=1}^{min(n,m)}f(\lfloor\frac{n}{T}\rfloor,\lfloor\frac{m}{T}\rfloor)\sum _{d|T}d\mu(\frac{T}{d})(\frac{T}{d})^2\\ &=\sum_{T=1}^{min(n,m)}f(\lfloor\frac{n}{T}\rfloor,\lfloor\frac{m}{T}\rfloor)T\sum _{d|T}\mu(d)*d\\ \end {align} \]
其中，定義\(f\)爲：
\[ \begin{align} f(n,m)&=\sum_{i=1}^n\sum_{j=1}^mi*j\\ &=\sum_{i=1}^n\frac{i*m*(m+1)}{2}\\ &=\frac{n*(n+1)*m*(m+1)}{4} \end{align} \]
\(f\)那部分數論分塊搞搞就好了。

而後考慮\(ans\)的最後一部分，設爲\(g\)，即：
\[ g(n)=n*\sum_{d|n}\mu(d)*d \]
因爲\(n\)較大，考慮線篩這個東西，設質數\(p\)，當\(p\nmid n\)時，展開\(g(n*p)\)：
\[ \begin{align} g(n*p)&=n*p*(\sum_{d|n}\mu(d)*d+\sum_{d|n}\mu(d*p)*d*p)\\ &=p*(g(n)+g(n)*(-p))\\ &=p*(1-p)*g(n)\\ &=g(p)*g(n) \end{align} \]
不然同理，可得\(g(n*p)=g(n)*p\)。

而後，很重要的一點是，1e8+9不是質數！！中間要用到的\(inv4\)能夠線篩或\(exgcd\)等等

線篩逆元：
\[ inv[i]=(mod-mod/i)*inv[mod\%i]\%mod; \]
時間複雜度\(O(n+q\sqrt{n})\)。

#include<bits/stdc++.h>
using namespace std;

#define int long long 

void read(int &x) {
    x=0;int f=1;char ch=getchar();
    for(;!isdigit(ch);ch=getchar()) if(ch=='-') f=-f;
    for(;isdigit(ch);ch=getchar()) x=x*10+ch-'0';x*=f;
}

void print(int x) {
    if(x<0) putchar('-'),x=-x;
    if(!x) return ;print(x/10),putchar(x%10+48);
}
void write(int x) {if(!x) putchar('0');else print(x);putchar('\n');}

const int maxn = 1e7+1;
const int mod = 100000009;

int pri[maxn],vis[maxn],f[maxn],tot,n,m,k,p[maxn],inv4,inv[maxn];

int qpow(int a,int x) {
    int res=1;
    for(;x;x>>=1,a=a*a%mod) if(x&1) res=res*a%mod;
    return res;
}

void sieve() {
    f[1]=1;
    for(int i=2;i<maxn;i++) {
        if(!vis[i]) pri[++tot]=i,f[i]=i*(1-i)%mod;
        for(int j=1;j<=tot&&i*pri[j]<maxn;j++) {
            vis[i*pri[j]]=1;
            if(!(i%pri[j])) {f[i*pri[j]]=f[i]*pri[j]%mod;break;}
            else f[i*pri[j]]=f[i]*f[pri[j]]%mod;
        }
    }inv[1]=1;
    for(int i=1;i<maxn;i++) {
        if(i!=1) inv[i]=(mod-mod/i)*inv[mod%i]%mod;
        f[i]=(f[i-1]+f[i])%mod;
    }
}

int calc(int n,int m) {return n*(n+1)%mod*m%mod*(m+1)%mod*inv4%mod;}

signed main() {
    sieve();int t;read(t);inv4=inv[4];
    while(t--) {
        int n,m;read(n),read(m);
        int T=1,ans=0;
        while(T<=n&&T<=m) {
            int pre=T;T=min(n/(n/T),m/(m/T));
            ans=(ans+calc(n/T,m/T)*(f[T]-f[pre-1])%mod)%mod;
            T++;
        }
        write((ans%mod+mod)%mod);
    }
    return 0;
}