hdu 4059 The Boss on Mars（容斥原理）

题意：求1-n中与n互质的数的4次方之和，即S=a1^4+a2^4+……; a1,a2……均小于等于n且与n互质。

先求出1^4+2^4+……n^4然后减去与n不互质的数的4次方。赤裸裸的容斥啊！！！！！！

必然要先要用到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，还得算一下30对MOD的逆元，也就是30^(MOD-2),（费马小定理）

接下来要减掉与n不互质的数4次方，将n质因子分解后运用容斥原理即可，就是减掉一个因子的倍数的4次方结果，加上两个因子乘积的倍数的4次方结果，减去……以此类推。

#include<iostream>
#include<cstdio>
#include<cstring>
#include<vector>
#define LL long long
#define MOD 1000000007
using namespace std;
LL res;  //30对MOD的逆元
int prime[10000],cnt=0,flag[10005]= {0};
vector<int>fact;
LL PowMod(LL a,LL b)
{
LL ret=1;
while(b)
{
if(b&1)
ret=(ret*a)%MOD;
a=(a*a)%MOD;
b>>=1;
}
return ret;
}
LL Sum(LL n)  //求an=n^4,的前n项和
{
LL ans=n;
ans=(ans*(n+1))%MOD;
ans=(ans*((2*n+1)%MOD))%MOD;
ans=(ans*(((3*n*n)%MOD+(3*n)%MOD-1+MOD)%MOD))%MOD;
ans=(ans*res)%MOD;
return ans;
}
LL Pow(LL n)  //求n^4
{
LL ans=n;
ans=(((((ans*n)%MOD)*n)%MOD)*n)%MOD;
return ans;
}
int t;
void Prime()   //筛选素数，便于后面的分解
{
for(int i=2; i<=10000; i++)
{
if(flag[i]) continue;
prime[cnt++]=i;
for(int j=2; j*i<=10000; j++)
flag[i*j]=1;
}
}
void Init()
{
res=PowMod(30,MOD-2);   //求30对MOD的逆元
Prime();
scanf("%d",&t);
}
LL dfs(int idx,LL n)     //容斥原理
{
LL ret=0,tmp;
for(int i=idx; i<fact.size(); i++)
{
tmp=fact[i];
ret=(ret+(Sum(n/tmp)*Pow(tmp))%MOD)%MOD;
ret=((ret-dfs(i+1,n/tmp)*Pow(tmp))%MOD+MOD)%MOD;
}
return ret%MOD;
}
int main()
{
LL n;
Init();
while(t--)
{
scanf("%I64d",&n);
fact.clear();
LL tmp=n;
for(int i=0; i<cnt&&prime[i]<=tmp; i++)
if(tmp%prime[i]==0)
{
fact.push_back(prime[i]);
while(tmp%prime[i]==0)
tmp/=prime[i];
}
if(tmp!=1)
fact.push_back(tmp);
LL sum=((Sum(n)-dfs(0,n))%MOD+MOD)%MOD;
printf("%I64d\n",sum);
}
return 0;
}