bzoj2154 Crash的数字表格

标签： 数学 莫比乌斯反演

---

题目描述

今天的数学课上，Crash小朋友学习了最小公倍数(Least Common Multiple)。对于两个正整数a和b，LCM(a, b)表示能同时被a和b整除的最小正整数。例如，LCM(6, 8) = 24。回到家后，Crash还在想着课上学的东西，为了研究最小公倍数，他画了一张NM的表格。每个格子里写了一个数字，其中第i行第j列的那个格子里写着数为LCM(i, j)。一个45的表格如下： 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

数据范围

N，M <= 10^7

题解

原题意相当与求下列式子:
$$\sum_{i=1}^{N}\sum_{j=1}^{M}lcm(i,j)$$
当然存在一个很简单的式子:$lcm(i,j)=\frac{ij}{gcd(i,j)}$
所以原式等于:$$\sum_{i=1}^{N}\sum_{j=1}^{M}\frac{ij}{gcd(i,j)}$$
然后开始更改枚举项（当时我差点傻逼的直接开始设函数了），枚举$gcd(i,j)$
$$\sum_{i = 1}^{min(n,m)}\sum_{x=1}^{\frac{N}{x}}\sum_{y=1}^{\frac{M}{y}}\frac{xy}{i}i^{2}[gcd(x,y)==1]=\sum_{i = 1}^{min(n,m)}i\sum_{x=1}^{\frac{N}{x}}\sum_{y=1}^{\frac{M}{y}}xy[gcd(x,y)==1]$$
到这里我们就可以开始设函数了:
$$f(m)=\sum_{x=1}^{\frac{N}{x}}\sum_{y=1}^{\frac{M}{y}}xy[gcd(x,y)==m]$$
看着变量有点乱，重新整理一下:
$$f(x)=\sum_{i=1}^{n}\sum_{j=1}^{m}ij[gcd(i,j)==x]$$
$$F(x)=\sum_{x|d}f(d)=\sum_{i=1}^{n}\sum_{j=1}^{m}ij[x|gcd(i,j)]$$
到这里肯定就是化简$F(x)$了，考虑直接枚举x的倍数:
$$F(x)=\sum_{i=1}^{\frac{n}{x}}\sum_{j=1}^{\frac{m}{x}}x^{2}ij=x^{2}\sum_{i=1}^{\frac{n}{x}}\sum_{j=1}^{\frac{m}{x}}ij$$
算到这里$F(x)$肯定能$O(1)$求出来了（就相当于一个等差数列求和嘛）
那么回到所求答案:
$$Ans=\sum_{p=1}^{min(n,m)}p*f(1)$$

$$Ans=\sum_{p=1}^{min(n,m)}p*\sum_{d=1}^{min(\frac{n}{p},\frac{m}{p})}u(d)F(d)$$

对于第一层的$\sum$我们可以用整除分块的（因为第二层$\sum$有下整除运算）
接着展开$F(d)$:
$$Ans=\sum_{p=1}^{min(n,m)}p*\sum_{d=1}^{min(\frac{n}{p},\frac{m}{p})}u(d)d^{2}\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}}ij$$

我们发现对于第二层$\sum $我们又可以用整除分块，因此复杂度由$O(n^2)$优化到$O(N*大常数)$

那么此题解法这里采用的是单变量函数$f(x)$，当然这里也有设双变量函数$f(x,y)$.

代码

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#define N 10000011
#define ll long long 
#define p 20101009
using namespace std;

int cnt = 0;
int pri[N], u[N];
ll sumi[N];
bool vis[N];

inline int minn(int a, int b) {return a < b ? a : b;}

inline void get(int t)
{
    sumi[0] = 0; u[1] = 1;
    for(int i = 2; i <= t; i ++)
    {
        if(!vis[i]) pri[++ cnt] = i, u[i] = -1;
        for(int j = 1; j <= cnt && i * pri[j] <= t; j ++)
        {
            vis[i * pri[j]] = true;
            if(i % pri[j]) u[i * pri[j]] = -u[i];
            else break;
        }
    }
    for(int i = 1; i <= t; i ++) sumi[i] = (sumi[i - 1] + 1ll * u[i] * i % p * i % p) % p;
}

inline ll F(ll x, ll y)
{
    x = (x * (x + 1) / 2) % p;
    y = (y * (y + 1) / 2) % p;
    return x * y % p;
}

inline ll calc(int x, int y)
{
    int d = minn(x, y); ll ans = 0;
    for(int l = 1, r; l <= d; l = r + 1)
    {
        r = minn(x / (x / l), y / (y / l));
        ans = (ans + 1ll * (sumi[r] - sumi[l - 1] + p) % p * F(1ll * x / l, 1ll * y / l)) % p;
    }
    return ans % p;
}

int main()
{
    int n, m;
    scanf("%d%d", &n, &m);
    get(n);
    int d = minn(n, m); ll ans = 0;
    for(int l = 1, r; l <= d; l = r + 1)
    {
        r = minn(n / (n / l), m / (m / l));
        ans = (ans + 1ll * (1ll * (l + r) * (r - l + 1) / 2 % p) * calc(n / l, m / l) % p) % p;
    }
    printf("%lld", ans);
    return 0;
} 

心得

表示自己又傻逼的在设函数那里崩掉了，当时化简原式子后直接设:$$f(x)=\sum_{i = 1}^{min(n,m)}i\sum_{x=1}^{\frac{N}{x}}\sum_{y=1}^{\frac{M}{y}}xy[gcd(x,y)==x]$$
我说tm怎么化简不了…
表示双变量函数解法也不错的，多积累多积累
还有温馨提示此题如果模数处理不好WA到哭，所以每次乘了以后都膜一膜吧…
最后温馨提示是先读入n,m后在根据数据大小处理u函数和其他前缀和，不要先直接处理个一千万的再读入（这里没有多组数据）否则你会T到屎…（我因为这个调了1h….o(╥﹏╥)o）