Description
今天的数学课上，Crash小朋友学习了最小公倍数(Least Common Multiple)。对于两个正整数a和b，LCM(a, b)表示能同时被a和b整除的最小正整数。例如，LCM(6, 8) = 24。回到家后，Crash还在想着课上学的东西，为了研究最小公倍数，他画了一张N*M的表格。每个格子里写了一个数字，其中第i行第j列的那个格子里写着数为LCM(i, j)。一个4*5的表格如下： 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的值。
HINT
N, M ≤ 107。Solution
依旧暴力写开
依旧枚举gcd
再令
继续暴力写出F
再令
于是得到
愉快地推完了//跪PoPoQQQ神推
1 #include<map> 2 #include<cmath> 3 #include<ctime> 4 #include<queue> 5 #include<stack> 6 #include<cstdio> 7 #include<climits> 8 #include<iomanip> 9 #include<cstring> 10 #include<cstdlib> 11 #include<iostream> 12 #include<algorithm> 13 14 #define maxp 10000000 15 #define maxn 10000000+5 16 #define mod 20101009 17 #define set(a,b) memset(a,(b),sizeof(a)) 18 #define fr(i,a,b) for(ll i=(a),_end_=(b);i<=_end_;i++) 19 #define rf(i,b,a) for(ll i=(a),_end_=(b);i>=_end_;i--) 20 #define fe(i,a,b) for(int i=first[(b)],_end_=(a);i!=_end_;i=s[i].next) 21 #define fec(i,a,b) for(int &i=cur[(b)],_end_=(a);i!=_end_;i=s[i].next) 22 23 using namespace std; 24 25 typedef long long ll; 26 27 ll f[maxn]; 28 int prime[maxn],pri[maxn],miu[maxn],tot=0; 29 ll ans; 30 int n,m; 31 32 void read() 33 { 34 #ifndef ONLINE_JUDGE 35 freopen("2154.in","r",stdin); 36 freopen("2154.out","w",stdout); 37 #endif 38 cin >> n >> m ; 39 } 40 41 void write() 42 { 43 cout << (ans+mod)%mod ; 44 } 45 46 void get() 47 { 48 miu[1]=1; 49 fr(i,2,min(n,m)){ 50 if( !prime[i] ) pri[++tot]=i,miu[i]=-1; 51 int j=1; 52 while( j<=tot && pri[j]*i<=min(n,m) ){ 53 prime[ pri[j]*i ]=1; 54 if( i%pri[j]==0 ){ 55 miu[pri[j]*i]=0; 56 break; 57 } 58 miu[pri[j]*i]=-miu[i]; 59 j++; 60 } 61 } 62 fr(i,1,min(n,m)) 63 f[i]=(f[i-1]+i*i*miu[i]%mod)%mod; 64 } 65 66 ll Sum(ll x,ll y) 67 { 68 return ((x+1)*x/2%mod)*((y+1)*y/2%mod)%mod; 69 } 70 71 ll F(ll x,ll y) 72 { 73 ll res=0,i=1,pos; 74 while( i<=x ){ 75 pos=min(x/(x/i),y/(y/i)); 76 res=(res+(f[pos]-f[i-1])*Sum(x/i,y/i)%mod)%mod; 77 i=pos+1; 78 } 79 return res; 80 } 81 82 ll calc(ll x,ll y) 83 { 84 if( x>y ) swap(x,y); 85 ll res=0,i=1,pos; 86 while( i<=x ){ 87 pos=min(x/(x/i),y/(y/i)); 88 res=(res+(pos-i+1)*(pos+i)/2%mod*F(x/i,y/i)%mod)%mod; 89 i=pos+1; 90 } 91 return res; 92 } 93 94 void work() 95 { 96 get(); 97 ans=calc(n,m); 98 } 99 100 int main() 101 { 102 read(); 103 work(); 104 write(); 105 return 0; 106 }