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 second

typedef 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