HDU 4059 The Boss on Mars-矩陣+容斥

錯了29遍，終成正果。。。。。

根據題意，很容易的可以想到容斥。

然後的問題就是如何求

sum（n）=1^4+2^4+3^4+....+n^4;

有三種道路：

很顯然：1^4+2^4+3^4+....+n^4=(n^5)/5+(n^4)/2+(n^3)/3-n/30；

則1,用java的大數去敲這個的代碼。

2，用c++敲，但是用到分數取模，求逆元。

3，用c++敲，但是不用這個公式，用矩陣去構造sum(n).

我用的是第三種。但是第三種有的缺陷，就是時間復雜度有點高。

接下來的問題就是如何優化時間復雜度。

可以預處理1~200W的sum，然後剩下的再用矩陣去構造。

具體構造方法看代碼。。

#include
#include
#include
#include
#include
#include
using namespace std;
#define LL __int64
#define MOD 1000000007
#define maxn 2000007
LL yu[maxn];
struct matrix
{
    LL mat[7][7];
    matrix()
    {
        memset(mat,0,sizeof(mat));
    }
    friend matrix operator *(matrix A,matrix B)
    {
        int i,j,k;
        matrix C;
        for(i=1; i<=6; i++)
        {
            for(j=1; j<=6; j++)
            {
                for(k=1; k<=6; k++)
                {
                    C.mat[i][j]=(C.mat[i][j]+(A.mat[i][k]*B.mat[k][j])%MOD)%MOD;
                }
                C.mat[i][j]=C.mat[i][j]%MOD;
            }
        }
        return C;
    }
} ONE,AA[30];
matrix powmul(matrix A,LL k)
{
    matrix B;
    for(int i=1; i<=6; i++)B.mat[i][i]=1;
    int l=1;
    while(k)
    {
        if(k&1)B=B*AA[l];
        l++;
        k>>=1;
    }
    return B;
}
vectorvec;
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=n)return 0;
    LL p=1;
    LL ns=4;
    while(ns--)
    {
        p=p*x;
        p=p%MOD;
    }
    ns=(n-1)/x;
    p=p*kan(ns);
    p=p%MOD;
    return p;
}
void dos(LL n)
{
    LL p=1;
    LL ans=0;
    int i,j,leap;
    LL m=n;
    vec.clear();
    for(i=2;i*i<=n;i++)//求因子的過程
    {
        if(n%i==0)
        {
            vec.push_back(i);
        }
        while(n%i==0)n=n/i;
    }
    if(n!=1)vec.push_back(n);
    n=m;
    int t=1<<(vec.size());
    for(i=1;i