題目描述

今天的數學課上，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只想知道表格裏所有數的和mod20101009的值。

輸入輸出格式

輸入格式：

輸入的第一行包含兩個正整數，分別表示N和M。

輸出格式：

輸出一個正整數，表示表格中所有數的和mod20101009的值。

輸入輸出樣例

輸入樣例#1：

4 5

輸出樣例#1：

122

說明

30%的數據滿足N, M≤ 10^3。

70%的數據滿足N, M≤ 10^5。

100%的數據滿足N, M≤ 10^7。

CTM啥破玩意兒

以上只是一只SB發的牢騷==

這題是反演題

然後我不會做==

首先題目要求的是\(\sum_{i = 1}^{n}\sum_{j=1}^{m}{\frac{i * j}{gcd(i,j)}}\)

先讓n < m

然後根據套路先考慮gcd

\(\sum_{d=1}^{n}\sum_{i=1}^{n}\sum_{j=1}^{m{\frac{i * j}{d}[gcd(i , j)== d]}}\)

然後繼續根據常規的套路我們再把那個d提出來

\(\sum_{d = 1}^{n}{d}\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}{i * j[gcd(i,j)==1]}\)

然後就可以設目標函數\(f(x) = \sum_{i=1}^{n/d}\sum_{j=1}^{m/d}{i * j[gcd(i,j)==x]}\)

然後再設\(g(x) = \sum_{x|t}{f(t)}\)

也就是\(g(x) = \sum_{i=1}^{n/d}\sum_{j=1}^{m/d}{i * j[x|gcd(i,j)]}\)

然後我們把這個x給提出來

\(g(x) = x^2 \sum_{i=1}^{n/dx}\sum_{j=1}^{m/dx}{i * j[1|gcd(i,j)]}\)

然後那個\([1|gcd(i,j)]\)沒啥用了

\(g(x)=x^2\sum_{i=1}^{n/dx}\sum_{j=1}^{m/dx}{i*j}\)

然後我們回到我們剛才化簡的答案式子

\(\sum_{d = 1}^{n}{d}\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}{i * j[gcd(i,j)==1]}\)

是不是可以寫成\(\sum_{d=1}{d}f(1)\)

然後我們對\(f(1)\)再反演一下

是不是就能得到\(\sum_{d=1}^{n}{d}\sum_{i=1}^{n}{μ(i)*g(i)}\)

然後我們再把這個\(g(i)\)展開

就得到了\(\sum_{d=1}^{n}{d}\sum_{i=1}^{n}\sum_{x=1}^{n/di}\sum_{y=1}^{m/di}{x*y}\)

然後我們需要預處理出\(\sum_{i=1}^{n}{μ(i)*i^2}\)

後面的那一坨是等差數列可以O(1)計算出來

所以可以兩次整除分塊O(n)求解

#include<cstdio>
#include<cstring>
#include<iostream>
#include<algorithm>
# define int long long
const int M = 10000005 ;
const int mod = 20101009 ;
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 ;
}
int mu[M] , p[M] , pnum ;
bool Notp[M] ;
int sum[M] , n , m ;
int fsum[M] , Ans ;
inline void Get_Mu() {
    Notp[1] = true ; mu[1] = 1 ;
    for(int i = 2 ; i <= n ; i ++) {
        if(!Notp[i]) p[++pnum] = i , mu[i] = -1 ;
        for(int j = 1 ; j <= pnum && i * p[j] <= n ; j ++) {
            Notp[i * p[j]] = true ;
            if(i % p[j] == 0) {
                mu[i * p[j]] = 0 ;
                break ;
            }
            mu[i * p[j]] = -mu[i] ;
        }
    }
    for(int i = 1 ; i <= n ; i ++) sum[i] = (sum[i - 1] + mu[i] + mod)%mod ;
}
inline int Calc(int n , int m) {
    int temp = 0 ;
    for(int l = 1 , r ; l <= n ; l = r + 1) {
        r = min(n / (n / l) , m / (m / l)) ;
        int val = ((1 + n / l) * (n / l) / 2 % mod)%mod * ((1 + m / l) * (m / l) / 2 % mod)%mod ;  
        temp = (temp + (fsum[r] - fsum[l - 1] + mod)%mod * val + mod)%mod ;
    }
    return (temp + mod)%mod ;
}
# undef int
int main() {
# define int long long
    n = read() ; m = read() ;
    if(n > m) swap(n , m) ;
    Get_Mu() ;
    for(int i = 1 ; i <= n ; i ++) fsum[i] = (fsum[i - 1] + mu[i] * i * i%mod + mod)%mod ;
    for(int l = 1 , r ; l <= n ; l = r + 1) {
        r = min(n / (n / l) , m / (m / l)) ;
        int val = ((l + r) * (r - l + 1LL) / 2LL + mod)%mod ;
        Ans = (Ans + Calc(n / l , m / l) * val%mod)%mod ;
    }
    printf("%lld\n",Ans) ;
    return 0 ;
}

[國家集訓隊]Crash的數字表格 / JZPTAB