HDU 4059 The Boss on Mars 2011大连现场赛

转载请注明出处，谢谢 http://blog.csdn.net/ACM_cxlove?viewmode=contents           by---cxlove

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

题目求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),证明请移步 http://blog.csdn.net/acm_cxlove/article/details/7422265

貌似前期的公式准备都已经结束了。接下来要减掉与n不互质的数4次方，将n质因子分解后运用容斥原理即可，就是减掉一个因子的倍数的4次方结果，加上两个因子乘积的倍

数的4次方结果，减去……以此类推。也可以通过状态压缩枚举,貌似最多9个质因子吧。

在运算的时候，注意各种相乘溢出就行了，类似计算几何的精度问题，数论的溢出也很纠结。

/*
ID:cxlove
PROB:HDU 4059 
HINT:数论+容斥原理
DATA:2012.4.7
*/
#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;
}