HDU4059 The Boss on Mars【容斥原理】【乘法逆元】【高次求和】

题目链接：

http://acm.hdu.edu.cn/showproblem.php?pid=4059

题目大意：

给一个整数N，求1~N中与N互质的数的4次方的和。

解题思路：

题目简单，过程有点复杂。理清思路就简单了。

利用公式1^4 + 2^4 + … + n^4 = n*(n+1)*(2*n+1)*(3*n^2+3*n-1)/30，可以求出

1^4 + 2^4 + … + n^4，除以30可以先求出30模M的逆元，然后将上式中除以30改为

乘以30的逆元。

再来求与n不互质的数的4次方和。将n质因数分解为 n = p1^a1* p2^a2 * … * pn^an

(其中p1、p2为质数)。

这样不互质的数就是n的质因数的倍数，用容斥原理求出不互质的数的4次方和，最终

(1^4 + 2^4 + … + n^4 - 不互质的数的4次方和)即为所求，具体参考代码。

AC代码：

#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cstring>
#define LL __int64
#define M 1000000007
using namespace std;

LL d,x,y,Inv,N;

void ExGcd(LL a,LL b,LL &d,LL &x,LL &y)
{
if(b == 0)
{
x = 1;
y = 0;
d = a;
}
else
{
ExGcd(b,a%b,d,y,x);
y -= x*(a/b);
}
}

LL ModInverse(LL a,LL n,LL &d,LL &x,LL &y) //a*a-1 = 1(mod n) 用来计算30模M的逆元
{
ExGcd(a,n,d,x,y);
if(1 % d != 0)
return -1;
LL ans = x/d < 0 ? x/d + n : x/d;
return ans;
}

LL Four(LL k) //计算1^4 + 2^4 + … + k^4( 公式为 k*(k+1)*(2*k+1)*(3*k*k+3*k-1)/30 )
{
LL ans = k;
ans = ans * (k+1) % M;
ans = ans * (2*k+1) % M;
ans = ans * ((3*k*k%M+3*k-1)%M) % M;
ans = (ans * Inv) % M; // 除以30 等于 乘以30模M的逆元
return ans;
// return (k%M) * ((k+1)%M) % M * ((2*k+1)%M) % M * ((3*k*k%M + (3*k-1)%M)%M)%M * Inv;
}

LL Factor[20],ct;

void Divide() //将N分解素因子
{
ct = 0;
LL n = N;
for(LL i = 2; i <= sqrt(n*1.0); ++i)
{
if(n % i == 0)
{
Factor[ct++] = i;
while(n % i == 0)
n /= i;
}
}
if(n != 1)
Factor[ct++] = n;
}

LL Cal(LL a)
//计算 质因数乘积为a的所有小于n的倍数的前四次方和，
//即a^4 + (2a)^4 + (3a)^4 + … + (k*a)^4 = a^4 * (1^4 + 2^4 + … + k^4)
{
LL k = N / a;
LL ans = 1;
for(int i = 1; i <= 4; ++i)
ans = ans * a % M;
ans = ans * Four(k) % M;
return ans;
}

LL Solve() //容斥原理求解与n不互质数的4次方和，最终结果为(1^4 + 2^4 + … + n^4 - 与n不互质数的4次方和)
{
Divide();
LL ans = 0;
for(LL i = 1; i < (1 << ct); ++i)
{
LL odd = 0;
LL tmp = 1;
for(LL j = 0; j < ct; ++j)
{
if((1<<j) & i)
{
tmp *= Factor[j];
odd++;
}
}
if(odd & 1)
ans = (ans%M + Cal(tmp)%M) % M;
else
ans = (ans%M - Cal(tmp)%M + M) % M;
}
return (Four(N) - ans + M)%M;
}

int main()
{
int T;
scanf("%d",&T);
while(T--)
{
scanf("%I64d",&N);
if(N == 1)
printf("0\n");
else
{
Inv = ModInverse(30,M,d,x,y); //求30模M的逆元
printf("%I64d\n",Solve());
}
}

return 0;
}