A sequence Sn is defined as:
Sn
=
⌈(a+b√)n⌉
%m
Where a, b, n, m are positive integers.┌x┐is the ceil of x. For example, ┌3.14┐=4. You are to calculate Sn.
You, a top coder, say: So easy!
Input
There are several test cases, each test case in one line contains four positive integers: a, b, n, m. Where 0< a, m < 215, (a-1)2< b < a2, 0 < b, n < 231.The input will finish with the end of file.
Output
For each the case, output an integer Sn.
Sample Input
2 3 1 2013
2 3 2 2013
2 2 1 2013
Sample Output
4
14
4
题目就是问给出a,b,n,m,按公式求解Sn。
直接求精度肯定有误差,我们可以发现(一般发现不了)
(a+b√)n
与
(a−b√)n
加起来就是答案。因为
(a−b√)n
是一个小数(a - 1 <
b2
< a)。
令
sn=(a+b√)n+(a−b√)n
c=a+b√
d=a−b√
即
sn=cn+dn
则可以得到
Sn+1=(cn+dn)(c+d)−cd(cn−1+dn−1)
Sn+1=(c+d)Sn−cdSn−1
即
Sn+1=2aSn−(a2−b)Sn−1
看到这个形式就知道接下来就可以构造矩阵
下面就是矩阵的快速计算了。
注意:结果可能为负的,就是说输出结果+mod再取mod。
#include<cstdio>
#include<cstring>
#include<iostream>
#include<queue>
#include<vector>
#include<algorithm>
#include<string>
#include<cmath>
#include<set>
#include<map>
#include<vector>
#include<stack>
#include<utility>
#include<sstream>
using namespace std;
typedef long long ll;
const int inf = 0x3f3f3f3f;
const int maxn = 1005;
ll m;
struct mat
{
int n,m;
ll a[4][4];
void clr()
{
memset(a,0,sizeof(a));
}
}ini;
mat mul(mat &a,mat &b)
{
mat c;
c.n = c.m = 2;
c.clr();
for(int i = 0;i < a.n;i++)
{
for(int k = 0;k < b.n;k++)
{
for(int j = 0;j < b.m;j++)
c.a[i][j] = (c.a[i][j] + a.a[i][k]*b.a[k][j])%m;
}
}
return c;
}
mat pow_mod(mat& a,ll n)
{
if(n == 0)return ini;
if(n == 1)return a;
mat temp = pow_mod(a,n/2);
temp = mul(temp,temp);
if(n%2)
temp = mul(temp,a);
return temp;
}
int main()
{
#ifdef LOCAL
freopen("C:\\Users\\巍巍\\Desktop\\in.txt","r",stdin);
//freopen("C:\\Users\\巍巍\\Desktop\\out.txt","w",stdout);
#endif // LOCAL
ll a,b,n;
ini.n = ini.m = 2;
ini.clr();
ini.a[0][0] = ini.a[1][1] = 1;
while(scanf("%lld%lld%lld%lld",&a,&b,&n,&m) != EOF)
{
if(n == 1)
{
printf("%lld\n",(2*a + m)%m);
continue;
}
mat t;
t.n = t.m = 2;
t.a[0][0] = 2*a;
t.a[0][1] = -a*a + b;
t.a[1][0] = 1;
t.a[1][1] = 0;
t = pow_mod(t,n - 1);
mat ans;
ans.n = 2;ans.m = 1;
ans.a[0][0] = 2*a;
ans.a[1][0] = 2;
ans = mul(t,ans);
printf("%lld\n",(ans.a[0][0]%m + m)%m);
}
return 0;
}