矩阵快速幂--(1+sqrt(2)) ^n 分解成 sqrt(m) +sqrt(m-1)

本文链接：https://blog.csdn.net/winycg/article/details/69808789

证明：对于 $(1+\sqrt{2})^n$ ,一定存在一个m,可以写成： $(1+\sqrt{2})^n=\sqrt{m}+\sqrt{m-1}$ 的形式,n=1,2,3.

$(1+\sqrt{2})^n=\sqrt{x}+\sqrt{y}$ 是必然成立的，那么 $(1-\sqrt{2})^n=\sqrt{x}-\sqrt{y}$ （因为奇数个 $\sqrt{2}$ 相乘的结果也是 $\sqrt{y}$ ）也成立

可得 $(1-\sqrt{2})^n(1+\sqrt{2})^n=(\sqrt{x}-\sqrt{y})(\sqrt{x}+\sqrt{y})=x-y=(-1)^{n}$

所以 $(1+\sqrt{2})^n=\sqrt{m}+\sqrt{m-1}$

$(1+\sqrt{2})^n=a+b*\sqrt{2}$

$(1-\sqrt{2})^n=a-b*\sqrt{2}$

$(1-\sqrt{2})^n(1+\sqrt{2})^n=a^2-2b^2=(-1)^n$

现需要凑出m为某个值，使得 $\sqrt{m}+\sqrt{m-1}=a+b*\sqrt{2}$

①当n为偶数， $a^2-2b^2=1$

可凑得 $m=a^2,m-1=a^2-1=2b^2$ 时，

原等式成立

②当n为奇数时， $a^2-2b^2=-1$

可凑得 $m=a^2+1=2b^2,m-1=a^2$

原等式成立

设 $(1+\sqrt{2})^{n}=a(n)+b(n)*\sqrt{2}$

$(1+\sqrt{2})^{n}=(1+\sqrt{2})^{n-1}(1+\sqrt{2})=a(n-1)+\sqrt{2}(b(n-1)+a(n-1))+2b(n-1)$

得到：

$a(n)=a(n-1)+2b(n-1)$

$b(n)=a(n-1)+b(n-1)$

$\begin{pmatrix} a(n)\\ b(n) \end{pmatrix} =\begin{pmatrix} 1 & 2\\ 1& 1 \end{pmatrix} \begin{pmatrix} a(n-1)\\ b(n-1) \end{pmatrix}$

#include<iostream>
#include<cmath>
#include<cstring>
#include<cstdio>
#include<vector>
#include<stack>
#include<queue>
#include<algorithm>
#include<sstream>
#define inf 0x3f3f3f3f
#define ll long long
using namespace std;
const int mod=1000000007;
struct matix
{
    ll mat[3][3];
};
matix multiply(matix a,matix b)
{
    matix c;
    memset(c.mat,0,sizeof(c.mat));
    for(int i=1;i<=2;i++)
        for(int j=1;j<=2;j++)
            for(int k=1;k<=2;k++)
                c.mat[i][j]=(c.mat[i][j]+(a.mat[i][k]*b.mat[k][j])%mod)%mod;
    return c;
}
int main()
{
    ll n;
    while(cin>>n)
    {
        matix ans;
        ans.mat[1][1]=1;
        ans.mat[1][2]=2;
        ans.mat[2][1]=1;
        ans.mat[2][2]=1;
        matix s;
        s.mat[1][1]=s.mat[2][2]=1;
        s.mat[1][2]=s.mat[2][1]=0;
        ll b=n-1;
        while(b)
        {
            if(b%2==1)
                s=multiply(s,ans);
            ans=multiply(ans,ans);
            b=b/2;
        }
        ll t=(s.mat[1][1]+s.mat[1][2])%mod;
        if(n%2==1)
            cout<<(ll)(t*t+1)%mod<<endl;
        else
            cout<<(ll)(t*t)%mod<<endl;
    }
    return 0;
}