【BZOJ 2154】Crash的数字表格

我们定义$sum(n,m)=\sum_{i=1}^{n}\sum_{j=1}^{m}i*j$
$\sum_{i=1}^{n}\sum_{j=1}^{m} lcm(i,j)=\sum_{i=1}^{n}\sum_{j=1}^{m}{i*j\over gcd(i,j)}$
$=\sum_{i=1}^{\lfloor {n\over d}\rfloor}\sum_{j=1}^{\lfloor {m\over d}\rfloor} {di*dj \over d}e(gcd(i,j))$
$=\sum_{d=1}^{min(n,m)}d\sum_{i=1}^{\lfloor {n\over d}\rfloor}\sum_{j=1}^{\lfloor {m\over d}\rfloor} i*j e(gcd(i,j))$
$=\sum_{d=1}^{min(n,m)}d \sum_{k=1}^{min(\lfloor {n\over d} \rfloor,\lfloor {m\over d} \rfloor)}\mu(k)*k^2*sum(\lfloor {\lfloor {n\over d} \rfloor \over k} \rfloor,\lfloor {\lfloor {m\over d} \rfloor \over k} \rfloor)$
由于$\lfloor {n\over d} \rfloor$与$\lfloor {\lfloor {m\over d} \rfloor \over k} \rfloor$各有$\sqrt n$与$\sqrt\lfloor {n\over d} \rfloor$个取值，因而我们可以分块统计，code：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#define P 20101009LL
using namespace std;
long long n,m,tot=0,mu[10000001],pri[700001];
long long pre[10000001],sum[10000001];
bool f[10000001];
long long get(long long n,long long m)
{
    long long i,last=0; long long ans=0;
    if (n>m) swap(n,m);
    for (i=1;i<=n;i=last+1)
      {
        last=min(n/(n/i),m/(m/i));
        ans=(ans+(long long)(((long long)((sum[last]-sum[i-1]+P)%P)*(long long)(((long long)((((m/i)*(m/i+1))/2)%P)*(long long)((((n/i)*(n/i+1))/2)%P))%P)))%P)%P;
        if (ans<0) cout<<'!'<<endl;
      }
    return ans;
}
long long work(long long n,long long m)
{
    long long i,last=0; long long ans=0;
    if (n>m) swap(n,m);
    for (i=1;i<=n;i=last+1)
      {
        last=min(n/(n/i),m/(m/i));
        ans=(ans+(long long)((get(n/i,m/i)%P*(long long)((pre[last]-pre[i-1]+P)%P)))%P)%P;
        if (ans<0) cout<<'!'<<endl;
      }
    return ans;
}
int main()
{
    long long i,j,maxn;
    scanf("%I64d%I64d",&n,&m);
    mu[1]=1;maxn=max(n,m);
    for (i=2;i<=maxn;++i)
      {
        if (!f[i])
          {
            mu[i]=-1;
            pri[++tot]=i;
          }
        for (j=1;j<=tot&&i*pri[j]<=maxn;++j)
          {
            f[i*pri[j]]=true;
            if (i%pri[j]==0)
              {
                mu[i*pri[j]]=0;
                break;
              }
            else mu[i*pri[j]]=-mu[i];
          }
      }
    for (i=1;i<=maxn;++i) {sum[i]=(sum[i-1]+(long long)((i*i)%P*mu[i])%P+P)%P; pre[i]=(long long)(pre[i-1]+i)%P;}
    printf("%I64d\n",work(n,m));
}