UVA11426 GCD - Extreme (II) (歐拉函數/莫比烏斯反演)

UVA11426 GCD - Extreme (II)

題目描述

PDF

輸入輸出格式

輸入格式：

輸出格式：

輸入輸出樣例

輸入樣例#1：

10
100
200000
0

輸出樣例#1：

67
13015
143295493160

Solution

這道題我用莫比烏斯反演和歐拉函數都寫了一遍,發現歐拉函數比莫比烏斯反演優秀?

求全部\(gcd=k\)的數對的個數,記做\(f[k],ans=\sum_{i=1}^{n}(f[i]-1)\),爲何還要-1,咱們注意到\(j=i+1\),本身與本身是不算的,再乘上這個數的大小\(k\)就能夠了

咱們發現要求\(\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}{gcd(i,j)}\),咱們令\(gcd(i,j)=k\),則必有\(gcd(a\times k,b\times k)=k\to gcd(a,b)=1\),咱們枚舉這兩個數中大的那個,另外一個數就有\(phi[i](1<=i <=n/k)\)個,因此\(f[n]=\sum_{i=1}^{n/k}\phi(i)\)

篩一下歐拉函數求前綴和就能夠了

那麼?莫比烏斯反演怎麼寫呢?
\(ans=\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}gcd(i,j)\)
\(ans=\sum_{i=1}^{n}\sum_{j=1}^{n}gcd(i,j)\)

我麼觀察這兩個式子的區別,一個是從j從1開始,另外一個是從i+1開始,這兩個式子是具備幾何意義的,上面的式子求得是一個三角形的答案(且不包括對角線),下面的是矩形.那麼咱們用下面的式子-對角線上的答案再除以2就是上面的答案了

(對角線上的答案就是\(gcd(i,i)=i\),因此一個等比數列求和就行了)

下面的式子很好反演,套路套路....
\[\sum_{d=1}^{n}d\sum_{i=1}^{n}\sum_{j=1}^{n}[gcd(i,j)==d]\]
\[\sum_{d=1}^{n}d\sum_{p=1,p|i,p|j}^{\lfloor\frac{n}{d}\rfloor}\mu(p) \lfloor\frac{n}{d\times p}\rfloor\lfloor\frac{n}{d\times p}\rfloor\]

線性篩莫比烏斯函數,而後套個整除分塊就行了(仍是沒有第一種方法快,若是莫比烏斯反演有比博主還快的方法,請告知我)

Code1

#include<bits/stdc++.h>
#define lol long long
#define il inline
#define rg register
#define Min(a,b) (a)<(b)?(a):(b)
#define Max(a,b) (a)>(b)?(a):(b)
#define NN 4000000

using namespace std;

const int N=4e6+10;
int n,tot;
lol phi[N],prime[N];
bool vis[N];

il void init() {
    phi[1]=1;
    for(rg int i=2;i<=NN;i++) {
        if(!vis[i]) prime[++tot]=i,phi[i]=i-1;
        for(rg int j=1;j<=tot && i*prime[j]<=NN;j++) {
            vis[i*prime[j]]=1;
            if(i%prime[j]==0) {
                phi[i*prime[j]]=phi[i]*prime[j];
            }
            else phi[i*prime[j]]=phi[i]*(prime[j]-1);
        }
    }
    for(rg int i=1;i<=NN;i++) phi[i]+=phi[i-1];
}
int main()
{
    ios::sync_with_stdio(0);
    init();
    while(cin>>n) {
        lol ans=0;
        if(n==0) break;
        for(rg int i=1;i<=n;i++) ans+=1ll*(phi[n/i]-1)*i;
        cout<<ans<<endl;
    }
}

Code2

#include<bits/stdc++.h>
#define in(i) (i=read())
#define il extern inline
#define rg register
#define mid ((l+r)>>1)
#define ll(x) (x<<1)
#define rr(x) (x<<1|1)
#define Min(a,b) ((a)<(b)?(a):(b))
#define Max(a,b) ((a)>(b)?(a):(b))
#define lol __int128
using namespace std;

const lol N=4e6+10;

lol read() {
    lol ans=0, f=1; char i=getchar();
    while (i<'0' || i>'9') {if(i=='-') f=-1; i=getchar();}
    while (i>='0' && i<='9') ans=(ans<<1)+(ans<<3)+(i^48), i=getchar();
    return ans*f;
}

lol cnt;
lol vis[N]={0,1},prime[N],mu[N]={0,1};

void init() {
    for (lol i=2;i<=N-10;i++) {
        if(!vis[i]) prime[++cnt]=i,mu[i]=-1;
        for (lol j=1;j<=cnt && prime[j]*i<=N-10;j++) {
            vis[i*prime[j]]=1;
            if(i%prime[j]==0) break;
            mu[i*prime[j]]=-mu[i];
        }
    }for (lol i=1;i<=N-10;i++) mu[i]+=mu[i-1];
}

lol work(lol d,lol n,lol ans=0) {
    n/=d;
    for (lol l=1,r;l<=n;l=r+1) {
        r=n/(n/l);
        ans+=(mu[r]-mu[l-1])*(n/l)*(n/l);
    }return ans;
}

int main()
{
    long long ans,n; init();
    while(scanf("%lld",&n)==1 && n) {
        ans=0;
        for (lol i=1;i<=n;i++) ans+=i*work(i,n);
        printf("%lld\n",(ans-(n+1)*n/2)/2);
    }
    return 0;
}