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

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

$$\sum_{i=1}^n\sum_{j=1}^mlcm(i,j)$$ $$=\sum_{i=1}^n\sum_{j=1}^m{ij\over (i,j)}$$ $$=\sum_{g=1}^n{1\over g}\sum_{i=1}^{\lfloor {n \over g} \rfloor}\sum_{i=1}^{\lfloor {m \over g} \rfloor }ij\sum_d\mu (d)[d|i][d|j]$$ (枚举gcd) $$=\sum_{g=1}^n{1 \over g}\sum_{d=1}^{\lfloor {n \over g} \rfloor}\mu (d)(gd)^2\lfloor {n \over dg}\rfloor \lfloor {m \over dg} \rfloor$$ $$=\sum_{i=1}^n\lfloor {n \over i} \rfloor \lfloor {m\over i}\rfloor i^2\sum_{d|i}\mu (d){i\over d}$$ (将g,d卷积) $$=\sum_{i=1}^n\lfloor {n \over i} \rfloor \lfloor {m\over i}\rfloor i\sum_{d|i}\mu (d)d$$
注意到 $f(n)=n\sum_{d|n}\mu (d)d$显然是一个积性函数，考虑n= $\prod_{i=1}^kp_i^{a_i}$，则 $f(n)=n\prod_{i=1}^k(1-p_i)$.所以就可以很容易的线筛出来了。然后记一下f(n)的前缀和，便可以分块统计了。所以对于bzoj2693（bzoj2154的多组询问版）便可以做到 $O(n+T\sqrt n)$
但是对于单组询问，既然后面的那个显然是积性函数，而且我们只要它的前缀和，那它是否应该可以杜教筛出来呢？
我一开始是想对这玩意儿直接做积性前缀和和普通前缀和。。然后我发现完全做不了。但是正如tangjz给的那个例题一样，应该先将式子化简。
$$\sum_{i=1}^ni\sum_{d|i}\mu (d)d$$ $$=\sum_{i=1}^n \sum_{d|i}\mu (d)d^2{i\over d}$$ $$=\sum_{i=1}^ni\sum_{j=1}^{\lfloor {n\over i}\rfloor }\mu (d)d^2$$ 所以问题就变成了求 $f(n)=\mu (n)n^2$的前缀和，这时我们再用杜教筛，将它与f(n)=n^2卷积。
$$\sum_{i=1}^n \sum_{d|i}\mu (d)d^2({i \over d})^2$$ $$=\sum_{i=1}^ni^2\sum_{d|i}\mu (d)$$ $$=1$$ $$=\sum_{i=1}^ni^2\sum_{j=1}^{\lfloor{n\over i}\rfloor}\mu (d)d^2$$
这样就可以杜教筛了！
但是还有一个问题，我们求出了 $\sum_{i=1}^n \mu (i)i^2$，可我们要求 $\sum_{i=1}^ni\sum_{j=1}^{\lfloor {n \over i}\rfloor}\mu (j)j^2$，所以我们还需要预处理出 $n^{2\over 3}$的这玩意儿，然后就可以在 $O(n^{2\over 3})$的时间复杂度搞出来了。
普通版：

#include<cstdio>
#include<iostream>
using namespace std;
#define Mod 20101009
typedef 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 120606054
bool 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;
}