【BZOJ 2154】【JZOJ 1938】【2011集训队出题】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的表格如下：

看着这个表格，Crash想到了很多可以思考的问题。不过他最想解决的问题却是一个十分简单的问题：这个表格中所有数的和是多少。当N和M很大时，Crash就束手无策了，因此他找到了聪明的你用程序帮他解决这个问题。由于最终结果可能会很大，Crash只想知道表格里所有数的和mod 20101009的值。

Solution

抽象题意：求：

$$\sum_{i=1}^n\sum_{j=1}^m\frac{i*j}{\gcd(i,j)}$$

用一个简单的繁衍，
设

$$f(d)=\sum_{i=1}^n\sum_{j=1}^m ij[gcd(i,j)=d]$$
同时设： $$g(d)=\sum_{i=1}^{\frac{n}{d}} f(id)$$
我们发现g很好求： $$g(d)=d^2*\frac{\lfloor\frac{n}{d}\rfloor\lfloor\frac{n}{d}+1\rfloor}{2}*\frac{\lfloor\frac{m}{d}\rfloor\lfloor\frac{m}{d}+1\rfloor}{2}$$
于是：
$$f(d)=\frac{1}{4}d^2*\sum_{i=1}^{\frac{n}{d}}\mu(i)i^2*\lfloor\frac{n}{id}\rfloor\lfloor\frac{n}{id}+1\rfloor\lfloor\frac{m}{id}\rfloor\lfloor\frac{m}{id}+1\rfloor$$
用分块搞，
$$ans=\sum_{d=1}^n \frac{f(d)}{d}$$
也用分块搞，

复杂度：$O(n)$

Code

有点小慢

#include<iostream>
#include<cstdio>
#include<cstdlib>
#define fo(i,a,b) for(int i=a;i<=b;i++)
using namespace std;
typedef long long LL;
const int N=1e7+10,mo=20101009;
const LL ny2=10050505,ny4=15075757;
int n,m;
int pr[N],mu[N];
bool prz[N];
int f[N];
LL ans;
void pre(int n)
{
    mu[1]=1;
    fo(i,2,n)
    {
        if(!prz[i])pr[++pr[0]]=i,mu[i]=-1;
        fo(j,1,pr[0])
        {
            int t=i*pr[j];
            if(t>n)break;
            prz[t]=1;
            if(i%pr[j]==0){mu[t]=0;break;}
            mu[t]=-mu[i];
        }
    }
    fo(i,1,n)f[i]=((LL)mu[i]*i%mo*i%mo+f[i-1])%mo;
}
LL fk(int d)
{
    LL ans=0;
    LL l=1,nd=n/d,md=m/d;
    while(l<=nd)
    {
        LL r=min(nd/(nd/l),md/(md/l));
        LL t=(nd/l)*(nd/l+1)%mo*(md/l)%mo*(md/l+1)%mo;
        ans=(ans+t*(f[r]-f[l-1])%mo)%mo;
        l=r+1;
    }
    return ans*ny4%mo;
} 
int main()
{
    scanf("%d%d",&n,&m);
    if(n>m)swap(n,m);
    pre(n);
    LL l=1;
    while(l<=n)
    {
        LL r=min(n/(n/l),m/(m/l));
        ans=((LL)ans+fk(l)*((l+r)*(r-l+1)%mo*ny2%mo)%mo)%mo;
        l=r+1;
    }
    printf("%lld\n",(ans+mo)%mo);
    return 0;
}