(color{#0066ff}{ 题目描述 })
这是一道非常简单的数学题。
最近 LzyRapxLzyRapx 正在看 mathematics for computer science 这本书，在看到数论那一章的时候， LzyRapxLzyRapx 突然想到这样一个问题。
设
[F(n)=sum_{i=1}^nsum_{j=1}^ifrac{mathrm{lcm}(i,j)}{mathrm{gcd}(i,j)}
]
其中，(mathrm{lcm}(a,b)) 表示 (a) 和 (b) 的最小公倍数，(mathrm{gcd}(a,b)) 表示 (a) 和 (b) 的最大公约数。
给定 (n) ，让你求： (F(n) mod1000000007)。
(L) 同学太菜啦，QAQ，并不会做这道简单题，所以他想请你帮他解决这个问题。
(color{#0066ff}{输入格式})
输入一行，一个正整数 (n\,(1 le n le 10^9))。
(color{#0066ff}{输出格式})
输出 (F(n)) ，对 (10^9 + 7) 取模。
(color{#0066ff}{输入样例})
5
(color{#0066ff}{输出样例})
84
(color{#0066ff}{数据范围与提示})
对于所有数据，(1 le n le 10^9)。
(color{#0066ff}{ 题解 })
开始推式子
题目要求
[F(n)=sum_{i=1}^nsum_{j=1}^ifrac{mathrm{lcm}(i,j)}{mathrm{gcd}(i,j)}
]
可以考虑求
[F(n)=sum_{i=1}^nsum_{j=1}^nfrac{mathrm{lcm}(i,j)}{mathrm{gcd}(i,j)}
]
然后发现我们要求的答案是下图S1+S3，上面的式子是整个正方形
于是答案可以等于（正方形+对角线）的一半
所以问题转为了下面的式子，顺便改写为gcd形式
[F(n)=sum_{i=1}^nsum_{j=1}^nfrac{i*j}{mathrm{gcd}(i,j)^2}
]
于是，枚举gcd
[sum_{d=1}^nsum_{i=1}^nsum_{j=1}^n [gcd(i,j)==d]frac{i*j}{d^2}
]
后面把d弄出来
[sum_{d=1}^nsum_{i=1}^{lfloor frac n d
floor}sum_{j=1}^{lfloor frac n d
floor} [gcd(i,j)==1]i*j
]
利用(mu * i = e)套进去
[sum_{d=1}^n sum_{i=1}^{lfloor frac n d
floor}sum_{j=1}^{lfloor frac n d
floor} sum_{k|gcd(i,j)} mu(k)*i*j
]
即
[sum_{d=1}^n sum_{k=1}^{lfloorfrac n d
floor}sum_{k|i} sum_{k|j} mu(k)*i*j
]
改为枚举倍数
[sum_{d=1}^n sum_{k=1}^{lfloorfrac n d
floor} mu(k)*k^2sum _{i=1}^{lfloorfrac n {kd}
floor} sum _{j=1}^{lfloorfrac n {kd}
floor}i*j
]
然后pd换q的套路
[sum_{q=1}^nsum_{k|q} mu(k)* k^2 (sum _{i=1}^{lfloorfrac n {q}
floor}i)^2
]
[令f(i)=sum_{k|i}mu(k) * k^2, g(i)=i^2,t(i)=mu(i)*i^2
]
我们现在要求的(h=f*g)这样才能杜教筛
考虑将f拆开，可以发现(f=t*1)
于是(h=t*1*g=t*g*1)
然后可以神奇的发现
[(t*g)(n)=n^2*sum_{d|n} mu(d)=e
]
然后。。。
[h=t*g*1=e*1=1
]
额，(h(n)=1)
不过杜教筛的前半部分是要线性筛的，这东西咋弄啊？
可以发现，(mu(d)*d^2) 是一个积性函数*完全积性函数
所以(mu(d)*d^2)肯定是记性函数
而f是这东西卷上一个1，所以肯定是积性函数，随便筛就行了qwq
#include<bits/stdc++.h>
#define LL long long
LL in() {
char ch; LL x = 0, f = 1;
while(!isdigit(ch = getchar()))(ch == '-') && (f = -f);
for(x = ch ^ 48; isdigit(ch = getchar()); x = (x << 1) + (x << 3) + (ch ^ 48));
return x * f;
}
const int mod = 1e9 + 7;
const int maxn = 4e6 + 10;
std::map<int, LL> mp;
LL pri[maxn], mu[maxn], tot, six;
bool vis[maxn];
LL ksm(LL x, LL y) {
LL re = 1LL;
while(y) {
if(y & 1) re = re * x % mod;
x = x * x % mod;
y >>= 1;
}
return re;
}
void predoit() {
six = ksm(6, mod - 2);
mu[1] = 1;
vis[1] = true;
for(LL i = 2; i < maxn; i++) {
if(!vis[i]) pri[++tot] = i, mu[i] = 1 - (LL)i * i;
for(int j = 1; j <= tot && (LL)i * pri[j] < maxn; j++) {
vis[i * pri[j]] = true;
if(i % pri[j] == 0) {
mu[i * pri[j]] = mu[i];
break;
}
else mu[i * pri[j]] = mu[i] * mu[pri[j]];
}
}
for(int i = 2; i < maxn; i++) mu[i] = (mu[i] + mu[i - 1]) % mod;
}
LL gettot(LL x) {
x %= mod;
LL ans = (x * (x + 1) >> 1) % mod;
return ans * ans % mod;
}
LL getpfh(LL x) {
x %= mod;
LL ans = x;
ans = (ans * (x + 1)) % mod;
ans = (ans * ((x << 1LL) + 1)) % mod;
return ans * six % mod;
}
LL getmu(LL n) {
if(n < maxn) return mu[n];
if(mp.count(n)) return mp[n];
LL ans = n;
for(LL l = 2, r; l <= n; l = r + 1) {
r = n / (n / l);
LL tot1 = ((getpfh(r) - getpfh(l - 1)) % mod + mod) % mod;
tot1 = (tot1 * getmu(n / l)) % mod;
ans = ((ans - tot1) % mod + mod) % mod;
}
return mp[n] = ans;
}
LL work(LL n) {
LL ans = 0;
for(LL l = 1, r; l <= n; l = r + 1) {
r = n / (n / l);
LL tot1 = ((getmu(r) - getmu(l - 1)) % mod + mod) % mod;
tot1 = tot1 * gettot(n / l) % mod;
ans = (ans + tot1) % mod;
}
return ((ans + n) % mod) * ksm(2, mod - 2) % mod;
}
int main() {
predoit();
printf("%lld
", work(in()));
return 0;
}