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

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=4059

题意：给出n，求所有小于等于n的数中与n互质的数字的四次方之和，模mod=100000007

思路：总的思路，求出1到n的所有数字的四次方之和，减去和n不互质的数字的四次方之和。因此，本题需要解决以下几个问题：

（1）求1到n的四次方和的通项公式f4(n)：f4(n)=n*(n+1)*(n*2+1)*(3*n*n+3*n-1)/30。推导如下：

(x+1)^5=x^5+5x^4+10x^3+10x^2+5x+1,

因此：1^5=1

      2^5=(1+1)^5=1^5+5*1^4+10*1^3+10*1^2+5*1+1

      3^5=(2+1)^5=2^5+5*2^4+10*2^3+10*2^2+5*2+1

      ……

  (n+1)^5=(n+1)^5=n^5+5*n^4+10*n^3+10*n^2+5*n+1

相加得：f5(n+1)=f5(n)+5f4(n)+10f3(n)+10f2(n)+5f1(n)+n+1

即：(n+1)^5=5f4(n)+10f3(n)+10f2(n)+5f1(n)+n+1

其中 f3(n)=n^2*(n+1)^2/4

     f2(n)=n*(n+1)*(2n+1)/6

     f1(n)=n*(n+1)/2

带入整理得：f4(n)=n*(n+1)*(n*2+1)*(3*n*n+3*n-1)/30

 

（2）求出30对mod的逆元m:m=30^(mod-2)%mod。推到如下：

设b=30*a，则b%mod=a%mod*30%mod,有费马小定理x^(p-1)%p=1(p是素数且x不能整除p)可得:

30^(mod-1)%mod=1,所以：a%mod=a*1%mod=a*30^(mod-1)%mod=a*30*30^(mod-2)%mod=b*30^(mod-1)%mod

即：a=b*m%mod。

 

（3）容斥原理求解与n不互质的所有数的四次方之和.假设n有三个质因数2,3,5，则DFS求解的基本过程如下：

sum=0;

sum+=2的倍数-2,3的倍数+2,3,5的倍数-2,5的倍数；

sum+=3的倍数-3,5的倍数；

sum+=5的倍数。

 

 #include <iostream>
 #include <cmath>
 #include <cstdio>
 #include <cstring>
 #define int64 __int64
 using namespace std;
 
 const int64 mod=1000000007;
 int64 factor[35],factorNum,n,m;
 
 void init1()
 {
     int64 x=30,y=mod-2;
     m=1;
     while(y)
     {
         if(y&1) m=m*x%mod;
         x=x*x%mod;
         y>>=1;
     }
 }
 
 int64 prime[60000],primeNum,isprime[60000];
 
 void init2()
 {
     int64 i,j;
     for(i=2;i<60000;i++)
     {
         if(!isprime[i]) prime[++primeNum]=i;
         for(j=1;j<=primeNum&&prime[j]*i<60000;j++)
         {
             isprime[prime[j]*i]=1;
             if(i%prime[j]==0) break;
         }
     }
 }
 
 int C;
 
 void cal_factor()
 {
     int64 i,p=n;
     factorNum=0;
     for(i=1;i<=primeNum&&prime[i]<=p;i++) if(p%prime[i]==0)
     {
         factor[++factorNum]=prime[i];
         while(p%prime[i]==0) p/=prime[i];
     }
     if(p>1) factor[++factorNum]=p;
 }
 
 inline int64 P(int64 x)
 {
     int64 ans=x%mod;
     ans=ans*(x+1)%mod;
     ans=ans*(x+x+1)%mod;
     ans=ans*((3*x*x+3*x-1)%mod)%mod;
     ans=ans*m%mod;
     return ans;
 }
 
 inline int64 Q(int64 x)
 {
     return x*x%mod*x%mod*x%mod;
 }
 
 int64 DFS(int64 n,int64 start,int64 end,int64 factor[])
 {
     int64 ans=0,i,x,t;
     for(i=start;i<=end;i++)
     {
         t=factor[i];
         ans=(ans+Q(t)*P(n/t)%mod)%mod;
         x=Q(t)*DFS(n/t,i+1,end,factor)%mod;
         ans=((ans-x)%mod+mod)%mod;
     }
     return ans;
 }
 
 int main()
 {
     init1();
     init2();
     for(scanf("%d",&C);C--;)
     {
         scanf("%I64d",&n);
         cal_factor();
         int64 ans;
         ans=P(n)-DFS(n,1,factorNum,factor);
         ans=(ans%mod+mod)%mod;
         printf("%I64d\n",ans);
     }
     return 0;
 }

 

　　

 

转载于:https://www.cnblogs.com/jianglangcaijin/archive/2012/10/01/2709778.html