求出(sum_{i = 1}^n [gcd(n, i) == 1] i ^ 4)
即求出([1,n])里面与(n)互质的数的四次方和
考虑互质,想到唯一分解求出n的所有质因子,然后利用这些质因子容斥求出所有不符合的数字的四次方和
然后求出所有的和,减去不符合的值就是答案。
[sum_{i=1}^{n} i^{4}=frac{n(n+1)(2 n+1)left(3 n^{2}+3 n-1
ight)}{30}=frac{1}{5} n^{5}+frac{1}{2} n^{4}+frac{1}{3} n^{3}-frac{1}{30} n
]
在求一个数字的倍数的四次方和时,比如2的倍数,把(2^4)提取出来就行了
(
2 4 6 8 \
2^4 + (2*2)^4 + (2*3)^4 + (2*4)^4\
= 2^4 * (1^4 + 2^4 + 3^4 + 4^4)
)
#include <iostream>
#include <cstdio>
#include <vector>
#include <cmath>
#define ll long long
using namespace std;
const int MOD = 1e9 + 7;
std::vector<int> divi(int n){
std::vector<int> v;
for(int i = 2; i * i <= n; i++) {
if(n % i == 0) {
v.push_back(i);
while(n % i == 0) n /= i;
}
}
if(n > 1) v.push_back(n);
return v;
}
int gcd(int a, int b){
return b == 0 ? a : gcd(b, a % b);
}
int lcm(int a, int b){
return a / gcd(a, b) * b;
}
ll pow(ll a, ll b, ll p){
ll ans = 1; a %= p;
while(b){
if(b & 1) ans = ans * a % p;
a = a * a % p;
b >>= 1;
}
return ans;
}
ll inv(ll a, ll p){
return pow(a, p - 2, p);
}
ll cal(ll n){
return n * (n + 1) % MOD * (2 * n + 1) % MOD * (3 * n * n % MOD + 3 * n - 1) % MOD * inv(30, MOD) % MOD;
}
namespace IO {
template <typename T>
inline void w(T x) { if (x > 9) w(x / 10); putchar(x % 10 + 48); }
template <typename T>
inline void write(T x, char c) { if(x < 0) putchar('-'), x = -x; w(x); putchar(c); }
template <typename T>
inline void read(T &x) {
x = 0; T f = 1; char c = getchar();
for (; !isdigit(c); c = getchar()) if (c == '-') f = -1;
for (; isdigit(c); c = getchar()) x = (x << 1) + (x << 3) + (c ^ 48);
x *= f;
}
};
void solve(){
int n; scanf("%d", &n);
ll ans = 0, ans2 = cal(n);
std::vector<int> v = divi(n);
for(int i = 1; i < (1 << v.size()); i++) {
int now = 1, k = 0;
for(int j = 0; j < v.size(); j++) {
if(i & (1 << j)) {
now = now * v[j]; k++;
}
}
now = pow(now, 4, MOD) * cal(n / now) % MOD;
if(k % 2 == 0) {
ans = ((ans - now) % MOD + MOD) % MOD;
}else {
ans = ((ans + now) % MOD + MOD) % MOD;
}
}
IO::write(((ans2 - ans) % MOD + MOD) % MOD, '
');
}
int main(){
int t; scanf("%d", &t);
for(int i = 1; i <= t; i++) solve();
return 0;
}