BZOJ 3529 [Sdoi2014]数表

#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long ll;
#define enter putchar('
')
#define space putchar(' ')
template <class T>
void read(T &x){
    char c;
    bool op = 0;
    while(c = getchar(), c > '9' || c < '0')
	if(c == '-') op = 1;
    x = c - '0';
    while(c = getchar(), c >= '0' && c <= '9')
	x = x * 10 + c - '0';
    if(op) x = -x;
}
template <class T>
void write(T x){
    if(x < 0) putchar('-'), x = -x;
    if(x >= 10) write(x / 10);
    putchar('0' + x % 10);
}
const int N = 10000000, P = 20101009;
int n, m, mu[N + 5], prime[N + 5], tot;
ll sum[N + 5], ans;
bool notprime[N + 5];
#define mod(x) (((x) % P + P) % P)
void init(){
    mu[1] = 1;
    for(int i = 2; i <= m; i++){
	if(!notprime[i]) prime[++tot] = i, mu[i] = -1;
	for(int j = 1; j <= tot && prime[j] * i <= m; j++){
	    notprime[i * prime[j]] = 1;
	    if(i % prime[j]) mu[i * prime[j]] = -mu[i];
	    else{
		mu[i * prime[j]] = 0;
		break;
	    }
	}
    }
    for(ll i = 1; i <= m; i++)
	sum[i] = (sum[i - 1] + i * i * mu[i] % P) % P;
}

ll G(ll x, ll y){
    return (x * (x + 1) / 2 % P) * (y * (y + 1) / 2 % P) % P;
}
ll F(ll x, ll y){
    ll ret = 0;
    if(x > y) swap(x, y);
    for(int i = 1, last; i <= x; i = last + 1){
	last = min(x / (x / i), y / (y / i));
	ret = (ret + G(x / i, y / i) * (sum[last] - sum[i - 1]) % P) % P;
    }
    return ret;
}
int main(){
    read(n), read(m);
    if(n > m) swap(n, m);
    init();
    for(ll i = 1, last; i <= n; i = last + 1){
	last = min(n / (n / i), m / (m / i));
	ans = (ans + (i + last) * (last - i + 1) / 2 % P * F(n / i, m / i) % P) % P;
    }
    write((ans + P) % P), enter;
    return 0;
}