$$\begin{gathered} \sum _{a_1=1}^m\sum _{a_2=1}^m\sum _{a_3=1}^m\cdots \sum _{a_{n-1}}^m\sum _{a_n}^m[gcd(a_1,a_2,a_3,\dots ,a_{n-1},a_n)==d](a_1,a_2,a_3,\dots ,a_{n-1},a_n{)}^k \\ =d^{kd}\sum _{a_1=1}^{\frac{m}{d}}\sum _{a_2=1}^{\frac{m}{d}}\sum _{a_3=1}^{\frac{m}{d}}\cdots \sum _{a_{n-1}}^{\frac{m}{d}}\sum _{a_n}^{\frac{m}{d}}[gcd(a_1,a_2,a_3,\dots ,a_{n-1},a_n)==1](a_1,a_2,a_3,\dots ,a_{n-1},a_n{)}^k \\ =d^{kd}\sum _{i=1}^{\frac{m}{d}}{i}^{kd}\mu (i)\sum _{a_1=1}^{\frac{m}{id}}\sum _{a_2=1}^{\frac{m}{id}}\sum _{a_3=1}^{\frac{m}{id}}\cdots \sum _{a_{n-1}=1}^{\frac{m}{id}}\sum _{a_n=1}^{\frac{m}{id}}(\prod _{j=1}^n{a}_i{)}^k \\ 这是一个多项式 \\ =d^{kd}\sum _{i=1}^{\frac{m}{d}}{i}^{kd}\mu (i)(\sum _{j=1}^{\frac{m}{id}}{i}^K{)}^n \end{gathered}$$

到这里这道题目就化简完成了，只需要通过简单的数列求和加欧拉降幂即可得到答案。

/*Author : lifehappy
*/
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#include <bits/stdc++.h>using namespace std;typedef long long ll;const int inf = 0x3f3f3f3f;
const double eps = 1e-6;const int N = 1e5 + 10, mod = 59964251, phi = 59870352;int mu[N], prime[N], cnt;ll m, d, k, n, sum[N];bool st[N];char str[N];ll quick_pow(ll a, int n, int mod = 59964251) {ll ans = 1;while(n) {if(n & 1) ans = ans * a % mod;a = a * a % mod;n >>= 1;}return ans;
}void init() {memset(st, 0, sizeof st);cnt = 0;mu[1] = 1, sum[1] = 1;for(int i = 2; i < N; i++) {if(!st[i]) {mu[i] = -1;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]] = sum[i] * sum[prime[j]] % mod;if(i % prime[j] == 0) break;mu[i * prime[j]] = -mu[i];}}   for(int i = 1; i < N; i++) {sum[i] = (sum[i] + sum[i - 1]) % mod;}
}ll solve(ll m) {ll ans = 0;for(ll i = 1; i <= m; i++) {ans = ans + 1ll * mu[i] * quick_pow(i, k * n % phi + phi) % mod * quick_pow(sum[m / i], n + phi) % mod;}return (ans % mod + mod) % mod;
}int main() {// freopen("in.txt", "r", stdin);// freopen("out.txt", "w", stdout);// ios::sync_with_stdio(false), cin.tie(0), cout.tie(0);int T;scanf("%d", &T);while(T--) {scanf("%s %lld %lld %lld\n", str + 1, &m, &d, &k);init();int len = strlen(str + 1);n = 0;for(int i = 1; i <= len; i++) {n = n * 10 + str[i] - '0';n %= phi;}ll ans = quick_pow(d, k * n % phi + phi) * solve(m / d) % mod;printf("%lld\n", ans);}return 0;
}