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;
}