Easy Math Problem

推式子

$$\begin{gathered} \sum _{i=1}^n\sum _{j=1}^{n}gcd(i,j{)}^{K}lcm(i,j)[gcd(i,j)\in prime] \\ \sum _{i=1}^n\sum _{j=1}^{n}gcd(i,j{)}^{K-1}ij[gcd(i,j)\in prime] \\ \sum _{d\in prime}^n{d}^{K+1}\sum _{i=1}^{\frac{n}{d}}\sum _{j=1}^{\frac{n}{d}}ij[gcd(i,j)==1] \\ 对\sum _{i=1}^n\sum _{j=1}^{n}ij[gcd(i,j)==1]化简 \\ 2(\sum _{i=1}^{n}i\frac{i\varphi (i)+[i==1]}{2})-1=\sum _{i=1}^n{i}^2\varphi (i) \\ \sum _{d\in prime}^n{d}^{K+1}\sum _{i=1}^{\frac{n}{d}}{i}^2\varphi (i) \end{gathered}$$

接下来就是快了的数论分块了，首先我们用min_25得到区间$[l,r]$部分的质数的贡献，然后再通过杜教筛得到后面，两者相乘然后累加即可，对于$\sum _{i=1}^n{i}^{k+1}$这一部分求和，显然我们可以通过拉格朗日插值来求解。

代码

/*Author : lifehappy
*/
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#include <bits/stdc++.h>#define mp make_pair
#define pb push_back
#define endl '\n'
#define mid (l + r >> 1)
#define lson rt << 1, l, mid
#define rson rt << 1 | 1, mid + 1, r
#define ls rt << 1
#define rs rt << 1 | 1using namespace std;typedef long long ll;
typedef unsigned long long ull;
typedef pair<int, int> pii;const double pi = acos(-1.0);
const double eps = 1e-7;
const int inf = 0x3f3f3f3f;inline ll read() {ll f = 1, x = 0;char c = getchar();while(c < '0' || c > '9') {if(c == '-')    f = -1;c = getchar();}while(c >= '0' && c <= '9') {x = (x << 1) + (x << 3) + (c ^ 48);c = getchar();}return f * x;
}const int N = 1e6 + 10, mod = 1e9 + 7, inv2 = (mod + 1) >> 1, inv6 = (mod + 1) / 6;namespace Djs {int prime[N], cnt;ll phi[N];bool st[N];void init() {phi[1] = 1;for(int i = 2; i < N; i++) {if(!st[i]) {prime[++cnt] = i;phi[i] = i - 1;}for(int j = 1; j <= cnt && 1ll * i * prime[j] < N; j++) {st[i * prime[j]] = 1;if(i % prime[j] == 0) {phi[i * prime[j]] = phi[i] * prime[j];break;}phi[i * prime[j]] = phi[i] * (prime[j] - 1);}}for(int i = 1; i < N; i++) {phi[i] = (phi[i - 1] + 1ll * i * i % mod * phi[i] % mod) % mod;}}ll calc1(ll n) {n %= mod;return (n * (n + 1) % mod * inv2 % mod) * (n * (n + 1) % mod * inv2 % mod) % mod;}ll calc2(ll n) {n %= mod;return n * (n + 1) % mod * (2 * n + 1) % mod * inv6 % mod;}unordered_map<ll, ll> ans_s;ll S(ll n) {if(n < N) return phi[n];if(ans_s.count(n)) return ans_s[n];ll ans = calc1(n);for(ll l = 2, r; l <= n; l = r + 1) {r = n / (n / l);ans = ((ans - (calc2(r) - calc2(l - 1)) * S(n / l) % mod) % mod + mod) % mod;}return ans_s[n] = ans;}
}namespace Lagrange {const int N = 110;ll fac[N], pre[N], suc[N], inv[N], prime[N], sum[N], mu[N], n, k, cnt;bool st[N];ll quick_pow(ll a, int n) {ll ans = 1;while(n) {if(n & 1) ans = ans * a % mod;a = a * a % mod;n >>= 1;}return ans;}void init1() {fac[0] = inv[0] = 1;for(int i = 1; i < N; i++) {fac[i] = 1ll * fac[i - 1] * i % mod;}inv[N - 1] = quick_pow(fac[N - 1], mod - 2);for(int i = N - 2; i >= 1; i--) {inv[i] = 1ll * inv[i + 1] * (i + 1) % mod;}}void init() {for(int i = 1; i < N; i++) {st[i] = 0;}cnt = 0;sum[1] = 1;for(int i = 2; i < N; i++) {if(!st[i]) {prime[cnt++] = i;sum[i] = quick_pow(i, k);}for(int j = 0; j < cnt && i * prime[j] < N; j++) {st[i * prime[j]] = 1;sum[i * prime[j]] = 1ll * sum[i] * sum[prime[j]] % mod;if(i % prime[j] == 0) break;}}for(int i = 1; i < N; i++) {sum[i] = (sum[i] + sum[i - 1]) % mod;}}ll solve(ll n) {n %= mod;ll ans = 0;pre[0] = suc[k + 3] = 1;for(int i = 1; i <= k + 2; i++) pre[i] = 1ll * pre[i - 1] * (n - i) % mod;for(int i = k + 2; i >= 1; i--) suc[i] = 1ll * suc[i + 1] * (n - i) % mod;for(int i = 1; i <= k + 2; i++) {ll a = 1ll * pre[i - 1] * suc[i + 1] % mod, b = 1ll * inv[i - 1] * inv[k + 2 - i] % mod;if((k + 2 - i) & 1) b *= -1;ans = ((ans + 1ll * sum[i] * a % mod * b % mod) % mod + mod) % mod;}return (ans - 1 + mod) % mod;}
}namespace Min_25 {int prime[N], id1[N], id2[N], cnt, m, k, T;ll g[N], sum[N], a[N], n;bool st[N];int ID(ll x) {return x <= T ? id1[x] : id2[n / x];}void init(ll x, int y) {n = x, k = y;cnt = 0, m = 0;T = sqrt(n + 0.5);for(int i = 2; i <= T; i++) {if(!st[i]) {prime[++cnt] = i;sum[cnt] = (sum[cnt - 1] + Lagrange::quick_pow(i, k)) % mod;}for(int j = 1; j <= cnt && 1ll * i * prime[j] <= T; j++) {st[i * prime[j]] = 1;if(i % prime[j] == 0) {break;}}}for(ll l = 1, r; l <= n; l = r + 1) {r = n / (n / l);a[++m] = n / l;if(a[m] <= T) id1[a[m]] = m;else id2[n / a[m]] = m;g[m] = Lagrange::solve(a[m]);}for(int j = 1; j <= cnt; j++) {for(int i = 1; i <= m && 1ll * prime[j] * prime[j] <= a[i]; i++) {g[i] = ((g[i] - Lagrange::quick_pow(prime[j], k) * (g[ID(a[i] / prime[j])] - sum[j - 1]) % mod) % mod + mod) % mod;}}for(int i = 1; i <= T; i++) {st[i] = 0;}}ll solve(ll x) {if(x <= 1) return 0;return g[ID(x)];}
}int main() {// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);Djs::init();Lagrange::init1();int T = read();while(T--) {ll n = read(), k = read() + 1;Lagrange::k = k;Lagrange::init();Min_25::init(n, k);ll ans = 0;for(ll l = 1, r; l <= n; l = r + 1) {r = n / (n / l);ans = (ans + (Min_25::solve(r) - Min_25::solve(l - 1)) * Djs::S(n / l) % mod) % mod;}printf("%lld\n", (ans % mod + mod) % mod);}return 0;
}