bzoj2154 Crash的数字表格

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

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

Solution

---

留下来的坑得一个一个填-_-||
答案要我们求这个

$$\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j)$$
只会做和gcd有关的东西，于是稍作变换 $$ans=\sum_{i=1}^{n}\sum_{j=1}^{m}\frac{ij}{gcd(i,j)}$$
套路一波枚举gcd得到 $$ans=\sum_{d=1}^{n}\frac{1}{d}\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}di*dj*[gcd(i,j)=1]$$
看到辣个=1就开心，强行套μ有 $$ans=\sum_{d=1}^{n}d\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}ij\sum_{x|gcd(i,j)}\mu(x)$$
把x拉出来，注意这里i、j枚举的是dx的倍数，那么后面实际要多乘上两个x $$ans=\sum_{d=1}^{n}d\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\mu(x)\sum_{i=1}^{\lfloor\frac{n}{dx}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{dx}\rfloor}i*j*x^2$$
设 $$S(x)=\sum_{i=1}^{n}i$$
那么 $$ans=\sum_{d=1}^{n}d\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\mu(x)x^2S(\lfloor\frac{n}{dx}\rfloor)S(\lfloor\frac{m}{dx}\rfloor)$$
这里令$T=dx$然后枚举T
$$ans=\sum_{T=1}^{n}S(\lfloor\frac{n}{T}\rfloor)S(\lfloor\frac{m}{T}\rfloor)T\sum_{x|T}\mu(x) x$$
这里的的$\sum_{x|T}\mu(x)x$虽然和狄利克雷卷积长得很像但是并不是狄利克雷卷积但是并不影响它是积性函数（绕

今天终于完全搞懂线性筛+筛积性函数的真谛惹（弱）
筛积性函数主要是得求出$f(p^k)$p是质数这样的函数值，然后根据线性筛里面的pj和i的关系分类搞一搞就OK了

一开始跑得巨慢是每次都跑了N

Code

---

#include <stdio.h>
#include <string.h>
#include <algorithm>
#define rep(i,st,ed) for (int i=st;i<=ed;++i)

typedef long long LL;
const LL ny2=10050505;
const int MOD=20101009;
const int N=10000005;

int prime[N/10],low[N+1];
short mu[N+1];
bool not_prime[N+1];
LL f[N+1];

void pre_work(int n) {
    mu[1]=f[1]=low[1]=1;
    rep(i,2,n) {
        if (!not_prime[i]) {
            prime[++prime[0]]=i;
            low[i]=i;
            mu[i]=-1;
            f[i]=1-i;
        }
        for (int j=1;i*prime[j]<=n&&j<=prime[0];j++) {
            LL x=i*prime[j];
            not_prime[x]=1;
            if (i%prime[j]==0) {
                low[x]=low[i]*prime[j];
                mu[x]=0;
                if (i==low[i]) {
                    f[x]=f[i];
                } else {
                    f[x]=(f[i/low[i]]*f[low[i]*prime[j]])%MOD;
                }
                break;
            }
            low[x]=prime[j];
            f[x]=(f[i]*f[prime[j]])%MOD;
            mu[x]=-mu[i];
        }
    }
}

LL S(LL n) {
    LL ret=n*(n+1)%MOD*ny2%MOD;
    return ret;
}

void solve(LL n,LL m) {
    LL ans=0;
    rep(i,1,n) f[i]=f[i]*i%MOD;
    rep(i,1,n) f[i]=(f[i]+f[i-1])%MOD;
    for (LL i=1,j;i<=n;i=j+1) {
        j=std:: min(n/(n/i),m/(m/i));
        LL sum=(f[j]-f[i-1]+MOD)%MOD;
        LL a=S(n/i); LL b=S(m/i);
        LL tot=a*b%MOD*sum%MOD;
        ans=(ans+tot)%MOD;
    }
    if (ans<0) ans=(ans+MOD)%MOD;
    printf("%lld\n", ans);
}

int main(void) {
    LL n,m; scanf("%lld%lld",&n,&m);
    if (m<n) std:: swap(n,m);
    pre_work(n);
    solve(n,m);
    return 0;
}