【 HDU 4059 The Boss on Mars】 差分序列+容斥
题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=4059
题目大意：给你一个整数n，让你求前面n个数中与n互素的数E（x^4）的和。
解题思路：
可以先求出Sum（n）：表示1^4+2^4+3^4+……+n^4;
然后再求Sum（n）-E(x^4) ,x表示与 n不互素的数。
求与n不互素的数可以利用容斥定理，求Sum（n）可以利用差分。
简单介绍一下差分序列：
给你一列数 a[i][1]，a[i][2]，a[i][3]，a[i][4]，a[i][5]……
那么a[i][j]=a[i-1][j+1]-a[i-1][j]， 即后一行是上一行相邻两项的差（第一行除外）。
如果给你一个多项式， 比如 f（x）=(x+1)*(x+2)*……*（x+p）,即多项式最高项系数为p。
则这个差分表有如下性质：
1、将 0、1、2、……、p-1、p分别代入多项式中，则他们组成多项式的第一行，后面的差分表可以由上一行推出。第p+1行以后差分表的值都为0。
2、我们这里要用的差分序列是他的第0行对角线的数。 我们设他们为c0、c1、c2、……cp；
==> 第n项的值：f（n）=c0*C（n，0）+c1*C(n，1)+c2*C（n，2）+……+cp*（n，p）；
==> 前n项的值：Sum（n）=c0*C（n+1，1）+c1*C(n+1，2)+c2*C（n+1，3）+……+cp*（n+1，p+1）；
这题藐视搞得我蛋疼，我是把求前n项和组合公式给化简出来（Sum（n）=(n^5)/5+(n^4)/2+(n^3)/3-n/30，最后发现求逆元取模的时候问题很多。最后又多写了几个函数，直接求组合数的值，分母和模求逆元的时候打表先求好。
还有要注意的是：
1、这题不打表先求素数会超时。
2、运用容斥定理的时候取模遇减号要加mod，还有容斥定理求出的数可以用到Sum（n）：比如2^4+4^4+6^4+8^4 = 2^4*(1^4+2^4+3^4+4^4) 看懂了木有。
1 #include<iostream> 2 #include<cstdio> 3 #include<cmath> 4 #include<cstring> 5 #include<vector> 6 #define lld long long 7 #define MOD 1000000007 8 using namespace std; 9 int prime[10000],cnt=0,flag[10005]= {0}; 10 vector<int>vt; 11 12 lld ext_gcd(lld a, lld b, lld &x, lld &y) 13 { 14 if( !b ) return x=1,y=0,a; 15 lld d= ext_gcd(b, a%b, y, x); 16 y-= a/b*x; 17 return d; 18 } 19 lld Inv(int a, int m) ///求逆元 20 { 21 lld d,x,y; 22 d= ext_gcd(a, m, x, y); 23 return x<0 ? x+m : x; 24 } 25 26 lld base[8]= {0,0,1,14,36,24}; ///差分表第0行对角线的值 27 lld dd[10]= {0, 1, 500000004, 166666668, 41666667, 808333339}; ///打表求出需要的逆元 28 29 lld Comb(lld n, lld k) ///组合 30 { 31 if(k==n) return 1; 32 if(k>n) return 0; 33 lld s= 1; 34 for(lld i=n,j=1; j<=k; --i,++j) 35 s= s*i%MOD; 36 return s*dd[k]%MOD; 37 } 38 39 lld Sum(lld n) ///求n^4的前n项和 40 { 41 if(n==1) return 1; 42 lld sum= 0; 43 for(int i=2; i<=5; ++i) ///差分 44 { 45 sum+= ( Comb( n+1, i ) * base[i] )%MOD; 46 if( sum>=MOD ) sum-= MOD; 47 } 48 return sum; 49 } 50 51 lld Pow(lld n) ///求n^4 52 { 53 lld ans=n; 54 ans=(((((ans*n)%MOD)*n)%MOD)*n)%MOD; 55 return ans; 56 } 57 58 void Prime() ///筛选素数，便于后面的分解 59 { 60 for(int i=2; i<=10000; i++) 61 { 62 if(flag[i]) continue; 63 prime[cnt++]=i; 64 for(int j=2; j*i<=10000; j++) 65 flag[i*j]=1; 66 } 67 } 68 69 lld dfs(int idx,lld n) ///容斥原理 70 { 71 lld ret=0,tmp; 72 for(int i=idx; i<vt.size(); i++) 73 { 74 tmp=vt[i]; 75 ret=(ret+(Sum(n/tmp)*Pow(tmp))%MOD)%MOD; 76 ret=((ret-dfs(i+1,n/tmp)*Pow(tmp))%MOD+MOD)%MOD; 77 } 78 return ret%MOD; 79 } 80 int main() 81 { 82 int t; 83 lld n; 84 Prime(); 85 scanf("%d",&t); 86 while(t--) 87 { 88 scanf("%I64d",&n); 89 vt.clear(); 90 lld tmp=n; 91 for(int i=0; i<cnt&&prime[i]<=tmp; i++) ///求n的素因子 92 if(tmp%prime[i]==0) 93 { 94 vt.push_back(prime[i]); 95 while(tmp%prime[i]==0) 96 tmp/=prime[i]; 97 } 98 if(tmp!=1) 99 vt.push_back(tmp); 100 101 lld sum=((Sum(n)-dfs(0,n))%MOD+MOD)%MOD; 102 printf("%I64d\n",sum); 103 } 104 return 0; 105 }
