[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}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}\)

那么\(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^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_{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^2\frac{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\n",(Ans % mod + mod) % mod) ;
    }
    return 0 ;
}