程式師世界 >> 編程語言 >> C語言 >> C++ >> C++入門知識 >> hdu 4059 The Boss on Mars(容斥)

hdu 4059 The Boss on Mars(容斥)

編輯：C++入門知識

hdu 4059 The Boss on Mars(容斥)

定義S = 1^4 + 2^4 + 3^4+.....+n^4，現在減去與n互質的數的4次方，問共減少了多少。

容斥原理，可以先把與n不互質的數的4次方求出來。那就先對n進行質因子分解，對質因子的組合運用容斥原理，質因子個數為奇數就加，偶數就減。其實與求[1,n]內與n互質的數的個數類似，該題重點是計算，防止乘法溢出。

對於求解1^4 + 2^4 + 3^4+.....+n^4，可以先類比1^2+2^2+...+n^2的求法，那麼求4次方，

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

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

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

........

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

將上述所有等式相加，兩邊抵消相同項，得到(n+1)^5 = 5*(1^4+2^4+……n^4)+10*(1^3+2^3+……+n^3)+10*(1^2+2^2+……+n^2)+5*(1+2+……+n)+n+1，

將1^3+2^3+……+n^3 = (n+1)^2*n^2/4和1^2+2^2+……+n^2 = (n*(n+1)*(2*n+1))/6帶入上式，化簡得到：

1^4+2^4+……n^4 = (n*(n+1)*(2n+1)*(3*n*n+3*n-1))/30。

因為要取余，要求30對1000000007的逆元，用擴展歐幾裡得即可。

#include 
#include 
#include
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#include 
#define LL __int64
//#define LL long long
#define eps 1e-9
#define PI acos(-1.0)
using namespace std;
const int maxn = 10010;
const LL mod = 1000000007;
LL n;
int fac[maxn];
int facCnt;
int prime[maxn];
LL ni,nii;

//求30對mod的逆元。
LL extend_gcd(LL a, LL b, LL &x, LL &y)
{
    if(b == 0)
    {
        x = 1;
        y = 0;
        return a;
    }
    LL d = extend_gcd(b,a%b,x,y);
    LL t = x;
    x = y;
    y = t-a/b*y;
    return d;
}

LL pow_4(LL t)
{
    LL anw =( ((t*(t+1))%mod*(2*t+1)%mod) * (((3*t*t)%mod+(3*t)%mod-1+mod)%mod )%mod*ni)%mod;
    return anw;
}

LL cal(LL m)
{
    LL t = n/m;
    LL anw1 = m;
    anw1 = (anw1*m)%mod;
    anw1 = (anw1*m)%mod;
    anw1 = (anw1*m)%mod;
    LL anw2 = pow_4(t);
    LL anw = (anw1*anw2)%mod;
    return anw;
}

void getPrime()
{
    bool flag[maxn];
    memset(flag,false,sizeof(flag));
    prime[0] = 0;
    for(int i = 2; i < maxn; i++)
    {
        if(flag[i] == false)
        {
            prime[++prime[0]] = i;
            for(int j = 1; j <= prime[0]&&prime[j]*i < maxn; j++)
            {
                flag[prime[j]*i] = true;
                if(i%prime[j] == 0)
                    break;
            }
        }
    }
}

void getFac()
{
    facCnt = 0;
    LL tmp = n;
    for(int i = 1; i <= prime[0] && prime[i]*prime[i] <= tmp; i++)
    {
        if(tmp % prime[i] == 0)
        {
            fac[facCnt++] = prime[i];
            while(tmp % prime[i] == 0)
                tmp /= prime[i];
        }
        if(tmp == 1) break;
    }
    if(tmp > 1)
        fac[facCnt++] = tmp;
}

int main()
{
    int test;
    scanf(%d,&test);
    getPrime();
    extend_gcd(30,mod,ni,nii);
    while(test--)
    {
        scanf(%I64d,&n);
        getFac();
        LL ans = 0;

        for(int i = 1; i < (1<