BZOJ2154: Crash的数字表格

BZOJ2154

题目要求的是

$$\sum_{i=1}^n\sum_{j=1}^m\frac{i*j}{gcd(i,j)}$$
那么我们枚举 $gcd(i,j)$，不妨令 $n<m$
原式变成了：
$$\sum_{d=1}^n\frac{\sum_{i=1}^n\sum_{j=1}^mi*j[gcd(i,j)==d]}{d}$$
$$=\sum_{d=1}^n\frac{\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}}i*j*d^2[gcd(i,j)==1]}{d}$$
$$=\sum_{d=1}^nd*{\sum_{i=1}^{\frac{n}{d}}}\sum_{j=1}^{\frac{m}{d}}i*j[gcd(i,j)==1]$$

定义$F(x,y)=\sum_{i=1}^n\sum_{j=1}^mi*j[gcd(i,j)==1]$
那么原式$=\sum_{d=1}^nd*F(\frac{n}{d},\frac{m}{d})$
定义$S(x,y)=\sum_{i=1}^x\sum_{j=1}^yi*j$
$G(x,y,z)=\sum_{i=1}^x\sum_{j=1}^y\sum_{z|gcd(i,j)}i*j=z^2*S(\frac{x}{z},\frac{y}{z})$

$$F(x,y)=G(x,y,1)-(G(x,y,2)+G(x,y,3)+G(x,y,5)+……)+(G(x,y,6)+(G(x,y,10)……)-(G(x,y,30)+……)$$
$$F(x,y)=\sum_{t=1}^xμ(t)t^2S(\frac{x}{t},\frac{y}{t})$$
那么 $$ans=\sum_{d=1}^nd*\sum_{t=1}^xμ(t)t^2S(\frac{n}{td},\frac{m}{td})$$
令 $T=td$
$$ans=\sum_{T=1}^nS(n/T,m/T)\sum_{d|T}\frac{T}{d}μ(d)d^2$$

令$f(T)=\sum_{d|T}\frac{T}{d}μ(d)d^2$
$f(T)$是积性函数，可以线性筛。暴力枚举好像不行。。反正我自己试了T掉了。

【代码】

#include <cstdio>
#include <iostream>
#include <algorithm>
#define N 10000005
#define INF 0x7fffffff
#define mod 20101009
using namespace std;
typedef long long ll;
typedef pair<int,int> pa;
int read()
{
    int x=0,f=1;char ch=getchar();
    while(!isdigit(ch)){if(ch=='-') f=-1;ch=getchar();}
    while(isdigit(ch)){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
    return x*f;
}

int n,m;
int Miu[N],p[N];
ll sum[N];
bool Not_Prime[N];

void Get_Miu()
{
    Miu[1]=sum[1]=1;
    for(register int i=2;i<=m;i++)
    {
        if(!Not_Prime[i]) p[++p[0]]=i,Miu[i]=-1,sum[i]=1LL*i*(1-i);
        for(register int j=1;j<=p[0]&&i*p[j]<=m;j++)
        {
            Not_Prime[i*p[j]]=1;
            if(i%p[j]!=0) Miu[i*p[j]]=-Miu[i],sum[i*p[j]]=sum[i]*sum[p[j]];
            else {sum[i*p[j]]=sum[i]*p[j];break;}
        }
    }
    for(register int i=1;i<=n;i++) sum[i]=(sum[i-1]+sum[i]+mod)%mod;
}

ll Get_Sum(int x,int y){
    return (1LL*x*(x+1)/2)%mod*((1LL*y*(y+1)/2)%mod)%mod;
}

int main()
{
    n=read(),m=read();if(n>m) swap(n,m);
    Get_Miu();
    int pos;ll ans=0;
    for(register int i=1;i<=n;i=pos+1)
    {
        pos=min(n/(n/i),m/(m/i));
        ans=(ans+Get_Sum(n/i,m/i)*(sum[pos]-sum[i-1]+mod))%mod;
        ans=(ans+mod)%mod;
    }
    printf("%lld\n",ans);
    return 0;
}