HDU

pro： T次询问，每次给出N（N<1e8），求所有Σi^4 (i<=N，且gcd(i,N)==1) ;

sol：  因为N比较小，我们可以求出素因子，然后容斥。  主要问题就是求1到P的4次幂和。  我们知道K次幂和是一个K+1次多项式。 

这里有公式Pre=P*(P+1)*(2P+1)*(3P^2+3P-1)/30; 在不知道公式的情况下，我们可以用拉格朗日差值法求。

1，下面给出DFS容斥的代码 。

#include<bits/stdc++.h>
#define ll long long
#define rep(i,a,b) for(int i=a;i<=b;i++)
using namespace std;
const int maxn=2000010;
const ll Mod=1e9+7;
int p[maxn],tot,N,T,ans,Rev=233333335;
void init()
{
    ll tN=N; tot=0; ans=0;
    for(int i=2;i<=tN/i;i++){
        if(tN%i==0){
            p[++tot]=i;
            while(tN%i==0) tN/=i;
        }
    }
    if(tN>1) p[++tot]=tN;
}
int get(ll y)
{
    ll res=y*(y+1)%Mod;
    res=(y*2+1)%Mod*res%Mod;
    res=(y*y%Mod*3%Mod+(y*3-1)%Mod)*res%Mod;
    res=res*Rev%Mod;
    return res;
    return 1LL*y*(y+1)%Mod*(y*2%Mod+1)%Mod*((3LL*y%Mod*y%Mod+y*3%Mod-1+Mod)%Mod)%Mod*Rev%Mod;
}
void dfs(int pos,int opt,ll sum)
{
    if(pos==tot+1){
        ll t=1LL*sum*sum%Mod*sum%Mod*sum%Mod;
        ll x=get(N/sum);
        x=x*t%Mod;
        if(opt==1) ans=(ans+x)%Mod;
        else ans=((ans-x)%Mod+Mod)%Mod;
        return ;
    }
    dfs(pos+1,opt,sum);
    dfs(pos+1,-opt,sum*p[pos]);
}
int main()
{
    scanf("%d",&T);
    while(T--){
        scanf("%d",&N);
        init();
        dfs(1,1,1LL);
        printf("%d
",ans);
    }
    return 0;
}

2，这里的DFS还可以小小的优化一下， 想象一下，DFS的过程是遍历一颗二叉树，那么它遍历了所有的节点， 而且遍历的过程是老老实实一步步走下去的，所以还可以优化一下。   假设不操作则加入左儿子，有操作进入右儿子。 由于每次我遍历到叶子节点才进行计算， 所以很多时候，我向左走其实进行了一些没有价值的访问，会浪费一些时间。

 而我现在可以在非叶子节点就进行计算，非叶子节点代表的是一直左走代表的叶子（即用节点代表对应的叶子，减少了不必要的访问）。  这样的话，不会存在浪费。 （虽然在此题中未体现出优劣性，但是去年深圳热身赛有一道搜索题就需要这样才能过。

#include<bits/stdc++.h>
#define rep(i,a,b) for(int i=a;i<=b;i++)
using namespace std;
const int maxn=2000010;
const int Mod=1e9+7;
int p[maxn],tot,N,T,ans,Rev=233333335;
void init()
{
    int tN=N; tot=0; ans=0;
    for(int i=2;i*i<=tN;i++){
        if(tN%i==0){
            p[++tot]=i;
            while(tN%i==0) tN/=i;
        }
    }
    if(tN>1) p[++tot]=tN;
}
void dfs(int pos,int opt,int sum) //稍微优化后的DFS
{
    int t=1LL*sum*sum%Mod*sum%Mod*sum%Mod;
    int y=N/sum;
    int x=1LL*y*(y+1)%Mod*(y*2+1)%Mod*(1LL*y*y*3%Mod+y*3-1)%Mod*Rev%Mod;
    if(opt==1) (ans+=1LL*x*t%Mod)%=Mod;
    else ((ans-=1LL*x*t%Mod)+=Mod)%=Mod;
    for(int i=pos+1;i<=tot;i++){
        dfs(i,-opt,sum*p[i]);
    }
}
int main()
{
    scanf("%d",&T);
    while(T--){
        scanf("%d",&N);
        init();
        dfs(0,1,1);
        printf("%d
",ans);
    }
    return 0;
}

3，以及把公式改为拉格朗日差值来求的代码，4次多项式，前缀和为5次多项式，可以通过求前6项来求：（当然这里的拉格朗日还可以优化为O（N），我懒得改了

#include<bits/stdc++.h>
#include<bits/stdc++.h>
#define rep(i,a,b) for(int i=a;i<=b;i++)
using namespace std;
const int maxn=2000010;
const int Mod=1e9+7;
int X[maxn]={0,1,2,3,4,5,6},Y[maxn]={0,1,17,98,354,979,2275};
int p[maxn],tot,N,T,ans,Rev=233333335;
int qpow(int a,int x)
{
    int res=1; while(x){
        if(x&1) res=1LL*res*a%Mod;
        x>>=1; a=1LL*a*a%Mod;
    }
    return res;
}
int Lange(int K)
{
    int res=0;
    rep(i,1,6) {
        int tmp=Y[i];
        rep(j,1,6) {
            if(j==i) continue;
            tmp=1LL*tmp*(K-X[j])%Mod*qpow(X[i]-X[j],Mod-2)%Mod;
            tmp=(tmp+Mod)%Mod;
        }
        (res+=tmp)%=Mod;
    }
    return res;
}
void init()
{
    int tN=N; tot=0; ans=0;
    for(int i=2;i*i<=tN;i++){
        if(tN%i==0){
            p[++tot]=i;
            while(tN%i==0) tN/=i;
        }
    }
    if(tN>1) p[++tot]=tN;
}
void dfs(int pos,int opt,int sum) //稍微优化后的DFS
{
    int t=1LL*sum*sum%Mod*sum%Mod*sum%Mod;
    int y=N/sum;
    int x=Lange(y);
    if(opt==1) (ans+=1LL*x*t%Mod)%=Mod;
    else ((ans-=1LL*x*t%Mod)+=Mod)%=Mod;
    for(int i=pos+1;i<=tot;i++){
        dfs(i,-opt,sum*p[i]);
    }
}
int main()
{
    scanf("%d",&T);
    while(T--){
        scanf("%d",&N);
        init();
        dfs(0,1,1);
        printf("%d
",ans);
    }
    return 0;
}