1847 奇怪的数学题

推式子

$$\begin{gathered} \sum _{i=1}^n\sum _{j=1}^{n}sgcd(i,j{)}^k \\ \sum _{d=1}^{n}sgcd(d{)}^k\sum _{i=1}^{\frac{n}{d}}\sum _{j=1}^{\frac{n}{d}}[gcd(i,j)=1] \\ \sum _{d=1}^{n}sgcd(d{)}^k(2\sum _{i=1}^{\frac{n}{d}}\varphi (i)-1) \end{gathered}$$

接下来我们只要求得$sgcd(d{)}^k$的前缀和就可以快乐地数论分块了，

考虑到这个是求最小质因子，Min_25筛在筛g函数的时候有一步是枚举最小质因子，

$sgcd=\frac{sgcd}{min(p)}$，对应的就是这一步，所以我们只要舍弃最小质因子，然后统计贡献就行了，

当然最后我们要加上质数的贡献，也就是质数的个数。

考虑到要求自然幂数前n项和，这题模数是个头疼的地方，无法用拉格朗日插值，只能乖乖的写第二类斯特林数了。

代码

/*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;
}#define uint unsigned intconst int N = 1e6 + 10;uint quick_pow(uint a, int n) {uint ans = 1;while(n) {if(n & 1) ans = ans * a;a = a * a;n >>= 1;}return ans;
}namespace Sum {uint S[60][60], k;uint Sumk(int n) {uint ret = 0;for(int i = 1; i <= k; ++i) {uint s = 1;for(int j = 0; j <= i; ++j)if((n + 1 - j) % (i + 1)) s *= n + 1 - j;else s *= (n + 1 - j) / (i + 1);ret += S[k][i] *s;}return ret;}void init() {S[0][0] = 1;for(int i = 1; i <= k; i++)for(int j = 1; j <= i; j++)S[i][j] = S[i-1][j-1] + j * S[i-1][j];}}namespace Min_25 {uint prime[N], id1[N], id2[N], m, cnt, k, T;uint a[N], g1[N], sum1[N], g2[N], sum2[N], f[N], ans[N], n;bool st[N];int ID(ll x) {return x <= T ? id1[x] : id2[n / x];}void init() {T = sqrt(n + 0.5);for(int i = 2; i <= T; i++) {if(!st[i]) {prime[++cnt] = i;sum1[cnt] = sum1[cnt - 1] + 1;f[cnt] = quick_pow(i, k);sum2[cnt] = sum2[cnt - 1] + f[cnt];}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;g1[m] = a[m] - 1;g2[m] = Sum::Sumk(a[m]) - 1;}for(int j = 1; j <= cnt; j++) {for(int i = 1; i <= m && 1ll * prime[j] * prime[j] <= a[i]; i++) {g1[i] -= (g1[ID(a[i] / prime[j])] - sum1[j - 1]);g2[i] -= f[j] * (g2[ID(a[i] / prime[j])] - sum2[j - 1]);ans[i] += g2[ID(a[i] / prime[j])] - sum2[j - 1];}}for(int i = 1; i <= m; i++) {ans[i] += g1[i];}}uint solve(int x) {if(x <= 1) return 0;return ans[ID(x)];}
}namespace Djs {uint prime[N], phi[N], cnt;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];}}unordered_map<int, uint> ans_s;uint S(int n) {if(n < N) return phi[n];if(ans_s.count(n)) return ans_s[n];uint ans = 1ll * n * (n + 1) / 2;for(uint l = 2, r; l <= n; l = r + 1) {r = n / (n / l);ans -=(r - l + 1) * S(n / l);}return ans_s[n] = ans;}
}int main() {// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);int n = read(), k = read();Sum::k = k;Sum::init();Djs::init();Min_25::n = n, Min_25::k = k;Min_25::init();uint ans = 0;for(uint l = 1, r; l <= n; l = r + 1) {r = n / (n / l);ans += (Min_25::solve(r) - Min_25::solve(l - 1)) * (2 * Djs::S(n / l) - 1);}cout << ans << endl;return 0;
}