《洛谷P1829 [国家集训队]Crash的数字表格 / JZPTAB》

第一次差不多算自己推出来的吧，代码也写出来了。

但是块的右边界写错了，真有我的.jpg。

看了一眼题，我就得出了$ans = sum_{i = 1}^{n}sum_{j = 1}^{m} i * j / gcd(i,j)$

然后就开始愉快（bushi）地推式子了。

$ans = sum_{i = 1}^{n}sum_{j = 1}^{m}sum_{d = 1}^{min(n,m)}i * j / d * |~d~|~gcd(frac{i}{d},frac{j}{d})~| = sum_{d = 1}^{min(n,m)}sum_{i = 1}^{[frac{n}{d}]}sum_{j = 1}^{[frac{m}{d}]}(id) * (jd) / d * |t | gcd(i,j) = 1|$

$ =  sum_{d = 1}^{min(n,m)}sum_{i = 1}^{[frac{n}{d}]}sum_{j = 1}^{[frac{m}{d}]}i * j * d * sum_{t | gcd(i,j)}^{} mu (t)$

$ = sum_{d = 1}^{min(n,m)} d sum_{t = 1}^{min(n / d,m /d))} mu (t)  sum_{i = 1}^{[frac{n}{dt}]}sum_{j = 1}^{[frac{m}{dt}]}it * jt$

$ =  sum_{d = 1}^{min(n,m)} d sum_{t = 1}^{min(n / d,m /d))} mu (t) * t^{2}sum_{i = 1}^{[frac{n}{dt}]}i sum_{j = 1}^{[frac{m}{dt}]}j$

后面的式子显然可以整除分块处理，然后也要预处理一下关于t的项。

然后一开始想的枚举d，但是这样肯定超时，后面仔细看后面的上限是n / d,m / d，那么可以再套一个整除分块。

这样的话就是最高复杂度也就是预处理：O(max(n,m)）。

细节就是与取模要防负数。

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
typedef pair<int,int> pii;
const int N = 1e7+5;
const int M = 1e7+5;
const LL Mod = 20101009;
#define pi acos(-1)
#define INF 1e9
#define CT0 cin.tie(0),cout.tie(0)
#define IO ios::sync_with_stdio(false)
#define dbg(ax) cout << "now this num is " << ax << endl;
namespace FASTIO{
    inline int read()
    {
        int x = 0,f = 1;char c = getchar();
        while(c < '0' || c > '9'){if(c == '-')f = -1;c = getchar();}
        while(c >= '0' && c <= '9'){x = (x << 1) + (x << 3) + (c ^ 48);c = getchar();}
        return x * f;
    }
    void print(int x){
        if(x < 0){x = -x;putchar('-');}
        if(x > 9) print(x/10);
        putchar(x%10+'0');
    }
}
using namespace FASTIO;

bool vis[N];
int mu[N],prime[N],tot = 0;
LL sum[N],f[N];
void init()
{
    mu[1] = 1;
    for(int i = 2;i < N;++i)
    {
        if(!vis[i])
        {
            prime[++tot] = i;
            mu[i] = -1;
        }
        for(int j = 1;j <= tot && i * prime[j] < N;++j)
        {
            vis[i * prime[j]] = 1;
            if(i % prime[j] == 0) break;
            else mu[i * prime[j]] = -mu[i];
        }
    }
    for(int i = 1;i < N;++i) sum[i] = (sum[i - 1] + (mu[i] * ((1LL * i * i) % Mod) % Mod) + Mod) % Mod;
    for(int i = 1;i < N;++i) f[i] = (f[i - 1] + i) % Mod;
}
LL solve(int x,int y)
{
    LL ans = 0;
    for(int L = 1,r = 0;L <= min(x,y);L = r + 1)
    {
        r = min(x / (x / L),y / (y / L));
        ans = (ans + ((sum[r] - sum[L - 1]) % Mod + Mod) % Mod * f[x / L] % Mod * f[y / L] % Mod) %Mod;
    }
    return ans;
}
int main()
{
    init();
    int n,m;n = read(),m = read();
    LL ans = 0;
    for(int L = 1,r = 0;L <= min(n,m);L = r + 1)
    {
        r = min(n / (n / L),m / (m / L));
        LL ma = solve(n / L,m / L);
        ans = (ans + ((f[r] - f[L - 1]) %Mod + Mod) % Mod * ma % Mod) % Mod;
    }
    printf("%lld
",ans);
    system("pause");
    return 0;
}