数论 + 容斥 - 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;
}
