# [BZOJ2693]:jzptab

## Description

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

## Output

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

1
4 5

122

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}d\sum_{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^2\sum_{i=1}^{x/d}\sum_{j=1}^{y/d}{ij}$

$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^2\sum_{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_{T=1}^{n}{B(n/T)B(m/T)\sum_{i|T}{\mu(i)i^2\frac{T}i}}$

## 代码

#include<cstdio>
#include<algorithm>
# define LL long long
const int M = 10000005 ;
const int mod = 100000009 ;
using namespace std ;
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\n",(Ans % mod + mod) % mod) ;
}
return 0 ;
}

posted @ 2019-02-13 22:00  beretty  阅读(98)  评论(0编辑  收藏  举报