hdu4059TheBossonMars(差分+容斥原理)

题意：

求小于n (1 ≤ n ≤ 10^8)的数中，与n互质的数的四次方和。

知识点：

    差分：

         一阶差分： 设 

                         则    为一阶差分。

         二阶差分：

           n阶差分：     且可推出   

         性质： 1.

                       2.

    差分序列：

         给你一列数 a[i][1]，a[i][2]，a[i][3]，a[i][4]，a[i][5]……

         那么a[i][j]=a[i-1][j+1]-a[i-1][j]， 即后一行是上一行相邻两项的差（第一行除外）。

         如果给你一个多项式， 比如 f（x）=(x+1)*(x+2)*……*（x+p）,即多项式最高项指数为p。

         则得到的差分序列有如下性质：

         1. f(0),f(1)…f(p)组成多项式的第一行，后面的差分序列可以由上一行推出。第p+1行以后差分序列的值都为0。

         2.我们这里要用的差分序列是其第0行对角线的数。 我们设他们为c0、c1、c2、……cp；   则：

            第n项的值：f（n）=c0*C(n,0)+c1*C(n,1)+c2*C(n,2)+……+cp*C(n,p)；

            前n项的值：Sum（n）=c0*C(n+1,1)+c1*C(n+1,2)+c2*C(n+1,3)+……+cp*C(n+1,p+1)；

            把求前n项和组合公式给化简出来Sum（n）=(n^5)/5+(n^4)/2+(n^3)/3-n/30

                                                                       =(n*(n+1)*(2n+1)*(3*n*n+3*n-1))/30

                  参考文献

题解：反面考虑，容斥原理，sum(n)=1^4+2^4+…n^4=(n*(n+1)*(2n+1)*(3*n*n+3*n-1))/30,减掉与n不互质的数4次方，将n质因子分解后减掉一个因子的倍数的4次方结果，加上两个因子乘积的倍数的4次方结果，减去……以此类推。

其中还涉及逆元，因为MOD为素数，用费马小定理求。

#include
#include
#include
using namespace std;
typedef long long LL;
const LL MOD=1e9+7;
const LL NN=1e8+5;
const int N=1e4+5;
const LL ni=233333335;  //30 mod MOD 的逆
LL n,syz[15],ans;
int ycnt;

LL mutisum(LL n)
{
    LL ans1=1;
    //long long 范围<18,446,744,073,709,551,616 约10^20    30*MOD> LL范围，此法不可，实力被坑
    //LL mod=30*MOD;
    //ans1=(((((n*(n+1))%mod)*(2*n+1))%mod)*((3*n*n+3*n-1)%mod))%mod;
    ans1=(((((((n*(n+1))%MOD)*(2*n+1))%MOD)*((3*n*n+3*n-1)%MOD))%MOD)*ni)%MOD;
    //ans1/=30;
    return ans1%MOD;
}

int prime[N];
bool vis[N];
int pcnt;
void is_prime()
{
    pcnt=0;
    memset(vis,0,sizeof(vis));
    for(int i=2;i    {
        if(!vis[i])
        {
            prime[pcnt++]=i;
            for(int j=i+i;j                vis[j]=1;
        }
    }
}

void fenjie(LL n1)
{
    ycnt=0;
    for(int i=0;i    {
        if(n1%prime[i]==0)
            syz[ycnt++]=prime[i];
        while(n1%prime[i]==0)
            n1/=prime[i];
    }
    if(n1>1)
        syz[ycnt++]=n1;
}

void dfs(int c,int cur,int j,LL ans1)      //dfs(c,1,0,1);
{
    if(cur==c+1)
    {
        LL nn=(n-1)/ans1,as1=ans1%MOD;
        if(c&1)
            ans-=(((((((mutisum(nn)*as1)%MOD)*as1)%MOD)*as1)%MOD)*as1)%MOD;
        else
            ans+=(((((((mutisum(nn)*as1)%MOD)*as1)%MOD)*as1)%MOD)*as1)%MOD;
        ans%=MOD;
        return;
    }
    for(;j    {
        dfs(c,cur+1,j+1,ans1*syz[j]);
    }
}

void test()
{
    for(int i=0;i        cout<    cout<
}

int main()
{
    int t;
    scanf("%d",&t);
    is_prime();
    while(t--)
    {
        scanf("%lld",&n);
        if(n==1)
        {
            printf("1\n");
            continue;
        }
        fenjie(n);
        ans=mutisum(n-1);
        for(int c=1;c<=ycnt;c++)
            dfs(c,1,0,1);
        if(ans<0)
            ans=(ans+MOD)%MOD;
        printf("%lld\n",ans);
       // test();
    }
}

 Case1: sum=1+3*3*3*3=82 Case2: sum=1+2*2*2*2+3*3*3*3+4*4*4*4=354