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

The Boss on Mars

为了高效求出数列的项ai=∑i=1ni4=14+24+34+...+n4=? a_i=\sum\limits^n_{i=1} i^4=1^4+2^4+3^4+...+n^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
0	1	4	6	4
0	0	1	3	3
0	0	0	1	2
0	0	0	0	1
0	0	0	0	0

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
1	0	0	0	0
1	0	0	0	0
1	0	0	0	0
1	0	0	0	0
1	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];
}