bzoj 2154
收获：
1、当一个东西的取值范围很小时，或者感觉它很麻烦时，就枚举它
2、熟悉mobius函数、euler函数的和函数，以及euler函数用mobius函数的表示。
3、下取整分块理解更深了。
1 /************************************************************** 2 Problem: 2154 3 User: idy002 4 Language: C++ 5 Result: Accepted 6 Time:5584 ms 7 Memory:157916 kb 8 ****************************************************************/ 9 10 #include <cstdio> 11 #include <iostream> 12 #define M 20101009 13 using namespace std; 14 15 typedef long long dnt; 16 17 int prm[100000], isnot[10000010], mu[10000010], ptot; 18 dnt mds[10000010]; 19 20 void init( int n ) { 21 mu[1] = 1; 22 for( int i=2; i<=n; i++ ) { 23 if( !isnot[i] ) { 24 prm[++ptot] = i; 25 mu[i] = -1; 26 } 27 for( int j=1; j<=ptot && i*prm[j]<=n; j++ ) { 28 isnot[i*prm[j]] = true; 29 if( i%prm[j]==0 ) { 30 mu[i*prm[j]] = 0; 31 break; 32 } 33 mu[i*prm[j]] = -mu[i]; 34 } 35 } 36 for( dnt d=1; d<=n; d++ ) { 37 mds[d] = (mu[d]*d*d)%M; 38 mds[d] += mds[d-1]; 39 if( mds[d]>=M ) mds[d]-=M; 40 if( mds[d]<0 ) mds[d]+=M; 41 } 42 } 43 44 inline dnt S( dnt x, dnt y ) { 45 return (((1+x)*x/2%M)*(((1+y)*y)/2%M)%M); 46 } 47 dnt F( dnt x, dnt y ) { 48 if( x>y ) swap(x,y); 49 dnt rt = 0; 50 for( dnt d=1; d<=x; d++ ) { 51 dnt dd = min( x/(x/d), y/(y/d) ); 52 rt += S(x/d,y/d) * (mds[dd]-mds[d-1]) % M; 53 if( rt<0 ) rt += M; 54 if( rt>=M ) rt -= M; 55 d = dd; 56 } 57 return rt; 58 } 59 dnt calc( dnt n, dnt m ) { 60 if( n>m ) swap(n,m); 61 dnt rt = 0; 62 for( dnt d=1; d<=n; d++ ) { 63 dnt dd=min( n/(n/d), m/(m/d) ); 64 rt += ((d+dd)*(dd-d+1)/2 % M) * F( n/d, m/d ) % M; 65 if( rt<0 ) rt+=M; 66 if( rt>=M ) rt-=M; 67 d = dd; 68 } 69 return rt; 70 } 71 int main() { 72 int n, m; 73 scanf( "%d%d", &n, &m ); 74 if( n>m ) swap(n,m); 75 init(n); 76 printf( "%lld\n", calc(n,m) ); 77 }
