数论 + 容斥 - HDU 4059 The Boss on Mars

The 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}对应一个五阶线性差分方程，只需要求出这个五阶线性差分方程的系数即可.

　　有关数列差分求幂数和通项的知识，click here.

　　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;
     while(b)
     {
           if(b&1)
           {
                 ans=mulMod(ans,a,m);
                 b--;
           }
           b>>=1;
           a=mulMod(a,a,m);
     }
     return ans;
}
// 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;
     return ans;
}

LL solve(LL n)
{
     primeFactorization(n);
     LL ans = 0;
     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;
     }
     return ans;
}

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;
}