hdu 45092011年大连现场赛I题（容斥原理）

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=4059；
题意：求小于n与n互质的数的四次方和。
分析：我们要先求解四次方和公式，简单的推一下：

必然要先要用到4次方的求和公式。接下来简单的证明一下，这里前提是你知道3次方的公式，如果不会照下面的模式可以利用2次公式推出3次公式

（x+1)^5=x^5+5*x^4+10*x^3+10*x^2+5*x+1;

则 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+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^4)=(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的逆元，我是手算的。。。正常用一次扩展欧几里得就可以了。

对于与n互质的数，就先求与n不互质的数，就是其素因子的倍数和素因子积的倍数。这里要用容斥求解。

容斥原理详见：容斥原理；

代码如下：

#include <set>
#include <map>
#include <stack>
#include <queue>
#include <math.h>
#include <vector>
#include <string>
#include <utility>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <iostream>
#include <algorithm>
#include <functional>

using namespace std;
long long ans;
const long long mod=1000000007;
const long long maxn=50;
long long fac[maxn];
long long cnt;
void fen(long long n){
    for(long long i=2; i*i<=n; i++){
        if(n%i==0){
            fac[cnt]=i;
            while(n%i==0)n/=i;
            cnt++;
        }
    }
    if(n>1)fac[cnt++]=n;
}//素因子分解
void dfs(long long id,long long k,long long w,long long sumsum){
    if(id==cnt+1)return;
    long long temp;
    for(long long i=w;i<cnt;i++){
        temp=sumsum*fac[i];
        long long m=k/temp;
        long long tmp=(m*(m+1))%mod;
        tmp=(tmp*(2*m+1))%mod;
        tmp=(tmp*(3*m*m%mod+3*m%mod-1+mod))%mod;
        tmp=((tmp*233333335ll)%mod+mod)%mod;
        tmp=(tmp*((((temp*temp%mod)*temp%mod)*temp)%mod))%mod;
        if(id&1)ans=((ans-tmp)%mod+mod)%mod;
        else ans=(ans+tmp)%mod;
        dfs(id+1,k,i+1,temp%mod);
    }
}//递归求解容斥原理
int main()
{
    long long t;
    scanf("%I64d",&t);
    while(t--){
        long long n;
        scanf("%I64d",&n);
        cnt=0;
        fen(n);
        ans=(n*(n+1))%mod;
        ans=(ans*(2*n+1))%mod;
        ans=(ans*(3*n*n%mod+3*n%mod-1+mod))%mod;
        ans=((ans*233333335ll)%mod+mod)%mod;
        dfs(1,n,0,1);
        printf("%I64d\n",(ans%mod+mod)%mod);
    }
    return 0;
}