hdu 4059 The Boss on Mars 容斥

题目链接

求出ai^4+a2^4+......an^4的值， ai为小于n并与n互质的数。

用容斥做， 先求出1^4+2^4+n^4的和的通项公式， 显然是一个5次方程， 然后6个方程6个未知数， 我gauss消元解的（雾

然后筛出n所以的因子， 容斥就好。

容斥的时候， 每一个数的贡献是x^4 * getSum(x), getSum就是上面求出的通项公式。

#include <iostream>
#include <vector>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <map>
#include <set>
#include <string>
#include <queue>
#include <stack>
#include <bitset>
using namespace std;
#define pb(x) push_back(x)
#define ll long long
#define mk(x, y) make_pair(x, y)
#define lson l, m, rt<<1
#define mem(a) memset(a, 0, sizeof(a))
#define rson m+1, r, rt<<1|1
#define mem1(a) memset(a, -1, sizeof(a))
#define mem2(a) memset(a, 0x3f, sizeof(a))
#define rep(i, n, a) for(int i = a; i<n; i++)
#define fi first
#define se second
typedef pair<int, int> pll;
const double PI = acos(-1.0);
const double eps = 1e-8;
const int mod = 1e9+7;
const int inf = 1061109567;
const int dir[][2] = { {-1, 0}, {1, 0}, {0, -1}, {0, 1} };
vector <int> p;
ll pow(ll a, int b) {
    ll ret = 1;
    while(b) {
        if(b&1)
            ret = ret*a%mod;
        a = a*a%mod;
        b>>=1;
    }
    return ret;
}
ll getsum(ll n) {
    ll ret = 6*pow(n, 5)+15*pow(n, 4)+10*pow(n, 3)-n;
    ret %= mod;
    ll inv = pow(30LL, mod-2)%mod;
    return ret*inv%mod;
}
void getFactor(ll n) {
    p.clear();
    for(int i = 2; i*i<=n; i++) {
        if(n%i==0) {
            p.pb(i);
            while(n%i==0)
                n/=i;
        }
    }
    if(n>1)
        p.pb(n);
}
int main()
{
    int t;
    ll n;
    cin>>t;
    while(t--) {
        cin>>n;
        ll sum = getsum(n), d = 0;
        getFactor(n);
        int num = p.size();
        for(int i = 1; i<(1<<num); i++) {
            ll multi = 1, one = 0;
            for(int j = 0; j<num; j++) {
                if((1<<j)&i) {
                    one++;
                    multi = multi*p[j]%mod;
                }
            }
            ll tmp = pow(multi, 4)%mod;
            tmp*=getsum(n/multi);
            tmp%=mod;
            if(one&1)
                d += tmp;
            else
                d -= tmp;
            d = (d+mod)%mod;
        }
        cout<<(sum-d+mod)%mod<<endl;
    }
    return 0;
}