题目链接
思路
首先我们将原式化简：
[egin{aligned}
&sumlimits_{l_1=1}^{n}sumlimits_{l_2=1}^{n}dotssumlimits_{l_k=1}^{n}gcd(l_1,l_2,dots,l_k)^2&\
=&sumlimits_{d=1}^{n}d^2sumlimits_{l_1=1}^{n}sumlimits_{l_2=1}^{n}dotssumlimits_{l_k=1}^{n}[gcd(l_1,l_2,dots,l_k)=d]&\
=&sumlimits_{d=1}^{n}d^2sumlimits_{l_1=1}^{frac{n}{d}}sumlimits_{l_2=1}^{frac{n}{d}}dotssumlimits_{l_k=1}^{frac{n}{d}}[gcd(l_1,l_2,dots,l_k)=1]&
end{aligned}
]
后面那一堆我们用经典反演套路进行反演得到：
[egin{aligned}
sumlimits_{d=1}^{n}d^2sumlimits_{p=1}^{lfloorfrac{n}{d}
floor}mu(p)lfloorfrac{n}{dp}
floor^k
end{aligned}
]
我们枚举(T=dp)得：
[egin{aligned}
sumlimits_{T=1}^{n}lfloorfrac{n}{k}
floor^ksumlimits_{t|T}mu(t)lfloorfrac{T}{t}
floor^2
end{aligned}
]
因此题目要求的式子最后为：
[egin{aligned}
&sumlimits_{cnt=2}^{k}sumlimits_{T=1}^{n}lfloorfrac{n}{k}
floor^{cnt}sumlimits_{t|T}mu(t)lfloorfrac{T}{t}
floor^2&\
=&sumlimits_{T=1}^{n}(sumlimits_{cnt=2}^{k}lfloorfrac{n}{k}
floor^{cnt})sumlimits_{t|T}mu(t)lfloorfrac{T}{t}
floor^2&\
=&sumlimits_{T=1}^{n}(frac{lfloorfrac{n}{T}
floor imes(lfloorfrac{n}{T}
floor^k-1)}{lfloorfrac{n}{T}
floor-1}-lfloorfrac{n}{T}
floor)sumlimits_{t|T}mu(t)lfloorfrac{T}{t}
floor^2& ext{等比数列求和}
end{aligned}\
]
比赛的时候我写到这里就不会了，主要是后面那个不知道该怎么卷积进行杜教筛，赛后看了这篇博客才懂，下面思路基本上都是参考这位大佬的。
首先我们知道(mu)为积性函数，(id^2)为积性函数，两者相乘还是积性函数，因此(f(n)=sumlimits_{d|n}mu(d)lfloorfrac{n}{d} floor^2=sumlimits_{d|n}mu(d)id(frac{n}{d})^2)。
对于这个(f)的前缀和我们可以用杜教筛进行求解，设(I*h=I*f=I*mu *id^2)，将其化简：
[egin{aligned}
&sumlimits_{t|n}I(frac{n}{t})f(t)&\
=&sumlimits_{t|n}I(frac{n}{t})sumlimits_{d|t}mu(t)id(frac{t}{d})^2&\
=&sumlimits_{d|n}id(frac{n}{d})^2sumlimits_{x|frac{n}{d}}I(frac{n}{xd})mu(frac{xd}{d})&\
=&sumlimits_{d|n}id(frac{n}{d})^2sumlimits_{x|frac{n}{d}}mu(x)&\
=&sumlimits_{d|n}id(frac{n}{d})^2[frac{n}{d}=1]&\
=&n^2&
end{aligned}
]
我当时看那篇博客时一直看不懂第二步到第三步是怎么来的，后面发现其实就是枚举(x=frac{t}{d})，然后那篇博客里面仍然用(t)来表示就导致看起来很迷。
然后进行杜教筛，设(S(n)=sumlimits_{i=1}^{n}f(i)=sumlimits_{i=1}^{n}i^2-sumlimits_{i=2}^{n}S(frac{n}{i})=frac{n(n+1)(2n+1)}{6}-sumlimits_{i=2}^{n}S(frac{n}{i}))。
在杜教筛之前我们需要先预处理出(nleq1000000)的情况，我们可以发现(f(1)=1,f(p^c)=p^{2c}-p^{2c-2},f(p_1^{c_1}p_2^{c_2}dots p_n^{c_n})=f(p_1^{c_1})f(p_2^{c_2})dots f(p_n^{c_n}))，然后就可以在线性筛时一并进行预处理掉辣~
对于最终答案里面的(sumlimits_{T=1}^{n}(frac{lfloorfrac{n}{T} floor imes(lfloorfrac{n}{T} floor^k-1)}{lfloorfrac{n}{T} floor-1}-lfloorfrac{n}{T} floor))我们可以用数论分块进行处理，注意处理公比为(1)的情况，然后这题就做完了。
代码
#include <set>
#include <map>
#include <deque>
#include <queue>
#include <stack>
#include <cmath>
#include <ctime>
#include <bitset>
#include <cstdio>
#include <string>
#include <vector>
#include <cassert>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <algorithm>
#include <unordered_map>
using namespace std;
typedef long long LL;
typedef pair<LL, LL> pLL;
typedef pair<LL, int> pLi;
typedef pair<int, LL> pil;;
typedef pair<int, int> pii;
typedef unsigned long long uLL;
#define lson rt<<1
#define rson rt<<1|1
#define lowbit(x) x&(-x)
#define name2str(name) (#name)
#define bug printf("*********
")
#define debug(x) cout<<#x"=["<<x<<"]" <<endl
#define FIN freopen("D://Code//in.txt","r",stdin)
#define IO ios::sync_with_stdio(false),cin.tie(0)
const double eps = 1e-8;
const int mod = 1000000007;
const int maxn = 1000000 + 7;
const double pi = acos(-1);
const int inf = 0x3f3f3f3f;
const LL INF = 0x3f3f3f3f3f3f3f3fLL;
bool v[maxn];
char s[100007];
int _, n, k0, k, cnt, inv;
int p[maxn], sum[maxn], mu[maxn];
unordered_map<int,int> dp;
int add(int x, int y) {
x += y;
if(x >= mod) x -= mod;
if(x < 0) x += mod;
return x;
}
void init() {
sum[1] = 1;
for(int i = 2; i < maxn; ++i) {
if(!v[i]) {
sum[i] = add(1LL * i * i % mod, -1);
p[cnt++] = i;
}
for(int j = 0; j < cnt && i * p[j] < maxn; ++j) {
v[i*p[j]] = 1;
if(i % p[j] == 0) {
sum[i*p[j]] = 1LL * sum[i] * p[j] % mod * p[j] % mod;
break;
}
sum[i*p[j]] = 1LL * sum[i] * add(1LL * p[j] * p[j] % mod, -1) % mod;
}
}
for(int i = 2; i < maxn; ++i) {
sum[i] = add(sum[i], sum[i-1]);
}
}
LL qpow(LL x, int n) {
LL res = 1;
while(n) {
if(n & 1) res = res * x % mod;
x = x * x % mod;
n >>= 1;
}
return res;
}
LL get_sum(int n, int k) {
LL ans = 1LL * n * add(qpow(n, k), -1) % mod;
int tmp = qpow(add(n, -1), mod - 2);
ans = ans * tmp % mod;
ans = add(ans, -n);
return ans;
}
int dfs(int x) {
if(x < maxn) return sum[x];
if(dp.count(x)) return dp[x];
LL ans = 1LL * x * (x + 1) % mod * (2 * x + 1) % mod * inv % mod;
for(int l = 2, r; l <= x; l = r + 1) {
r = x / (x / l);
ans = add(ans, -(1LL * (r - l + 1) * dfs(x / l) % mod));
}
return dp[x] = ans;
}
int main() {
#ifndef ONLINE_JUDGE
FIN;
#endif // ONLINE_JUDGE
init();
inv = qpow(6, mod - 2);
for(scanf("%d", &_); _; _--) {
scanf("%d%s", &n, s + 1);
int len = strlen(s + 1);
k0 = k = 0;
for(int i = 1; i <= len; ++i) {
k = (10LL * k + (s[i] - '0')) % (mod - 1);
k0 = (10LL * k0 + (s[i] - '0')) % mod;
}
k0 = (k0 - 1 + mod) % mod;
LL ans = 0;
for(int l = 1, r; l <= n; l = r + 1) {
r = n / (n / l);
int x = n / l;
LL tmp;
if(x == 1) tmp = k0;
else tmp = get_sum(x, k);
ans = add(ans, tmp * add(dfs(r), -dfs(l - 1)) % mod);
}
printf("%lld
", ans);
}
return 0;
}