[bzoj2154]crash的数字表格 解题报告

借这题理解了tangjz的例题。。这题做法与那道题基本相同。
先来看一下最普通的做法：（以下均设n≤m）

∑i=1n∑j=1mlcm(i,j)

=∑i=1n∑j=1mij(i,j)

=∑g=1n1g∑i=1⌊ng⌋∑i=1⌊mg⌋ij∑dμ(d)[d|i][d|j]

(枚举gcd)

=∑g=1n1g∑d=1⌊ng⌋μ(d)(gd)2⌊ndg⌋⌊mdg⌋

=∑i=1n⌊ni⌋⌊mi⌋i2∑d|iμ(d)id

(将g,d卷积)

=∑i=1n⌊ni⌋⌊mi⌋i∑d|iμ(d)d

注意到f(n)=n∑d|nμ(d)d显然是一个积性函数，考虑n=∏ki=1paii，则f(n)=n∏ki=1(1−pi).所以就可以很容易的线筛出来了。然后记一下f(n)的前缀和，便可以分块统计了。所以对于bzoj2693（bzoj2154的多组询问版）便可以做到O(n+Tn√)
但是对于单组询问，既然后面的那个显然是积性函数，而且我们只要它的前缀和，那它是否应该可以杜教筛出来呢？
我一开始是想对这玩意儿直接做积性前缀和和普通前缀和。。然后我发现完全做不了。但是正如tangjz给的那个例题一样，应该先将式子化简。

∑i=1ni∑d|iμ(d)d

=∑i=1n∑d|iμ(d)d2id

=∑i=1ni∑j=1⌊ni⌋μ(d)d2

所以问题就变成了求f(n)=μ(n)n2的前缀和，这时我们再用杜教筛，将它与f(n)=n^2卷积。

∑i=1n∑d|iμ(d)d2(id)2

=∑i=1ni2∑d|iμ(d)

=1

=∑i=1ni2∑j=1⌊ni⌋μ(d)d2

这样就可以杜教筛了！
但是还有一个问题，我们求出了∑ni=1μ(i)i2，可我们要求∑ni=1i∑⌊ni⌋j=1μ(j)j2，所以我们还需要预处理出n23的这玩意儿，然后就可以在O(n23)的时间复杂度搞出来了。
普通版：

#include<cstdio>#include<iostream>using namespace std;#define Mod 20101009typedef long long LL;LL pow(int a,int x){    LL ans=1,prod=a;    for(;x;x>>=1,prod=prod*prod%Mod)        if(x&1)            ans=ans*prod%Mod;    return ans;}LL sf[10000005];int prime[1000005];bool p[10000005];int main(){    int i,j;    sf[1]=1;    for(i=2;i<=10000000;++i){        if(!p[i])sf[i]=1-(prime[++prime[0]]=i);        for(j=1;j<=prime[0]&&i*prime[j]<=10000000;++j){            p[i*prime[j]]=1;            if(i%prime[j])sf[i*prime[j]]=sf[i]*(1-prime[j]);            else{                sf[i*prime[j]]=sf[i];                break;            }        }    }    //for(i=1;i<=20;++i)printf("%d:%I64d\n",i,sf[i]);    for(i=1;i<=10000000;++i)sf[i]=(sf[i-1]+sf[i]*i)%Mod;    int N,M,ans;    scanf("%d%d",&N,&M);    if(N>M)swap(N,M);    ans=0;    for(i=1;i<=N;){        j=i;        i=min(N/(N/i),M/(M/i))+1;        ans=(ans+(sf[i-1]-sf[j-1])*(N/j)%Mod*(N/j+1)%Mod*(M/j)%Mod*(M/j+1))%Mod;    }    cout<<(ans*pow(4,Mod-2)%Mod+Mod)%Mod;}

杜教版：

#include<cstdio>#include<iostream>using namespace std;#include<algorithm>#include<cstring>const int S=50000;#define Mod 20101009#define Mod6 120606054bool p[50005];int prime[10005];typedef long long LL;LL F[2][205],inis[50005],sg[50005];LL pow(int a,int x){    LL ans=1,prod=a;    for(;x;x>>=1,prod=prod*prod%Mod)        if(x&1)            ans=ans*prod%Mod;    return ans;}LL cal(int x){    return (LL)x*(x<<1|1)%Mod6*(x+1)%Mod6/6;}void work(int N,LL F[]){    F[0]=1;    while(N/F[0]>S)++F[0];    ++F[0];    for(int i=F[0],j,x,r;--i;){        x=N/i;        F[i]=1;        for(j=2;i*j<F[0];++j)F[i]=(F[i]-j*j*F[i*j])%Mod;        for(;j<=x;j=r+1){            r=x/(x/j);            F[i]=(F[i]-(cal(r)-cal(j-1))%Mod*inis[x/j])%Mod;        }    }    //for(int i=1;i<F[0];++i)cout<<"F("<<N<<"/"<<i<<"="<<N/i<<")="<<F[i]<<endl;}int N,M;LL query(int x){    if(x==0)return 0;    //printf("----query(%d)-----\n",x);    int o,k,i,j,n;    if(N/(N/x)==x){        o=0;        k=N/x;    }    else{        o=1;        k=M/x;    }    LL ans=0;    //cout<<o<<" "<<k<<endl;    if(x>S){        for(i=1;i*k<F[o][0];i=j+1){            j=x/(x/i);            ans=(ans+((LL)(i+j)*(j-i+1)>>1)%Mod*F[o][i*k])%Mod;            //cout<<"Get("<<i<<"):"<<((LL)(i+j)*(j-i+1)>>1)%Mod*F[o][i*k]%Mod<<endl;        }        //puts("");        for(;i<=x;i=j+1){            j=x/(x/i);            ans=(ans+((LL)(i+j)*(j-i+1)>>1)%Mod*inis[x/i])%Mod;            //cout<<"Get("<<i<<"):"<<((LL)(i+j)*(j-i+1)>>1)%Mod*inis[x/i]%Mod<<endl;        }    }    else        for(i=1;i<=x;i=j+1){            j=x/(x/i);            ans=(ans+((LL)(i+j)*(j-i+1)>>1)%Mod*inis[x/i])%Mod;        }    //cout<<"ans="<<ans<<endl;    return ans;}int main(){    int i,j;    inis[1]=sg[1]=1;    for(i=2;i<=50000;++i){        if(!p[i]){            prime[++prime[0]]=i;            inis[i]=-(LL)i*i%Mod;            sg[i]=1-i;        }        for(j=1;j<=prime[0]&&i*prime[j]<=50000;++j){            p[i*prime[j]]=1;            if(i%prime[j]){                inis[i*prime[j]]=inis[i]*inis[prime[j]]%Mod;                sg[i*prime[j]]=sg[i]*(1-prime[j]);            }            else{                sg[i*prime[j]]=sg[i];                break;            }        }    }    //for(i=1;i<=20;++i)printf("%d:%I64d %I64d\n",i,inis[i],sg[i]*i);    for(i=2;i<=50000;++i){        inis[i]=(inis[i]+inis[i-1])%Mod;        sg[i]=(sg[i-1]+sg[i]*i)%Mod;    }    //for(i=1;i<=105;++i)printf("%d:%I64d\n",i,sg[i]);    scanf("%d%d",&N,&M);    if(N>M)swap(N,M);    work(N,F[0]),work(M,F[1]);    LL ans=0;    if(S<N){        for(i=1,j=min(N/(N/i),M/(M/i));j<=S;i=j+1,j=min(N/(N/i),M/(M/i))){            ans=(ans+(sg[j]-sg[i-1])%Mod*(N/i)%Mod*(N/i+1)%Mod*(M/i)%Mod*(M/i+1))%Mod;            //cout<<"Get("<<i<<")="<<(sg[j]-sg[i-1])%Mod*(N/i)%Mod*(N/i+1)%Mod*(M/i)%Mod*(M/i+1)%Mod<<endl;        }        for(;i<=N;i=j+1){            j=min(N/(N/i),M/(M/i));            ans=(ans+(query(j)-query(i-1))%Mod*(N/i)%Mod*(N/i+1)%Mod*(M/i)%Mod*(M/i+1))%Mod;            //cout<<"Get("<<i<<")="<<(query(j)-query(i-1))%Mod*(N/i)%Mod*(N/i+1)%Mod*(M/i)%Mod*(M/i+1)%Mod<<endl;        }    }    else        for(i=1;i<=N;i=j+1){            j=min(N/(N/i),M/(M/i));            ans=(ans+(sg[j]-sg[i-1])%Mod*(N/i)%Mod*(N/i+1)%Mod*(M/i)%Mod*(M/i+1))%Mod;            //cout<<"Get("<<i<<")="<<(sg[j]-sg[i-1])%Mod*(N/i)%Mod*(N/i+1)%Mod*(M/i)%Mod*(M/i+1)<<endl;        }    cout<<(ans*pow(4,Mod-2)%Mod+Mod)%Mod<<endl;}