需要用到概率论的容斥定理以及计算1 ^ 4 + 2 ^ 4 + ……+ n ^ 4的计算公式1^4+2^4+……+n^4=n(n+1)(2n+1)(3n^2+3n-1)/30

#pragma comment(linker,"/STACK:1024000000,1024000000")
#include <cstdio>
#include <cmath>#define LL long long
const LL mod = 1e9 + 7;
#define MAX 10000int len, prime[MAX], num[30];
bool vis[MAX + 5];
LL n, sum, pi;void get_prime(){len = 0;for(int i = 2; i<=MAX; ++i){if(!vis[i]) prime[len++] = i;for(int j = i * i; j <= MAX; j+=i) vis[j] = 1;}
}LL power(LL x, LL y){if(y == 0) return 1;if(y == 1) return x;LL v = power(x, y / 2);v = v * v % mod;if(y % 2 == 1) v = v * x % mod;return v;
}LL cal(LL v){return v * (v  + 1) % mod * (v * 2 + 1) % mod * (v * v * 3 % mod + v * 3 - 1 + mod) % mod * pi % mod;
}void dfs(int cnt, int p, int pos, LL s){if(cnt % 2 == 1) sum = (sum + cal(n / s) * s % mod * s % mod * s % mod * s % mod) % mod;else sum = (sum - cal(n / s) * s % mod * s % mod * s % mod * s % mod + mod) % mod;for(int i = pos; i < p; ++i)dfs(cnt + 1, p, i + 1, s * num[i] % mod);
}int main ()
{//freopen ("in.txt", "r", stdin);get_prime();//for(int i = 0; i < len; ++i) printf("%d ", prime[i]);pi = power(30, mod - 2);int t;scanf ("%d", &t);while (t--){int x, p = 0;scanf("%d", &x);n = x;sum = cal(n);//printf("%d\n", n);for(int i = 0; i < len; ++i){int v = prime[i];if(v > x) break;if(x % v == 0) num[p++] = v;while(x % v ==0) x /= v;}if(x > 1) num[p++] = x;//for(int i = 0; i < p; ++i) printf("%d ", num[i]);for (int i = 0; i < p; ++i)dfs(0, p, i + 1, (LL)num[i]);printf("%lld\n", sum);}return 0;
}