[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

前置知识：莫比乌斯反演

推一下题目中的式子，莫比乌斯反演下，可得：

[egin{align} ans&=sum_{i=1}^nsum_{j=1}^mfrac{i*j}{gcd(i,j)}\ &=sum_{d=1}^{min(n,m)}dsum_{i=1}^{lfloorfrac{n}{d} floor}sum_{j=1}^{lfloorfrac{m}{d} floor}i*j[gcd(i,j)=1]\ &=sum_{d=1}^{min(n,m)}dsum_{i=1}^{lfloorfrac{n}{d} floor}sum_{j=1}^{lfloorfrac{m}{d} floor}i*j sum_{d^prime|i&d^prime|j}mu(d^prime) \ &=sum_{d=1}^{min(n,m)}dsum_{d^prime}mu(d^prime)f(lfloorfrac{n}{dd^prime} floor,lfloorfrac{m}{dd^prime} floor){d^prime}^2 \ &=sum_{T=1}^{min(n,m)}f(lfloorfrac{n}{T} floor,lfloorfrac{m}{T} floor)sum _{d|T}dmu(frac{T}{d})(frac{T}{d})^2\ &=sum_{T=1}^{min(n,m)}f(lfloorfrac{n}{T} floor,lfloorfrac{m}{T} floor)Tsum _{d|T}mu(d)*d\ end {align} ]

其中，定义(f)为：

[egin{align} f(n,m)&=sum_{i=1}^nsum_{j=1}^mi*j\ &=sum_{i=1}^nfrac{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 mid n)时，展开(g(n*p))：

[egin{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+qsqrt{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('
');}

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;
}