BZOJ-2154 Crash的数字表格(莫比乌斯反演 + 分块)

---

题目链接

---

Time limit 20000 ms Memory limit 265216 kB

Problem Description

今天的数学课上，Crash小朋友学习了最小公倍数(Least Common Multiple)。对于两个正整数a和b，LCM(a, b)表示能同时被a和b整除的最小正整数。例如，LCM(6, 8) = 24。回到家后，Crash还在想着课上学的东西，为了研究最小公倍数，他画了一张NM的表格。每个格子里写了一个数字，其中第i行第j列的那个格子里写着数为LCM(i, j)。一个45的表格如下： 1 2 3 4 5 2 2 6 4 10 3 6 3 12 15 4 4 12 4 20 看着这个表格，Crash想到了很多可以思考的问题。不过他最想解决的问题却是一个十分简单的问题：这个表格中所有数的和是多少。

当N和M很大时，Crash就束手无策了，因此他找到了聪明的你用程序帮他解决这个问题。由于最终结果可能会很大，Crash只想知道表格里所有数的和mod 20101009的值。

Input

输入的第一行包含两个正整数，分别表示N和M。

Output

输出一个正整数，表示表格中所有数的和mod 20101009的值。

Sample Input

4 5

Sample Output

122

【数据规模和约定】 100%的数据满足N, M ≤ 10^7。

---

题意:

中文题目！！！

思路:

莫比乌斯反演。

以下公式默认n < m;

$ans = \sum_{i = 1}^n \sum_{j=1}^mLCM(i, j)\\ = \sum_{i = 1}^n \sum_{j=1}^m \frac {i\ast j}{GCD(i, j)}\\ =\sum_{d = 1}^n \sum_{i = 1}^{n/d} \sum_{j=1}^{m/d} \frac{d\ast i \ast d \ast j}{d}[gcd(d\ast j, d \ast i) = d]\\ =\sum_{d = 1}^n d \sum_{i = 1}^{n/d} \sum_{j=1}^{m/d} i \ast j [gcd(i,j) = 1]\\$

现在我们要求的是后面的两重求和式子，感觉和普通的求gcd = 1的反演基本一致。

$设f(n, m, k) =\sum_{i = 1}^{n} \sum_{j=1}^{m} i \ast j [gcd(i,j) = k]\\ 设F(n, m, k) =\sum_{i = 1}^{n} \sum_{j=1}^{m} i \ast j [k | gcd(i,j) ]\\ = \sum_{i = 1}^{n/d} \sum_{j=1}^{m/d} i \ast k \ast j \ast k [k | gcd(i,j) ]\\ = k\ast k \sum_{i = 1}^{n/k} \sum_{j=1}^{m/k} i \ast j [ k | gcd(i,j) ]\\ 设sum(n, m) = \sum_{i = 1} ^ n \sum_{j = 1}^m i \ast j\\ = \frac{n \ast (n+1)}{2} \ast \frac{m*(m+1)}{2}\\ 那么对于F(n, m, k) 由反演我们可以得到\\ F(n, m, k) =\sum_{k | d}^{n} f(n, m, d)\\ f(n, m, k) = \sum_{k | d}^{n} F(n, m, d) \varphi(\frac{d}{k})\\ 那么结合一下上面 的公式，我们要求gcd(x, y) = 1的答案\\ f(n, m, 1) = \sum_{d = 1} ^n\varphi(d) \ast F(n, m, d)\\ = \sum_{d = 1} ^n \varphi(d) \ast d \ast d \ast sum(n/d,m/d)\\ 所以回到开始\\ ans = \sum_{d=1}^n d \ast f(n/d,m/d,1)\\ = \sum_{d = 1} ^n d \ast \sum_{k=1}^{n/d}\varphi(k) \ast k \ast k \ast sum(n/d/k,m/d/k)$

最后我们用两个分块来解决这个问题时间复杂度大概是O(n)

~~公式是真的难推，~~感觉头发都掉完了。

开始T到怀疑人生，后来发现这是单组数据，所以我们每次不要直接打1e7的表，根据输入数据打表才是最优解。

代码:

/*************************************************************************
    > File Name: bzoj-2154.cpp
    > Author: wood
    > Mail: cbcruel@gmail.com
    > Created Time: 2018年08月02日 星期四 10时30分53秒
 ************************************************************************/

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<queue>
#include<map>
#include<set>
#include<string>

using namespace std;
typedef long long LL;
#define sqr(x) (x)*(x)
const double eps = 1e-8;
const double C = 0.57721566490153286060651209;
const double pi = acos(-1.0);
const int mod = 20101009;
const int maxnn = 1e7 + 10;
int maxn;
int mu[maxnn],pri[maxnn], sum[maxnn];
bool vis[maxnn];
int cnt;

template <class T>
inline bool scan_d(T &ret) {
    char c; int sgn;
    if(c = getchar(),c == EOF) return 0; //EOF
    while(c != '-' && (c < '0' || c > '9')) c = getchar();
    sgn = (c == '-') ? -1 : 1;
    ret = (c == '-') ? 0 : (c-'0');
    while(c = getchar(),c >= '0' && c <= '9') ret = ret * 10 + (c - '0');
    ret *= sgn;
    return 1;
}
inline void out(int x) {
    if(x > 9) out( x / 10);
    putchar( x % 10 + '0');
}

void init(){
    sum[1] = mu[1] = 1;
    cnt = 0;
    for(int i = 2; i < maxn; ++i){
        if(!vis[i]){
            pri[++cnt] = i;
            mu[i] = -1; 
        }for(int j = 1; j <= cnt; ++j){
            if(i*pri[j] >= maxn) break;
            vis[i*pri[j]] = 1;
            if(i % pri[j] == 0){
                mu[i*pri[j]] = 0;
                break;
            }mu[i*pri[j]] = -mu[i];
        }sum[i] = (sum[i-1] + 1LL*i*i*mu[i])%mod;
        sum[i] = (sum[i] + mod)% mod;
    }
}

inline int cal(int x, int y){
    return (1LL * x*(x+1)/2 % mod) * (1LL * y *(y+1) /2%mod) % mod;
}

inline int solve(int n, int m){
    int ans = 0;
    for(int i = 1, j; i <= n; i=j+1){
        j = min(n/(n/i),m/(m/i));
        ans = (ans +1LL* (sum[j] - sum[i-1]) * cal(m/i,n/i) )% mod;
    }return (ans+mod)%mod;
}

int main(){
    int n, m;
    scan_d(n);
    scan_d(m);
    LL ans = 0;
    int tt;
    if(n > m) swap(n,m);
    maxn = n+10;
    init();
    for(int i = 1, j; i <= n; i = j+1){
        j = min((n/(n/i)), m / (m/i));
        ans = (ans +(i +j) *1LL* (j -i+1)/2 % mod * solve(n/i,m/i)) % mod;
    }out(ans);
    return 0;
}