[BZOJ2693]:jzptab

Description

(sum_{i=1}^{n}sum_{j=1}^{m}{lcm(i,j)})

Input

一个正整数T表示数据组数
接下来T行 每行两个正整数 表示N、M

Output

T行 每行一个整数 表示第i组数据的结果

Sample Input

1
4 5

Sample Output

122

HINT

T <= 10000
N, M<=10000000

---

题解

抄的
(sum_{i=1}^{n}sum_{j=1}^{m}{i*j/gcd(i,j)})
(=sum_{d=1}^{n}sum_{i=1}^{n}sum_{j=1}^{m}{i j / d[gcd(i,j)==d]})
(sum_{d=1}^{n}sum_{i=1}^{n/d}sum_{j=1}^{m/d}{ijd[gcd(i,j)==1]})
(sum_{d=1}^{n}dsum_{i=1}^{n/d}sum_{j=1}^{m/d}{ij[gcd(i,j)==1]})

设(x=n/d,y=m/d)
(f(d)=sum_{i=1}^{x}sum_{j=1}^{y}{ij[gcd(i,j)==d]})
(F(d)=sum_{i=1}^{x}sum_{j=1}^{y}{ij[d|gcd(i,j)]})
(F(d)=d^2sum_{i=1}^{x/d}sum_{j=1}^{y/d}{ij})

那么(Ans=sum_{d=1}^{n}df(1))

(f(1)=sum_{i=1}^{x}sum_{j=1}^{y}{ij[gcd(i,j)==1]})
(f(1)=sum_{t=1}^{x}mu(t)F(t))
(f(1)=sum_{t=1}^{x}mu(t)t^2sum_{i=1}^{x/t}sum_{j=1}^{y/t}{ij})
(f(1)=sum_{t=1}^{x}mu(t)t^2B(x/t)B(y/t))，(其中B(x)=frac{(1 + x) * x}2)
(Ans=sum_{d=1}^{n}{d}sum_{t=1}^{n/d}mu(t)t^2B(n/dt)B(m/dt))

然后单独看后头那些东西
(P(x,y)=sum_{t=1}^{min(x,y)}{mu(t)t^2B(x/t)B(y/t)})
可以除法分块

然后(Ans=sum_{d=1}^{min(n,m)}{d}P(n/d,m/d))
这玩意儿是不是也能除法分块
所以复杂度就是(O(n))

然后还能不能化简?

交换求和号

(Ans=sum_{T=1}^{n}{B(n/T)B(m/T)sum_{i|T}{mu(i)i^2frac{T}i}})
这样以后后面的东西是一个积性函数,可以通过线筛求出前缀和
即可(O(sqrt{n}))回答询问

代码

#include<cstdio>
#include<algorithm>
# define LL long long
const int M = 10000005 ;
const int mod = 100000009 ;
using namespace std ;
inline int read() {
	char c = getchar() ; int x = 0 , w = 1 ;
	while(c>'9'||c<'0') { if(c=='-') w = -1 ; c = getchar() ; }
	while(c>='0'&&c<='9') { x = x*10+c-'0' ; c = getchar() ; }
	return x*w ;
}
 
bool notp[M] ;
int p[M] , pnum , sum[M] ;
int n , m , Ans , f[M] ;
void Get_Sum(int n) {
    sum[1] = 1 ;
    for(int i = 2 ; i <= n ; i ++) {
        if(!notp[i]) {
            p[++pnum] = i ;
            sum[i] = ((i - 1LL * i * i % mod) % mod + mod) % mod ;
        }
        for(int j = 1 ; j <= pnum && 1LL * i * p[j] <= n ; j ++) {
            notp[i * p[j]] = true ;
            if(i % p[j] == 0) {
                sum[i * p[j]] = (1LL * sum[i] * p[j] % mod + mod) % mod ;
                break ;
            }
            sum[i * p[j]] = (1LL * sum[i] * sum[p[j]] % mod + mod) % mod ;
        }
    }
    for(int i = 1 ; i <= n ; i ++) {
        sum[i] = ((sum[i - 1] + sum[i]) % mod + mod) % mod ;
        f[i] = (f[i - 1] + i) % mod ;
    }
}
int main() {
    int T = read() ; Get_Sum(10000000) ;
    while(T --) {
        n = read() ; m = read() ; if(n > m) swap(n , m) ; Ans = 0 ;
        for(int l = 1 , r ; l <= n ; l = r + 1) {
            r = min(n / (n / l) , m / (m / l)) ;
            Ans = (Ans + 1LL * f[n / l] * f[m / l] % mod * (sum[r] - sum[l - 1] + mod) % mod + mod) % mod ;
        }
        printf("%d
",(Ans % mod + mod) % mod) ;
    }
    return 0 ;
}