题解 BZOJ-2154 莫比乌斯反演

2154: Crash的数字表格

Time Limit: 20 Sec Memory Limit: 259 MB

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的值。

Input

输入的第一行包含两个正整数，分别表示N和M。

Output

输出一个正整数，表示表格中所有数的和mod 20101009的值。

Sample Input

4 5

Sample Output

122

【数据规模和约定】

100%的数据满足N, M ≤ 10^7。

WA了因为#define long long int搜索替换的时候把LL换成了int，失策。

大概就是给定n, m, 求

具体看PPT第23页。

http://pan.baidu.com/s/1dFwyR7b

#include "cstdio"
#include "cctype"
#include "cstring"

#define max(a, b)  ((a) > (b) ? (a) : (b))
#define min(a, b)  ((a) < (b) ? (a) : (b))

#define rep(i, a, b)  for(register long long i = a; i <= b; ++i)
#define res(i, a, b)  for(register long long i = a; i >= b; --i)

template <class T>
inline bool readIn(T &x)  {
    T flag = 1;  char ch;
    while(!(isdigit(ch = (char) getchar())) && ch != EOF)  if( ch == '-' )  flag = -1;
    if(ch == EOF)  return false;
    for(x = ch - 48; isdigit(ch = (char) getchar()); x = (x << 1) + (x << 3) + ch - 48);
    x *= flag;
    return true;
}

template <class T>
inline void write(T x)  {
     if (x > 9)
       write(x / 10);
     putchar(x % 10 + 48);
}

template <class T>
inline bool writeIn(T x)  {
     if (x < 0)  {
         putchar('-');
         x = -x;
     }
     write(x);
     return true;
}

const long long MAXN = (long long) 10000005, MOD = 20101009;

long long mu[MAXN], prime[MAXN], s[MAXN], tot, m, n, t, ans;
bool mark[MAXN];

inline long long sum(long long x, long long y)  {   
    return ((x * (x + 1) >> 1) % MOD) * (( y * (y + 1) >> 1) % MOD) % MOD;
}

inline void init()  {
    long long t = min(n, m);
    mu[1] = 1;
    rep(i, 2, t)  {
        if( !mark[i] )  
            prime[++tot] = i, mu[i] = -1;
        for( long long j = 1; j <= tot && prime[j] * i <= t; ++j )  {
            mark[prime[j] * i] = true;
            mu[prime[j] * i] = -mu[i];
            if( i % prime[j] == 0 )  {mu[prime[j] * i] = 0; break;      }
        }
    }
    rep(i, 1, t)  
      s[i] = (s[i - 1] + ( i * i * mu[i] ) % MOD) % MOD;
}

inline long long F(long long x, long long y)  { 
    long long ans = 0, t= min(x, y);
    for(register long long i = 1, last; i <= t; i = last + 1)   
        last = min(x / (x / i), y / (y / i)),
        ans = (ans + (s[last] - s[i-1]) * sum(x / i, y / i) % MOD )% MOD;   
    return ans;
}

int main()  {
    readIn(n);readIn(m);
    init();
    for(register long long d=1, last;d<=min(n,m);d=last+1)
        last=min(n/(n/d),m/(m/d)),
        ans=(ans + ((d + last) * (last - d + 1) >> 1) % MOD * F(n / d, m / d) % MOD) % MOD;
    writeIn((ans + MOD) % MOD);
    putchar('\n');
    return 0;
}