51nod 1537 分解 (矩阵快速幂)

                                                                                                                                                   题目链接

1537 分解

基准时间限制：0.5 秒 空间限制：131072 KB 分值: 80  难度：5级算法题

 收藏

 取消关注

 问(1+sqrt(2)) ^n  能否分解成 sqrt(m) +sqrt(m-1)的形式 

如果可以 输出 m%1e9+7 否则 输出no

Input

一行，一个数n。（n<=10^18）

Output

一行，如果不存在m输出no，否则输出m%1e9+7

Input示例

2

Output示例

9

题解：

首先有个定理：

对于(sqrt(2)+1)^n,一定存在一个m,可以写成：(sqrt(2)+1)^n=sqrt(m)-sqrt(m-1)的形式

基本思路：如果(sqrt(2)-1)^n=sqrt(x)-sqrt(y),那么(sqrt(2)+1)^n=sqrt(x)+sqrt(y)……
于是((sqrt(2)-1)^n)*((sqrt(2)+1)^n)=x-y=1.
证明：把(sqrt(2)-1)^n展开,它可以表示成(b-a*sqrt(2))*(-1)^n.
同样,把(sqrt(2)+1)^n展开,它可以表示成b+a*sqrt(2).
n为奇数时,可令x=2*a^2,y=b^2.n为偶数时,可令x=b^2,y=2*a^2.
于是就和上面一样……
(sqrt(2)-1)^n=sqrt(x)-sqrt(y),(sqrt(2)+1)^n=sqrt(x)+sqrt(y).
x-y=((sqrt(2)-1)^n)*((sqrt(2)+1)^n)=1.

可以构造出：
假设（1-sqrt（2））N-1次为：A+B*sqrt2
则（1-sqrt（2））N次为：A+2*B +（A+B）*sqrt2

矩阵

|1 2 | 
|1 1 |

最后结果为sqrt（m）+sqrt（m-1）
要对m取模---
所以--A+B*sqrt2===sqrt（A*A）+sqrt(2*B*B)

当某一步为（A+10^9+7）+B*sqrt2时---
下一步为（A+10^9+7+2*B）+（B+A+10^9+7）*sqrt2--

等价于取模后---A+B*sqrt2----的
下一步为：（A+2*B）+（B+A）*sqrt2


即对A和B取模不会影响最后的结果。

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<vector>
using namespace std;
typedef long long ll;
const ll mod=1e9+7;

struct Mat {
    ll mat[2][2];
};
Mat operator * (Mat a, Mat b) {
    Mat c;
    memset(c.mat, 0, sizeof(c.mat));
    int i, j, k;
    for(k = 0; k < 2; ++k) {
        for(i = 0; i < 2; ++i) {
            for(j = 0; j < 2; ++j) {
                c.mat[i][j] = (c.mat[i][j] +a.mat[i][k] * b.mat[k][j])%mod;
            }
        }
    }
    return c;
}
Mat operator ^ (Mat a, ll k) {
    Mat c;
    int i, j;
    for(i = 0; i < 2; ++i)
        for(j = 0; j < 2; ++j)
            c.mat[i][j] = (i == j);    //初始化为单位矩阵

    for(; k; k >>= 1) {
        if(k&1) c = c*a;
        a = a*a;
    }
    return c;
}
int main()
{
	ll n;
	scanf("%lld",&n);
	if(n==0)
	{
		puts("1");
		return 0;
	}
	Mat a;
	a.mat[0][0]=a.mat[1][0]=a.mat[1][1]=1;
	a.mat[0][1]=2;
	Mat b=a^(n-1);
	ll aa=(b.mat[0][0]+b.mat[0][1])%mod;
	ll bb=(b.mat[1][0]+b.mat[1][1])%mod;
	if(n&1) printf("%lld\n",2*bb%mod*bb%mod);
	else printf("%lld\n",aa%mod*aa%mod);
	return 0;
}