數論 + 容斥 - HDU 4059 The Boss on MarsThe Boss on Mars Problem's Link

Mean: 

給定一個整數n，求1~n中所有與n互質的數的四次方的和.(1<=n<=1e8)

analyse:

看似簡單，倘若自己手動推公式的話，還是需要一定的數學基礎.

總的思路：先求出sum1=(1^4)+(2^4)+...(n^4)，再求出sum2=（1~n中與n不互質的數的四次方的和），answer=sum1-sum2.

如何求sum1呢？

有兩種方法：

　　1.數列差分.由于A={Sn}={a1^4+a2^4+...an^4}對應一個五階線性差分方程，隻需要求出這個五階線性差分方程的系數即可.

　　2.利用低次幂數和來遞推高次幂數和公式.

最終求得的公式為：Sn=(n*(n+1)*(2n+1)*(3*n*n+3*n-1))/30.

注意，上式中最後有除法，而我們的最終答案要對1e9+7取餘，是以需要求30對1e9+7的模逆元.

由于1e9+7是質數，是以可以直接使用結論：

　　a % m = (b/c)%m

　　a % m = b * c ^(m-2)%m ( m為素數 )

證明：

 　　 b = a * c % m;

則有：b = a % m * c %m;

根據費馬小定理:

　　　a^(p-1)= 1 %p;(p是素數)

可推出：

 　　a%m

　　= a*1%m = a * c^(m-1)%m

　　= a*c*c^(m-2)%m

　　= b*c^(m-2)%m;

-------------------------------------------------------------------------

 求sum2時需要用容斥，當然直接容斥暴力統計的話也會逾時.

注意到：

　　　　2^4+4^4+6^4+8^4 = 2^4*(1^4+2^4+3^4+4^4) .

是以再求sum2時仍然可以使用幂數求和公式，這樣一來時間複雜度就非常低了.

Time complexity: O(logn)

view code

/*

* this code is made by crazyacking

* Verdict: Accepted

* Submission Date: 2015-10-10-22.47

* Time: 0MS

* Memory: 137KB

*/

#include <queue>

#include <cstdio>

#include <set>

#include <string>

#include <stack>

#include <cmath>

#include <climits>

#include <map>

#include <cstdlib>

#include <iostream>

#include <vector>

#include <algorithm>

#include <cstring>

#define max(a,b) (a>b?a:b)

using namespace std;

typedef long long(LL);

typedef unsigned long long(ULL);

const double eps(1e-8);

const int MOD = 1000000007;

typedef long long LL;

int cnt=0;

int a[50];

LL n,inv;

// prime factor decomposition

void primeFactorization(int n)

{

     cnt=0;

     for(int i=2; i*i<=n; i++)

     {

           if(n%i==0)

           {

                 a[cnt++]=i;

                 while(n%i==0)

                       n/=i;

           }

     }

     if(n>1) a[cnt++]=n;

}

LL mulMod(LL a,LL b,LL m)

     LL ans = 0;

     while(b)

           if(b&1)

                 ans = (ans + a)%m;

                 b--;

           b>>=1;

           a=a*2;

           if(a>n) a%=m;

     return ans;

LL quickPow(LL a,LL b,LL m)

     LL ans = 1;

                 ans=mulMod(ans,a,m);

           a=mulMod(a,a,m);

// Exponential sum

LL f(LL n)

     LL ans=n;

     ans=(ans*(n+1))%MOD;

     ans=(ans*(2*n+1))%MOD;

     ans=(ans*((3*n*n+3*n-1)%MOD))%MOD;

     ans=(ans*inv)%MOD;

LL solve(LL n)

     primeFactorization(n);

     for(int i=1; i<(1<<cnt); i++)

           LL tmp = 1;

           int tt = 0;

           for(int j=0; j<cnt; j++)

                 if((1<<j)&i)

                 {

                       tmp=tmp*a[j];

                       tt++;

                 }

           LL q=n/tmp;

           LL t = mulMod(mulMod(tmp,tmp,MOD),mulMod(tmp,tmp,MOD),MOD);

           if(tt&1)

                 ans=(ans+mulMod(f(q),t,MOD))%MOD;

           else

                 ans=(ans-mulMod(f(q),t,MOD)+MOD)%MOD;

int main()

     int t;

     scanf("%d",&t);

     while(t--)

           scanf("%I64d",&n);

           if(n==1)

                 puts("0");

                 continue;

           inv = quickPow(30,MOD-2,MOD);

           LL tmp = f(n);

           LL ans = solve(n);

           ans=(tmp-ans+MOD)%MOD;

           printf("%I64d\n",ans);

     return 0;