hdu 4059 The Boss on Mars 容斥原理 数列通项公式

The Boss on Mars

为了高效求出数列的项 $ a_i=\sum\limits^n_{i=1} i^4 $，不使用通项公式是无法实现的。
另外实际分解所得的质因数的种类并不多，无需浪费大量时间、空间筛选大质数。

/** Aug 27, 2015 2:45:37 PM
 * PrjName:hdu4059
 * @author Semprathlon
 */
import java.io.*;
import java.util.*;
public class Main {
    static PrintWriter out=new PrintWriter(System.out);
    static int maxn=105;
    static Vector<Integer> vec=new Vector<Integer>();
    static Vector<Integer> get_prime_factor(int n){
        Vector<Integer> res=new Vector<Integer>();
        for(int i=2;i*i<=n;i++)
            if (n%i==0){
                res.add(i);
                while(n%i==0)
                    n/=i;
            }
        if (n>1) res.add(n);
        return res;
    }
    final static long mod=1000000007L;
    static long pow(long n,long m,long mod){
        long res=1L;
        while(m>0L){
            if ((m&1L)>0L) res=res*n%mod;
            n=n*n%mod;
            m>>=1;
        }
        return res;
    }
    static long div(long a,long b,long mod){
        return a*pow(b,mod-2,mod)%mod;
    }
    static long sum(int n,long mod){
        long res=n;
        res*=2*n+1;res%=mod;
        res*=n+1;res%=mod;
        res*=(3L*n*n+3*n-1)%mod;res%=mod;
        return div(res,30,mod);
    }
    static long solve(Vector<Integer> vec,int n){
        long res=0L;
        int m=vec.size();
        for(int i=1;i<(1L<<m);i++){
            boolean tag=false;
            int tmp=1;
            for(int j=0;j<m;j++)
                if (((1L<<j)&i)>0){
                    tmp*=vec.get(j);
                    tmp%=mod;
                    tag^=true;
                }
            res=res+(tag?1:-1)*(pow(tmp,4,mod)*sum(n/tmp,mod)%mod);
            res%=mod;
        }
        return res;
    }
    public static void main(String[] args) throws IOException{
        // TODO Auto-generated method stub
        InputReader in=new InputReader(System.in);
        int T=in.nextInt();
        while(T-->0){
            int n=in.nextInt();
            Vector<Integer> v=get_prime_factor(n);
            out.println((sum(n,mod)-solve(v,n)+mod)%mod);
        }
        out.flush();
        out.close();
    }

}

另有一个未知的通过矩阵快速幂计算该数列项的方法（求以下矩阵的n次方）。

1	1	4	6	4	1
0	1	4	6	4	1
0	0	1	3	3	1
0	0	0	1	2	1
0	0	0	0	1	1
0	0	0	0	0	1

void init()  
{  
    int i,j;  
    //构造矩阵  
    ONE.mat[1][1]=1;ONE.mat[1][2]=1;ONE.mat[1][3]=4;  
    ONE.mat[1][4]=6;ONE.mat[1][5]=4;ONE.mat[1][6]=1;  
    ONE.mat[2][1]=0;ONE.mat[2][2]=1;ONE.mat[2][3]=4;  
    ONE.mat[2][4]=6;ONE.mat[2][5]=4;ONE.mat[2][6]=1;  
    ONE.mat[3][1]=0;ONE.mat[3][2]=0;ONE.mat[3][3]=1;  
    ONE.mat[3][4]=3;ONE.mat[3][5]=3;ONE.mat[3][6]=1;  
    ONE.mat[4][4]=1;ONE.mat[4][5]=2;ONE.mat[4][6]=1;  
    ONE.mat[5][4]=0;ONE.mat[5][5]=1;ONE.mat[5][6]=1;  
    ONE.mat[6][6]=1;  
    yu[0]=0;  
    yu[1]=1;  
    LL x;  
    for(i=2;i<maxn;i++)//预处理1~200W的sum[x]  
    {  
        x=i%MOD;  
        x=x*i;  
        x=x%MOD;  
        x=x*i;  
        x=x%MOD;  
        x=x*i;  
        x=x%MOD;  
        yu[i]=yu[i-1]+x;  
        yu[i]=yu[i]%MOD;  
    }  
    AA[1]=ONE;  
    for(i=2;i<30;i++)AA[i]=AA[i-1]*AA[i-1];  
}  

最后以这种方式取得计算结果(再乘上以下矩阵)：

1	0	0	0	0	0
1	0	0	0	0	0
1	0	0	0	0	0
1	0	0	0	0	0
1	0	0	0	0	0
1	0	0	0	0	0

matrix A;  
LL kan(LL x)//求sum[x]  
{  
    if(x<maxn)return yu[x];  
    matrix OP;  
    for(int i=1;i<=6;i++)OP.mat[i][1]=1;  
    A=powmul(ONE,x-1);  
    OP=A*OP;  
    return OP.mat[1][1];  
}