HDU 5341 Gcd and Lcm

题意：给你一个N，让你求∑​i=1​n​​∑​j=1​n​​∑​k=1​n​​∑​l=1​n​​lcm(gcd(i,j),gcd(k,l))

参考论文：贾志鹏线性筛

首先，考虑子问题，求这样一个东西：∑​i=1​n​​∑​j=1​n​​[(i,j)=d]

（中括号是一个函数：若其中的命题为真，返回1，否则返回0）

转化一下，就是求1到n/d两两之间互素的数然后再乘以d就是gcd为d的数了。

即：∑​i=1​n​​∑​j=1​n​​[(i,j)=d] = ∑、/123​i=1​n/d​​∑​j=1​n/d​​[(i,j)=1]

那么如题意所说，窝们可以定义一个函数 f(x)=∑、/123​i=1​x​​∑​j=1​x​​[(i,j)=1]

那么对于固定的n，窝们求gcd为x的个数为f(n/x)。

那么问题就变成了

ans=∑、/123​i=1​n​​∑​j=1​n​​lcm(i,j)*f(n/i)*f(n/j)

令p=gcd(i,j)，则ans=∑​p=1​n​​∑​i=1​n/p​​∑​j=1​n/p​​p∗i∗j∗f(n/(i*p))∗f(n/(j*p))[(i,j)=1]

ans=∑​p=1​n​​∑​i=1​n/p​​∑​j=1​n/p​​p∗i∗j∗f(n/(i*p))∗f(n/(j*p))∑​d|gcd(i,j)​​​u(d)

=∑​p=1​n​​∑​i=1​n/p​​∑​j=1​n/p​​p∗i∗j∗f(n/(i*p))∗f(n/(j*p))∑​d|i&&d|j​​​u(d)

其中u(d)表示莫比乌斯函数，然后窝们想方设法将u(d)提到前面来：

=∑​p=1​n​​  p*∑​d​​​u(d)∑​i=1​n/p/d​​   (i*d)*f(n/(i*d*p))∑​j=1​n/p/d​​  (j*d)*f(n/(j*d*p))

可以将ij视为相同的两个东西，同时把d提出来放到前面：

=∑​p=1​n​​  ∑​d​​​u(d)*p*d*d*(∑​i=1​n/p/d​​   i*f(n/(i*d*p)))​2 

然后，利用类似于参考论文中贾志鹏处理的方法，令t=p*d

=∑、/123​t=1​n​​(∑​i=1​n/t​​   i*f(n/(i*t)))​2* t*∑​d|t​​​u(d)*d

h(t)=t*∑​d|t​​​u(d)*d 是一个积性函数，窝们可以用线性筛的思想去O(n)的预处理出所有的h。

g(x)=∑​i=1​x​​  i*f(x/i) 并不是一个积性函数但是x/i一共有大概sqrt(x)种取值。

在这里窝们可以用一种方法处理复杂度为sqrt(x)，详细可以看代码。

那么ans=∑、/123​t=1​n​​g(n/t)​2*h(t)

同理，n/t也是只有sqrt(n)种取值，所以复杂度最大也不超过sqrt(n)*sqrt(n)<=n,最后的复杂度不超过O(n)。

#include<stdio.h>#include<cstring>#include<cmath>#include<string>#include<algorithm>#include<iostream>#include<vector>#include<queue>#include<map>#include<set>#include<stack>#include<bitset>#include<time.h>#define Msn(x,y) (memset((x),0,sizeof((x[0]))*(y+1)))#define msn(x) (memset((x),0,sizeof((x))))#define msx(x,y) (memset((x),0x7f,sizeof((x[0]))*(y+3)))#define fuck(x) cerr << #x << " <- " << x << endl#define acer cout<<"sb"<<endl#define mkp(x,y) (make_pair(x,y))#define X first#define Y secondtypedef long long ll;typedef unsigned int ui;using namespace std;const int maxn=1e7+7;int mu[maxn];int phi[maxn];bool not_pr[maxn];int pr[maxn];int pr_num;ui g[maxn],f[maxn],sum[maxn],h[maxn],sum2[maxn];void get_mobius_and_euler(int sz){    pr_num=0;    mu[1]=f[1]=sum[1]=h[1]=sum2[1]=1;    for(int i=2;i<=sz;i++)    {        if(!not_pr[i])        {            pr[pr_num]=i;            mu[i]=-1;            phi[i]=i-1;            h[i]=1-i;            pr_num++;        }        for(int j=0;pr[j]*i<=sz&&j<pr_num;j++)        {            not_pr[pr[j]*i]=1;            if(i%pr[j]==0)            {                mu[pr[j]*i]=0;                phi[i*pr[j]]=phi[i]*pr[j];                h[i*pr[j]]=h[i];                break;            }            mu[pr[j]*i]=-mu[i];            phi[i*pr[j]]=phi[i]*(pr[j]-1);            h[i*pr[j]]=h[i]*h[pr[j]];        }        f[i]=f[i-1]+phi[i]*2;        sum[i]=sum[i-1]+i;        sum2[i]=sum2[i-1]+h[i]*i;    }}int n;bool vis[maxn];ui gg(int x){    if(vis[x])return g[x];    vis[x]=1;    g[x]=0;    int nxt;    for(int i=1;i<=x;i=nxt+1)    {        nxt=x/(x/i);        g[x]+=(sum[nxt]-sum[i-1])*f[x/i];    }    return g[x];}ui work(){    ui ans=0;    int nxt;    for(int p=1;p<=n;p=nxt+1)    {        nxt=n/(n/p);        ans+=(sum2[nxt]-sum2[p-1])*gg(n/p)*gg(n/p);    }    return ans;}int main(){    get_mobius_and_euler(maxn-3);    int T;    scanf("%d",&T);    while(T--)scanf("%d",&n),cout<<work()<<endl;    return 0;}

​ssfasdfasdfasdfasdfasdf