【数学&矩阵加速】hihocoder1555 四次方根

题目描述:

当 HungryBacon 还是个小学生的时候，体育老师就告诉过他，求四次方程的根是非常简单的。小 HungryBacon 很聪明，并算出了方程 $x^4+ax^3+bx^2+cx+d=0$的四个根$x_1, x_2, x_3, x_4$。紧接着，他就找他的密友小 TinkEngineer 对答案，他们对答案的方式是计算 $x_1^n + x_2^n + x_3^n + x_4^n$。然而，小 TinkEngineer 急着和女朋友约会，所以就让你来帮忙。
给出$a,b,c,d$求$(x_1^n + x_2^n + x_3^n + x_4^n)mod (10^9+7)$

---

分析：

这里需要用到韦达定理：

很显然我们需要求的就是$S_n$，由此，我们可以得到一个递推式：
$S_n=-a_{k-1}S_{k-1}-a_{k-2}S_{k-2}-a_{k-3}S_{k-3}-a_{k-4}S_{k-4}$
知道这些还不够，还必须求出前3项，前三项满足以下公式：

以下证明由数竞大佬周亮提供：
将原方程写为：$(x-x_1)(x-x_2)(x-x_3)(x-x_4)=0$
拆开括号，得$x^4-(x_1+x_2+x_3+x_4)x^3+(x_1x_2+x_1x_3+x_1x_4+x_2x_3+x_2x_4+x_3x_4)x^2-(x_1x_2x_3+x_1x_2x_4+x_1x_3x_4+x_2x_3x_4)x+x_1x_2x_3x_4=0$
首先得到$x_1+x_2+x_3+x_4=-a$
$(x_1+x_2+x_3+x_4)^2-2x_1x_2-2x_1x_3...=a^2-2b=x_1^2+x_2^2+x_3^2+x_4^2$
最后，$(x_1^2+x_2^2+x_3^2+x_4^2)(x_1+x_2+x_3+x_4)=x_1^3+x_2^3+x_3^3+x_4^3+x_1x_2^2+x_1x_3^2+x_1x_4^2+x_2x_1^2...$
由$x_1x_2^2+x_1x_3^2+x_1x_4^2+x_2x_1^2...=-ab+3c$
最终得到$x_1^3+x_2^3+x_3^3+x_4^3=(-a)(a^2-2b)+ab-3c=-a^3+3ab-3c$
综上：

$$x_1+x_2+x_3+x_4=-a$$
$$x_1^2+x_2^2+x_3^2+x_4^2=a^2-2b$$
$$x_1^3+x_2^3+x_3^3+x_4^3=-a^3+3ab-3c$$

知道这些开外挂获得的知识后，可以得到一个递推矩阵：

$$\begin{matrix} -a & 1 & 0 & 0\\ -b & 0 & 1 & 0\\ -c & 0 & 0 & 1\\ -d & 0 & 0 & 0 \end{matrix} \tag{1}$$
再将这个矩阵快速幂以求出答案。

#include<cstdio>
#include<cstring>
#include<cmath>
#define SF scanf
#define PF printf
#define MAXN 1010
#define MOD 1000000007
using namespace std;
struct node{
    long long a[5][5];
    node operator * (const node & x) const{
        node b;
        memset(b.a,0,sizeof b.a);
        for(int i=0;i<4;i++)
            for(int j=0;j<4;j++)
                for(int k=0;k<4;k++){
                    b.a[i][j]+=a[i][k]*x.a[k][j];
                    b.a[i][j]%=MOD;
                }
        return b;
    }
}st;
node fsp(long long x){
    if(x==1){
        return st;
    }
    node a=fsp(x/2ll);
    a=a*a;
    if(x%2ll==1)
        a=a*st;
    return a;
}
long long n,s[5],c[5];
int t;
int main(){
    SF("%d",&t);
    while(t--){
        SF("%lld",&n);
        for(int i=1;i<=4;i++)
            SF("%lld",&c[i]);
        c[1]=(MOD+c[1])%MOD;
        c[2]=(MOD+c[2])%MOD;
        c[3]=(MOD+c[3])%MOD;
        c[4]=(MOD+c[4])%MOD;
        for(int i=1;i<=4;i++){
            s[i]=0;
            s[0]=i;
            for(int j=1;j<=i;j++)
                s[i]=(s[i]+c[j]*s[i-j])%MOD;
            s[i]=MOD-s[i];
        }//之前看的一种写法，与题解给出的公式是等效的
        if(n<=4){
            PF("%lld\n",s[n]%MOD);
        }
        if(n>4){
            n-=4;
            memset(st.a,0,sizeof st.a);
            st.a[0][0]=MOD-c[1];
            st.a[0][1]=MOD-c[2];
            st.a[0][2]=MOD-c[3];
            st.a[0][3]=MOD-c[4];
            st.a[1][0]=1;
            st.a[2][1]=1;
            st.a[3][2]=1;
            node res=fsp(n);
            long long ans=0;
            for(int i=1;i<=4;i++)
                ans=(ans+s[i]*res.a[0][4-i])%MOD;
            PF("%lld\n",ans);
        }
    }
}