http://acm.hdu.edu.cn/showproblem.php?pid=4059

定义S = 1^4 + 2^4 + 3^4+.....+n^4。如今减去与n互质的数的4次方。问共降低了多少。

容斥原理。能够先把与n不互质的数的4次方求出来。那就先对n进行质因子分解，对质因子的组合运用容斥原理。质因子个数为奇数就加，偶数就减。事实上与求[1,n]内与n互质的数的个数类似，该题重点是计算，防止乘法溢出。

对于求解1^4 + 2^4 + 3^4+.....+n^4，能够先类比1^2+2^2+...+n^2的求法，那么求4次方，

首先(n+1)^5= n^5 + 5*n^4 + 10*n^3 + 10*n^2 + 5*n^1 + 1.

那么2^5 = (1+1)^5 = 1^5 + 5*1^4 + 10*1^3 + 10*1^2 + 5*1^1 + 1.

3^5 = (2+1)^5 = 2^5 + 5*2^4 + 10*2^3 + 10*2^2 + 5*2^1 + 1.

........

........

(n+1)^5 = n^5 + 5*n^4 + 10*n^3 + 10*n^2 + 5*n^1 + 1.

将上述全部等式相加。两边抵消同样项，得到(n+1)^5 = 5*(1^4+2^4+……n^4)+10*(1^3+2^3+……+n^3)+10*(1^2+2^2+……+n^2)+5*(1+2+……+n)+n+1，

将1^3+2^3+……+n^3 = (n+1)^2*n^2/4和1^2+2^2+……+n^2
= (n*(n+1)*(2*n+1))/6带入上式，化简得到：

1^4+2^4+……n^4
= (n*(n+1)*(2n+1)*(3*n*n+3*n-1))/30。

由于要取余，要求30对1000000007的逆元，用扩展欧几里得就可以。

#include <stdio.h>

#include <iostream>

#include <map>

#include <set>

#include <bitset>

#include <list>

#include <stack>

#include <vector>

#include <math.h>

#include <string.h>

#include <queue>

#include <string>

#include <stdlib.h>

#include <algorithm>

#define LL __int64

//#define LL long long

#define eps 1e-9

#define PI acos(-1.0)

using namespace std;

const int maxn = 10010;

const LL mod = 1000000007;

LL n;

int fac[maxn];

int facCnt;

int prime[maxn];

LL ni,nii;

//求30对mod的逆元。

LL extend_gcd(LL a, LL b, LL &x, LL &y)

{

    if(b == 0)

    {

        x = 1;

        y = 0;

        return a;

    }

    LL d = extend_gcd(b,a%b,x,y);

    LL t = x;

    x = y;

    y = t-a/b*y;

    return d;

}

LL pow_4(LL t)

{

    LL anw =( ((t*(t+1))%mod*(2*t+1)%mod) * (((3*t*t)%mod+(3*t)%mod-1+mod)%mod )%mod*ni)%mod;

    return anw;

}

LL cal(LL m)

{

    LL t = n/m;

    LL anw1 = m;

    anw1 = (anw1*m)%mod;

    anw1 = (anw1*m)%mod;

    anw1 = (anw1*m)%mod;

    LL anw2 = pow_4(t);

    LL anw = (anw1*anw2)%mod;

    return anw;

}

void getPrime()

{

    bool flag[maxn];

    memset(flag,false,sizeof(flag));

    prime[0] = 0;

    for(int i = 2; i < maxn; i++)

    {

        if(flag[i] == false)

        {

            prime[++prime[0]] = i;

            for(int j = 1; j <= prime[0]&&prime[j]*i < maxn; j++)

            {

                flag[prime[j]*i] = true;

                if(i%prime[j] == 0)

                    break;

            }

        }

    }

}

void getFac()

{

    facCnt = 0;

    LL tmp = n;

    for(int i = 1; i <= prime[0] && prime[i]*prime[i] <= tmp; i++)

    {

        if(tmp % prime[i] == 0)

        {

            fac[facCnt++] = prime[i];

            while(tmp % prime[i] == 0)

                tmp /= prime[i];

        }

        if(tmp == 1) break;

    }

    if(tmp > 1)

        fac[facCnt++] = tmp;

}

int main()

{

    int test;

    scanf("%d",&test);

    getPrime();

    extend_gcd(30,mod,ni,nii);

    while(test--)

    {

        scanf("%I64d",&n);

        getFac();

        LL ans = 0;

        for(int i = 1; i < (1<<facCnt); i++)

        {

            LL mul = 1;

            int cnt = 0;

            for(int j = 0; j < facCnt; j++)

            {

                if(i&(1<<j))

                {

                    cnt++;

                    mul *= fac[j];

                }

            }

            if(cnt&1)

                ans = (ans + cal(mul) )%mod;

            else

                ans = (ans - cal(mul) )%mod;

        }

        ans = ((pow_4(n) - ans)%mod+mod)%mod; //减法时取余

        printf("%I64d\n",ans);

    }

    return 0;

}