題目鏈接

題意

其實就是求

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^ngcd(i,j)\]

思路

建議先看一下此題的一個弱化版

推一下式子

\[\sum\limits_{i=1}^n\sum\limits_{j=1}^ngcd(i,j)\]
\[= \sum\limits_{k=1}^nk\sum\limits_{i=1}^n\sum\limits_{j=1}^n[gcd(i,j)=k]\]

\[=\sum\limits_{k=1}^nk\sum\limits_{i=1}^{\frac{n}{k}}\sum\limits_{j=1}^{\frac{n}{k}}[gcd(i,j)=1]\]

\[=\sum\limits_{k=1}^nk\sum\limits_{i=1}^{\frac{n}{k}}2\varphi(i)-1\]

設\(\phi(i)=\varphi(1)+\varphi(2)+...+\varphi(i)\)

則原式
\[=\sum\limits_{i=1}^ni(2\phi(\frac{n}{i})-1)\]

然後就可以數論分塊啦。

至於怎麽比較快的求\(\phi(i)\)，基本的杜教篩嘍。。

代碼

//loj6074
/*
* @Author: wxyww
* @Date:   2019-03-30 12:43:48
* @Last Modified time: 2019-03-30 19:43:10
*/
#include<cstdio>
#include<iostream>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<queue>
#include<map>
#include<vector>
#include<ctime>
using namespace std;
typedef long long ll;
const int mod = 1e9 + 7,N = 1000000 + 100,inv2 = (mod + 1) / 2;

ll read() {
    ll x=0,f=1;char c=getchar();
    while(c<'0'||c>'9') {
        if(c=='-') f=-1;
        c=getchar();
    }
    while(c>='0'&&c<='9') {
        x=x*10+c-'0';
        c=getchar();
    }
    return x*f;
}
map<ll,ll>ma;
ll n,sum[N];
int dis[N],vis[N],js;
int dls(ll x) {
    if(x <= 1000000) return sum[x];
    if(ma[x]) return ma[x];
    ll ret = 1ll * x % mod * ((x + 1) % mod) % mod * inv2 % mod;
    for(ll l = 2,r;l <= x;l = r + 1) {
        r = x / (x / l);
        ret -= 1ll * (r - l + 1) % mod * dls(x / l) % mod;
        ret = (ret + mod) % mod;
    }
    return ma[x] = ret;
}
void pre() {
    sum[1] = 1;vis[1] = 1;
    int NN = min(n,1000000ll);
    for(int i = 2;i <= NN;++i) {
        if(!vis[i]) {
            dis[++js] = i;
            sum[i] = i - 1;
        }
        for(int j = 1;j <= js && dis[j] * i <= NN;++j) {
            vis[dis[j] * i] = 1;
            if(i % dis[j] == 0) {
                sum[dis[j] * i] = sum[i] * dis[j] % mod;    break;
            }
            sum[dis[j] * i] = (dis[j] - 1) * sum[i] % mod;
        }
        sum[i] += sum[i - 1];
        sum[i] %= mod;
    }
    
}
signed main() {
    n = read();
    pre();
    ll ans = 0;
    for(ll l = 1,r;l <= n;l = r + 1) {
        r = n / (n / l);
        ans = (ans + (1ll * (r - l + 1) % mod * ((r + l) % mod) % mod * inv2 % mod) % mod * ((2ll * dls(n / l) % mod) - 1 + mod) % mod) % mod;  
    }
    cout<<ans;
    return 0;
}
 

51nod1237 最大公約數之和