hdu 4059 The Boss on Mars(纳入和排除)
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; }