以下转载于：
http://blog.csdn.net/xieshimao/article/details/6840731
数论题！
求与N不互质的数的K次方（K=4），反过来想若知道与N互质的K次方和，那所求就容易多了哦。
观察到与n互质的数的性质
比如12=2*2*3
那么与12不互质的数就有2,3,4,6,8,9,10,12
其实就是2的所有倍数，以及3的所有倍数
所以可以先求一个1到12的所有数的四次方和。这个有公式：
n*(n+1)*(2*n+1)*(3*n*n+3*n-1)/30
注意对与除以30可以看成是乘以30的逆元（对于1000000007的逆元）。
求的所有的四次方和之后当然要减去那些不互质的数的四次方
也就是说分别剪去了2和3的倍数的四次方，注意这里2和3的公倍数被多减去了，所以要加回来
所以容斥原理也要用
/*
* zoj3547.c
*
* Created on: 2011-10-3
* Author: bjfuwangzhu
*/
#include<stdio.h>
#include<string.h>
#include<math.h>
#define LL long long
#define nmod 1000000007
#define nmax 10010
#define nnum 1500
int prime[nnum], flag[nmax], pfactor[nnum];
int plen, pflen;
LL inverse, x, y;
void extend_gcd(LL a, LL b) {
if (b == 0) {
x = 1, y = 0;
return;
}
extend_gcd(b, a % b);
LL tx = x;
x = y, y = tx - a / b * y;
}
void init() {
int i, j;
memset(flag, -1, sizeof(flag));
for (i = 2, plen = 0; i < nmax; i++) {
if (flag[i]) {
prime[plen++] = i;
}
for (j = 0; (j < plen) && (i * prime[j] < nmax); j++) {
flag[i * prime[j]] = 0;
if (i % prime[j] == 0) {
break;
}
}
}
extend_gcd(30, nmod);
inverse = (x % nmod + nmod) % nmod;
}
LL getSumPow(int n) {
LL res, temp;
res = 1;
res = res * n % nmod;
res = res * (n + 1) % nmod;
res = res * (n * 2 + 1) % nmod;
temp = (LL) n;
temp = ((temp * temp % nmod * 3 + temp * 3 % nmod - 1) % nmod + nmod) % nmod;
res = res * temp % nmod;
res = res * inverse % nmod;
return res;
}
LL getPow(int n) {
LL res;
res = (LL) n;
res = res * res % nmod * res % nmod * res % nmod;
return res;
}
LL dfs(int start, int n) {
LL res;
int i, temp;
for (i = start, res = 0; i < pflen; i++) {
temp = pfactor[i];
res = (res + getSumPow(n / temp) * getPow(temp) % nmod) % nmod;
res = ((res - dfs(i + 1, n / temp) * getPow(temp) % nmod) % nmod + nmod)
% nmod;
}
return res;
}
int main() {
#ifndef ONLINE_JUDGE
freopen("data.in", "r", stdin);
#endif
int t, n, N, i, te;
LL res, temp;
init();
while (~scanf("%d", &t)) {
while (t--) {
scanf("%d", &n);
te = (int) sqrt(n * 1.0);
N = n;
for (i = 0, pflen = 0; (i < plen) && (prime[i] <= te); i++) {
if (n % prime[i] == 0) {
pfactor[pflen++] = prime[i];
while (n % prime[i] == 0) {
n /= prime[i];
}
}
}
if (n > 1) {
pfactor[pflen++] = n;
}
if (n == N) {
printf("%lld\n", getSumPow(n - 1));
continue;
}
temp = dfs(0, N);
res = ((getSumPow(N) - temp) % nmod + nmod) % nmod;
printf("%lld\n", res);
}
}
return 0;
}
以下转载于：
http://hi.baidu.com/tomspirit/blog/item/6b41263b4ef927ef3b87cecf.html
/*题意：找出比n小且与n互质的数的4次方和4: sum=1+3*3*3*3=82解法：
* 容斥定理+状态压缩枚举
* 比如12:sum = Sum(1^4+...12^4)%MOD - (所有非互质数的4次方和)%MOD
* 求解Sum(1^4+...+12^4)有公式：1/30*n*(n+1)(2n+1)(3n^2+3n-1)，
* 因为式子中有除，取余扩展欧几里得来做。现在的问题是所有非互质数的4次方和。
* 可以先求12的素因子：2,3。在非互质数中，
* ①含2素因子的数为2,4,6,8,10；含3素因子的数为3,6,9,12
* ②含2,3的数为6,12可以看出非互质数 = ① - ②；
* 满足容斥定理。对于2(含2因子)分别是2,4,6,8,10,12
* ，转化为题意,2^4+4^4+6^4+8^4+10^4+12^4 = 2^4*(1^4+2^4+3^4+4^4+5^4+6^4)
* 这样就可以利用上面的公式求解，但是这只是一种情况。我们可以用状态压缩来枚举所有的情况。
* 对于12，有4中情况，全不选，只选2，只选3和2,3都选,
* 那么相对于的和为0*(0)，2^4(1^4+...+6^4)，3^4(1^4+...+4^4)，6^4(1^4+...+2^4)
* 这样才把所有的非互质数的4次方和找到，并利用容斥定理（奇数个素因子为+,偶数个素因子为-），
* 答案就出来了。*/
#include<stdio.h>
#include<cstring>
#include<vector>
#include <map>
#include<iostream>
using namespace std;
const long long MOD = 1000000007;
const int N = 10010;
vector<long long> a;
bool f[N];
long long prim[N];
void Extend_euclid(long long a, long long b, long long &x, long long &y) {
if (b == 0) {
x = 1;
y = 0;
return;
}
long long tx, ty;
Extend_euclid(b, a % b, tx, ty);
x = ty;
y = tx - a / b * x;
}
void getPrime() {
long long i, j;
memset(f, 1, sizeof(f));
for (i = 2; i * i <= N; i++)
if (f[i])
for (j = i * i; j <= N; j += i)
f[j] = 0;
f[1] = 0;
long long pt = 0;
for (i = 2; i <= N; i++)
if (f[i])
prim[pt++] = i;
}
long long getCnt(long long nn) {
long long res = 0, i;
for (i = 0; prim[i] * prim[i] <= nn; ++i) {
if (nn % prim[i] == 0) {
a.push_back(prim[i]);
nn /= prim[i];
while (nn % prim[i] == 0) {
nn /= prim[i];
}
if (nn == 1)
break;
}
}
if (nn > 1)
a.push_back(nn);
return a.size();
}
long long difNum, n, sum;
bool bit[10];
long long pow4(long long aa) {
return aa * aa % MOD * aa % MOD * aa % MOD;
}
int main() {
long long t, x, y, len, nn;
int i, j;
getPrime();
Extend_euclid(30, MOD, x, y);
x = (x % MOD + MOD) % MOD;
scanf("%lld", &t);
while (t--) {
a.clear();
scanf("%lld", &n);
if (n == 1) {
printf("0\n");
continue;
}
nn = n;
sum = nn * (nn + 1) % MOD * (2 * nn % MOD + 1) % MOD
* (nn * nn % MOD * 3 % MOD + 3 * nn % MOD - 1) % MOD * x % MOD;
difNum = getCnt(n);
long long Limit = 1 << difNum;
for (i = 0; i < Limit; ++i) {
long long sum1 = 1, cnt = 0;
for (j = 0; j < difNum; ++j) {
if ((1 << j) & i) {
sum1 *= a[j];
++cnt;
}
}
if (cnt)
nn = n / sum1;
else
nn = 0;
if (cnt & 1)
sum = (sum
- pow4(sum1)
* (nn * (nn + 1) % MOD * (2 * nn % MOD + 1)
% MOD
* (nn * nn % MOD * 3 % MOD
+ 3 * nn % MOD - 1) % MOD * x
% MOD) % MOD + MOD) % MOD;
else
sum = (sum
+ pow4(sum1)
* (nn * (nn + 1) % MOD * (2 * nn % MOD + 1)
% MOD
* (nn * nn % MOD * 3 % MOD
+ 3 * nn % MOD - 1) % MOD * x
% MOD) % MOD + MOD) % MOD;
}
printf("%lld\n", sum);
}
return 0;
}