您的位置:首页 > 其它

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

2018-01-15 16:26 381 查看

题目描述:

当 HungryBacon 还是个小学生的时候,体育老师就告诉过他,求四次方程的根是非常简单的。小 HungryBacon 很聪明,并算出了方程 x4+ax3+bx2+cx+d=0的四个根x1, x2, x3, x4。紧接着,他就找他的密友小 TinkEngineer 对答案,他们对答案的方式是计算 xn1 + xn2 + xn3 + xn4。然而,小 TinkEngineer 急着和女朋友约会,所以就让你来帮忙。

给出a,b,c,d求(xn1 + xn2 + xn3 + xn4)mod(109+7)

分析:

这里需要用到韦达定理



很显然我们需要求的就是Sn

知道这些还不够,还必须求出前3项,前三项满足以下公式:

以下证明由数竞大佬周亮提供:

将原方程写为:(x−x1)(x−x2)(x−x3)(x−x4)=0

拆开括号,得x4−(x1+x2+x3+x4)x3+(x1x2+x1x3+x1x4+x2x3+x2x4+x3x4)x2−(x1x2x3+x1x2x4+x1x3x4+x2x3x4)x+x1x2x3x4=0

首先得到x1+x2+x3+x4=−a

(x1+x2+x3+x4)2−2x1x2−2x1x3...=a2−2b=x21+x22+x23+x24

最后,(x21+x22+x23+x24)(x1+x2+x3+x4)=x31+x32+x33+x34+x1x22+x1x23+x1x24+x2x21...

由x1x22+x1x23+x1x24+x2x21...=−ab+3c

最终得到x31+x32+x33+x34=(−a)(a2−2b)+ab−3c=−a3+3ab−3c

综上:x1+x2+x3+x4=−a

x21+x22+x23+x24=a2−2b

x31+x32+x33+x34=−a3+3ab−3c

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

−a−b−c−d100001000010(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
%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);
}
}
}
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: