bzoj 2154

收获：

　　1、当一个东西的取值范围很小时，或者感觉它很麻烦时，就枚举它

　　2、熟悉mobius函数、euler函数的和函数，以及euler函数用mobius函数的表示。

　　3、下取整分块理解更深了。

/**************************************************************
    Problem: 2154
    User: idy002
    Language: C++
    Result: Accepted
    Time:5584 ms
    Memory:157916 kb
****************************************************************/
 
#include <cstdio>
#include <iostream>
#define M 20101009
using namespace std;
 
typedef long long dnt;
 
int prm[100000], isnot[10000010], mu[10000010], ptot;
dnt mds[10000010];
 
void init( int n ) {
    mu[1] = 1;
    for( int i=2; i<=n; i++ ) {
        if( !isnot[i] ) {
            prm[++ptot] = i;
            mu[i] = -1;
        }
        for( int j=1; j<=ptot && i*prm[j]<=n; j++ ) {
            isnot[i*prm[j]] = true;
            if( i%prm[j]==0 ) {
                mu[i*prm[j]] = 0;
                break;
            }
            mu[i*prm[j]] = -mu[i];
        }
    }
    for( dnt d=1; d<=n; d++ ) {
        mds[d] = (mu[d]*d*d)%M;
        mds[d] += mds[d-1];
        if( mds[d]>=M ) mds[d]-=M;
        if( mds[d]<0 ) mds[d]+=M;
    }
}
 
inline dnt S( dnt x, dnt y ) {
    return (((1+x)*x/2%M)*(((1+y)*y)/2%M)%M);
}
dnt F( dnt x, dnt y ) {
    if( x>y ) swap(x,y);
    dnt rt = 0;
    for( dnt d=1; d<=x; d++ ) {
        dnt dd = min( x/(x/d), y/(y/d) );
        rt += S(x/d,y/d) * (mds[dd]-mds[d-1]) % M;
        if( rt<0 ) rt += M;
        if( rt>=M ) rt -= M;
        d = dd;
    }
    return rt;
}
dnt calc( dnt n, dnt m ) {
    if( n>m ) swap(n,m);
    dnt rt = 0;
    for( dnt d=1; d<=n; d++ ) {
        dnt dd=min( n/(n/d), m/(m/d) );
        rt += ((d+dd)*(dd-d+1)/2 % M) * F( n/d, m/d ) % M;
        if( rt<0 ) rt+=M;
        if( rt>=M ) rt-=M;
        d = dd;
    }
    return rt;
}
int main() {
    int n, m;
    scanf( "%d%d", &n, &m );
    if( n>m ) swap(n,m);
    init(n);
    printf( "%lld
", calc(n,m) );
}