hdu 4059 推公式+容斥

去年大连赛区现场赛的一道题。应该说一看就要用到容斥原理，筛除与n不互质的数即可。。问题是如何求公式1^4+2^4+3^4+...+n^4?这首先要明白1^2+2^2+3^2+...+n^2=1/6n(2n+1)(n+1)是如何来的（不是数学归纳法）。。

(1+1)^3=1^3+3*1*1+3*1*1+1

(2+1)^3=2^3+3*2^2*1+3*2*1+1

...

(n+1)^3=n^3+3*n^2+3*n+1

把上述n个式子累加，整理后得到上述的公式。。同理可推出1^3+2^3+3^3+...+n^3，最后推出1^4+2^4+3^4+...+n^4=(6n^5+15n^4+10n^3-n)/30

还有一个地方必须要注意，就是模的除法，不能像+,-,*那样直接处理，这个比较麻烦。。略掉，还有模的减法是(a-b+mod)%mod，不要把+mod忘掉，否则结果可能为负值

附代码：

#include <iostream>
#include <vector>
 
using namespace std;

#define mod 1000000007
#define pb push_back
#define LL long long
 
const int N=10010;

int a[N];
vector<int>b,c;
LL x1,y1;
const int d[6]={0,-1,0,10,15,6};

LL mul(int x,int t)
{
    LL res=1;
    for (int i=1;i<=t;i++) res=(res*x)%mod;
    return res;
}

void js1(int n,int t)
{
    int i;
    LL x=n/30,y=n%30;
    for (i=1;i<t;i++) 
    {
        x=(x*n)%mod;
        x=x+(y*n)/30;
        y=(y*n)%30;
    }
    x=(x*d[t]+y*d[t]/30)%mod;
    y=(y*d[t])%30;
    y1+=y;
    x1=(x1+x+y1/30)%mod;
    y1=y1%30;
}

LL js(int n)
{
    x1=0;y1=0;
    for (int i=5;i>=1;i--) 
    {
        if (i==2) continue;
        js1(n,i);
    }
    return x1;
}

void deal(int n)
{
    int n1=n,i;
    for (i=0;i<b.size();i++)
    {
        if (b[i]*b[i]>n) break;
        if (n%b[i]==0)
        {
            c.pb(b[i]);
            while (n1%b[i]==0) n1/=b[i];
            if (n1==1) break;
        }
    }
    if (n1!=1) c.pb(n1);
}

LL rc(int n,int i,int k,int r)
{
    int j,t;
    LL x,res=0;
    if (i>c.size()) return 0;
    for (j=k+1;j<c.size();j++)
    {
        t=r*c[j];
        x=mul(t,4)*js(n/t)%mod;
        if (i%2==0) x*=-1;
        res=(res+x+rc(n,i+1,j,t))%mod;
    } 
    return res;
}

int main()
{
    int n,i,j,T;
    b.clear();
    for (i=2;i<N;i++)
        if (!a[i])
        {
            for (j=i*2;j<N;j+=i) a[j]=1;
            b.pb(i);
        }
    scanf("%d",&T);
    while (T--)
    {
        scanf("%d",&n);
        LL zs=js(n);
        c.clear();
        deal(n);
        int res=(zs-rc(n,1,-1,1)+mod)%mod;
        printf("%d\n",res);
    }
    return 0;
}