[國家集訓隊]Crash的數字表格 / JZPTAB
題目描述
今天的數學課上,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