51nod 1575 Gcd and Lcm

先吐槽一下这题题解的不负责，完全就是什么都没讲

然后去膜sxt大爷

讨论区也给出了关于$\sum_{j=1}^i\sum_{k=1}^i[(i,j),(i,k)]$这个东西对于一个质数和质数的若干次幂的公式，然后就可以用洲阁筛了

以下是个人看了论文后的一些理解和总结（主要写给自己看）

部分摘自论文，侵删

以$\sqrt n$为界，因为$1 \backsim n$的数最多只有1个$>\sqrt n$的质因子，所以，设$F(x)$为$x$的函数值，$G(p)$为质数$p$的函数值

那么，有

$$\sum_{x=1}^nF(x)=\sum_{x<=n，x无>\sqrt n的质因子}^n*(1+\sum_{\sqrt n < p <=\lfloor \frac{n}{x} \rfloor p是质数}G(p))$$

以$\lfloor \frac{n}{x} \rfloor$将$x$和$p$分成$(O(\sqrt n))$组，考虑将他们筛出来

设$g[i][j]$表示$1 \backsim j$和前$i$个质数互质的数的$G$值之和，$g[\sqrt n内质数个数][p]$就是要用到的$\sum G(p)$
可得递推式$g[i][j]=g[i-1][j]-G(primei)*g[i-1][\lfloor \frac{j}{primei}\rfloor]$
由$\lfloor \dfrac{ \lfloor \frac{n}{a}\rfloor}{b} \rfloor = \lfloor \dfrac{n}{ab}\rfloor$，可得第二维只有$O(\sqrt n)$个状态

然后根据$g$的$j<primei^2$时的一些性质，优化了转移的复杂度

求出了$g$之后，再设$f[i][j]$表示$\lfloor\dfrac{n}{x}\rfloor=j$，$x$只含前$i$个质数为质因子的$F(x)$之和，同理第二维只有$O(\sqrt n)$种状态

然后与$g$类似的性质，也可以优化掉$j<primei^2$时状态，也是用一个数组存$j$上一次转移时的值，用到$j$的时候也可以$O(1)$从上一次转移时的值求出需要的$f[i][j]$，同时在某次更新到一个不会考虑的状态或者一个状态即将不再被考虑，要更新答案

$g$里面的1要注意去掉，$f$和 g <script type="math/tex" id="MathJax-Element-556">g</script>用滚动数组存储

感谢Drin_E在我学习这个东西的时候给予我的一些帮助

具体做法请参考论文和代码…
code：

#include<set>
#include<map>
#include<deque>
#include<queue>
#include<stack>
#include<cmath>
#include<ctime>
#include<bitset>
#include<string>
#include<vector>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<climits>
#include<complex>
#include<iostream>
#include<algorithm>
#define ll long long
#define un unsigned int
using namespace std;

inline void swap(int &x,int &y){x^=y;y^=x;x^=y;}
const int maxn = 2e5+1000;

un pri,p0[maxn],p[maxn],p2[maxn],pre[maxn][3],ps[maxn];
bool v[maxn];
un w[maxn],h[maxn],h1[maxn],num;
un g[2][maxn][3],la[maxn][3],tg[maxn][3],G[maxn];
un f[2][maxn];
un ret;
int n,m;

un count(int N,int t)
{
    if(!t) return N;
    un t1=N,t2=N+1,t3=2*N+1;
    if(!(t1&1)) t1>>=1; else t2>>=1;
    if(t==1) return t1*t2;
    if(t1%3==0) t1/=3;
    else if(t2%3==0) t2/=3;
    else t3/=3;
    return t1*t2*t3;
}

int main()
{
    int t; scanf("%d",&t);
    while(t--)
    {

        scanf("%d",&n); m=sqrt(n);

        memset(v,false,sizeof v); pri=0;
        for(un i=2;i<=m;i++)
        {
            if(!v[i])
            {
                p[++pri]=i;
                p2[pri]=i*i;
                pre[pri][0]=pre[pri-1][0]+1;
                pre[pri][1]=pre[pri-1][1]+i;
                pre[pri][2]=pre[pri-1][2]+p2[pri];
                ps[pri]=ps[pri-1]+3*(p2[pri]-i)+1;
            }
            for(int j=1;j<=pri;j++)
            {
                un k=p[j]*i;
                if(k>m) break;
                v[k]=true;
                if(i%p[j]==0) break;
            }
            p0[i]=pri;
        }

        num=0;
        for(int i=1;i<=n;)
        {
            int tmp=n/i,r=n/tmp;
            w[++num]=tmp;
            (tmp<=m)?h[tmp]=num:h1[r]=num;
            i=r+1;
        }

        int i1=1,i2=0,N=num;
        for(int i=1;i<=N;i++) 
        {
            v[i]=false;
            g[i1][i][0]=count(w[i],0);
            g[i1][i][1]=count(w[i],1);
            g[i1][i][2]=count(w[i],2);
            g[i2][i][0]=g[i2][i][1]=g[i2][i][2]=0;
        }
        for(int i=1;i<=pri;i++)
        {
            swap(i1,i2);
            for(;w[N]<p2[i];N--)
            {
                v[N]=true;
                la[N][0]=g[i2][N][0]+pre[i-1][0];
                la[N][1]=g[i2][N][1]+pre[i-1][1];
                la[N][2]=g[i2][N][2]+pre[i-1][2];
                tg[N][0]=la[N][0]-pre[pri][0];
                tg[N][1]=la[N][1]-pre[pri][1];
                tg[N][2]=la[N][2]-pre[pri][2];
            }
            for(int j=1;j<=N;j++)
            {
                int k=w[j]/p[i];
                (k<=m)?k=h[k]:k=h1[n/k];
                if(v[k])
                {
                    g[i2][k][0]=la[k][0]-pre[i-1][0];
                    g[i2][k][1]=la[k][1]-pre[i-1][1];
                    g[i2][k][2]=la[k][2]-pre[i-1][2];
                }
                g[i1][j][0]=g[i2][j][0]-g[i2][k][0];
                g[i1][j][1]=g[i2][j][1]-g[i2][k][1]*p[i];
                g[i1][j][2]=g[i2][j][2]-g[i2][k][2]*p2[i];
            }
        }
        for(int i=1;i<=N;i++) tg[i][0]=g[i1][i][0],tg[i][1]=g[i1][i][1],tg[i][2]=g[i1][i][2];
        for(int i=1;i<=num;i++)
        {
            if(w[i]>m) G[i]=3*(tg[i][2]-tg[i][1])+tg[i][0]-1;
            else G[i]=0;
        }

        ret=0;
        i1=1,i2=0;
        N=num; for(int i=1;i<=N;i++) v[i]=false,f[i1][i]=0;
        f[i1][1]=1;
        for(int i=1;i<=pri;i++)
        {
            swap(i1,i2);
            for(;w[N]<p2[i];N--)
            {
                v[N]=true;
                ret+=f[i2][N]*(1+G[N]);
                int u=w[N]<=m?p0[w[N]]:pri;
                if(u>i-1) ret+=f[i2][N]*(ps[u]-ps[i-1]);
            }
            for(int j=1;j<=N;j++) f[i1][j]=0;
            for(int j=1;j<=N;j++)
            {
                un s=w[j];
                un c=0,x1,x2,T=1;
                while(s)
                {
                    int k=(s<=m)?h[s]:h1[n/s];
                    if(v[k])
                    {
                        ret+=f[i2][j]*T*(1+G[k]);
                        int u=s>m?pri:p0[s];
                        if(u>i) ret+=f[i2][j]*T*(ps[u]-ps[i]);
                    }
                    else f[i1][k]+=f[i2][j]*T;
                    s/=p[i];
                    c++;
                    if(c==1) x1=1,x2=p[i]; else x1*=p[i],x2*=p2[i];
                    T=(2*c+1)*(x2*p[i]-x2)+x1;
                }
            }
        }
        for(int i=1;i<=N;i++) ret+=f[i1][i]*(1+G[i]);

        printf("%u\n",ret);
    }

    return 0;
}