HDU 4059 The Boss on Mars ( 容斥原理)

題目連結：

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

題意：

給定一個數n求小於n的與n互斥的數的四次方的和。

分析：

我們可以求出從1~n的所有數的四次方的和sum1，然後容斥求出1~n所有與n不互斥的數的四次方的和sum2；

ans =sum1 - sum2;

設f(n)表示從1~n的所有數的四次方的和 f(n)=1/30*n*(n+1)(2n+1)(3n^2+3n-1);

推倒如下：

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

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

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

...

...

...

2^5-1^5=5*1^4+10*1^3+10*1^2+5*1+1

需要用到的公式有

1^2+2^2+...+n^2=n*(n+1)*(2*n+1)/6;

1^3+2^3+...+n^3=n^2*(n+1)^2/4;

然後全部相加起來，經行化簡 就可以就可以得到f(n)

因為要取模，而且有除法因此我們需要求逆元

a^b=1(mod m)

這裡介紹兩種方法

1)

因為30與1000000007互斥 因此由費馬小定理可知這是的逆元為(30^(mod-2))%(mod);

2)

a/b mod m = (a mod (m*b) )/m;

程式碼如下：

#include <iostream>
#include <cstring>
#include <cstdio>
#include <algorithm>
using namespace std;

const int mod = 1000000007;
typedef long long LL;

int cnt=0;
int a[50];
LL n,inv;
void fen(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 multi(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 quick_mod(LL a,LL b,LL m)
{
    LL ans = 1;
    while(b){
        if(b&1){
            ans=multi(ans,a,m);
            b--;
        }
        b>>=1;
        a=multi(a,a,m);
    }
    return ans;
}

LL f(LL n)// 四次方和公式：1/30*n*(n+1)(2n+1)(3n^2+3n-1)
{
    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)
{
    fen(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 = multi(multi(tmp,tmp,mod),multi(tmp,tmp,mod),mod);
        if(tt&1)
            ans=(ans+multi(f(q),t,mod))%mod;
        else
            ans=(ans-multi(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 = quick_mod(30,mod-2,mod);
        LL tmp = f(n);
        LL ans = solve(n);
        ans=(tmp-ans+mod)%mod;
        printf("%I64d\n",ans);
    }
    return 0;
}