hdu 4059 数论+高次方求和+容斥原理

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

现场赛中通过率挺高的一道题 但是容斥原理不怎么会。。

参考了http://blog.csdn.net/acm_cxlove/article/details/7434864

1、求逆元   p=1e9+7是素数，所以由 a^(p-1)%p同余于1 可得a%p的逆元为a^(p-2)

2、segma(i^k)都可以通过推导得到求和公式 详见http://blog.csdn.net/acm_cxlove/article/details/7434864

3、容斥原理，还在恶补中  代码写的挺漂亮

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#include <set>
#include <cmath>
#include <vector>
using namespace std;
#define ll long long
#define IN(s) freopen(s,"r",stdin)
#define OUT(s) freopen(s,"w",stdout)

const int MOD  =  1000000007;
const int N = 10005;
const int M = 10050;
ll n,thr;//thr 30的逆元
vector<int>fac;

bool is[N]; int prm[M];
int getprm(int n){
    int i, j, k = 0;
    int s, e = (int)(sqrt(0.0 + n) + 1);
    memset(is, 1, sizeof(is));
    prm[k++] = 2; is[0] = is[1] = 0;
    for(i = 4; i < n; i += 2) is[i] = 0;
    for(i = 3; i < e; i += 2) if(is[i]) {
            prm[k++] = i;
            for(s = i * 2, j = i * i; j < n; j += s)
                is[j] = 0;
// 因为j是奇数，所以+奇数i后是偶数，不必处理！
        }
    for( ; i < n; i += 2) if(is[i]) prm[k++] = i;
    return k;  // 返回素数的个数
}

ll qmod(ll x,ll t)
{
    ll ret=1,base=x;
    while(t)
    {
        if(t&1)ret=(ret*base)%MOD;
        base=(base*base)%MOD;
        t/=2;
    }
    return ret;
}

ll sum(ll x)
{
    ll ret=1;
    ret=(ret*x)%MOD;
    ret=(ret*(x+1))%MOD;
    ret=(ret*((2*x+1)%MOD))%MOD;
    ret=(((3*x*x)%MOD+(3*x)%MOD-1+MOD)%MOD*ret)%MOD;
    return (ret*thr)%MOD;
}

inline ll four(ll x)
{
    return (((x%MOD)*x%MOD)*x%MOD)*x%MOD;
}

ll dfs(int cur, ll tmp)//容斥原理
{
    ll ret=0,f;
    for(int i=cur;i<fac.size();i++)
    {
        f=fac[i];
        ret=(ret+(sum(tmp/f)*four(f))%MOD)%MOD;
        ret=( (ret-dfs(i+1,tmp/f)*four(f))%MOD+MOD )%MOD;
    }
    return ret%MOD;
}

int main()
{
    //IN("hdu4059.txt");
    int ncase;
    ll s1,s2;
    scanf("%d",&ncase);
    int prmnum=getprm(N-1);
    thr=qmod(30,MOD-2);
    while(ncase--)
    {
        scanf("%I64d",&n);
        //int sn=(int)sqrt(n);
        fac.clear();
        ll tmp=n;
        for(int i=0;i<prmnum && prm[i]<=tmp;i++)
        {
            if(tmp%prm[i] == 0)
            {
                fac.push_back(prm[i]);
                while(tmp%prm[i] == 0)
                    tmp/=prm[i];
            }
        }
            //while(tmp%prm[i] == 0)fac.push_back(prm[i]),tmp/=prm[i];
        if(tmp!=1)fac.push_back(tmp);
        //cout << "FUck= " << sum(n) << endl;
        printf("%I64d\n",( (sum(n)- dfs(0,n)+MOD)%MOD + MOD)%MOD );
    }
    return 0;
}