hdu 4565 So Easy! /2013 ACM-ICPC 长沙赛区全国邀请赛A题矩阵乘法

最新推荐文章于 2019-07-27 18:29:09 发布

knownothing

最新推荐文章于 2019-07-27 18:29:09 发布

阅读量2.5k

点赞数

分类专栏：矩阵数论* 文章标签：矩阵乘法

本文链接：https://blog.csdn.net/a601025382s/article/details/10045229

版权

数论* 同时被 2 个专栏收录

113 篇文章 0 订阅

订阅专栏

矩阵

38 篇文章 0 订阅

订阅专栏

#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <vector>
#include <iostream>
using namespace std;
#define LL __int64
struct matrix{
    LL f[2][2];
};
LL m;
matrix mul(matrix a,matrix b)
{
    LL i,j,k;
    matrix c;
    memset(c.f,0,sizeof(c.f));
    for(i=0;i<2;i++)
        for(j=0;j<2;j++)
            for(k=0;k<2;k++)
            c.f[i][j]=(c.f[i][j]+a.f[i][k]*b.f[k][j])%m;
    return c;
}
matrix pow_mod(matrix e,LL b)
{
    matrix s;
    s.f[0][0]=s.f[1][1]=1;
    s.f[0][1]=s.f[1][0]=0;
    while(b)
    {
        if(b&1)
            s=mul(s,e);
        e=mul(e,e);
        b=b>>1;
    }
    return s;
}
int main()
{
    LL a,b,n;
    while(scanf("%I64d%I64d%I64d%I64d",&a,&b,&n,&m)!=EOF)
    {
        LL p,q,ans;
        matrix e;
        p=2*a;
        q=a*a-b;
        e.f[0][0]=p;e.f[0][1]=1;
        e.f[1][0]=-q;e.f[1][1]=0;
        e=pow_mod(e,n-1);
        ans=((p*e.f[0][0]+2*e.f[1][0])%m+m)%m;//注意可能是负的结果
        printf("%I64d\n",ans);
    }
    return 0;
}
/*
    因为(a-1)^2<b<a^2,所以0<a-√b<1得到0<(a-√b)^n<1

    因为(a+√b)^n+(a-√b)^n为整数，令x=a+√b,y=a-√b,则p=x+y=2*a,q=x*y=a*a-b所以ceil（(a+√b)^n）=
    (a+√b)^n+(a-√b)^n=x^n+y^n=(x+y)*(x^(n-1)+y^(n-1))-x*y*(x^(n-2)+y^(n-2))

    令f[i]=x^i+y^i，f[0]=2,f[1]=p得到矩阵
    |f[n] f[n-1]|=|f[1] f[0]|*|p  1|^(n-1)
                              |-q 0|
    之后就是矩阵快速幂了
*/