So Easy!

最新推荐文章于 2020-10-27 17:40:50 发布

原创最新推荐文章于 2020-10-27 17:40:50 发布 · 574 阅读

2 ·

CC 4.0 BY-SA版权

数学专栏收录该内容

13 篇文章

订阅专栏

本文介绍了一种高效的方法来解决复杂数学序列计算问题，通过矩阵快速计算技术简化了求解过程，避免了精度损失，并提供了实例代码演示。

部署运行你感兴趣的模型镜像

　A sequence Sn is defined as:

$S_{n}$ = $\lceil(a + \sqrt{b})^{n}\rceil$ %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 + \sqrt{b})^{n}$ 与 $(a - \sqrt{b})^{n}$ 加起来就是答案。因为 $(a - \sqrt{b})^{n}$ 是一个小数（a - 1 < $b^2$ < a)。
令 $s_n = (a + \sqrt{b})^{n} + (a - \sqrt{b})^{n}$
$c = a + \sqrt{b}$
$d = a - \sqrt{b}$
即 $s_n = c^{n} + d^{n}$
则可以得到
$S_{n + 1} = (c^{n} + d^{n})(c + d) - cd(c^{n - 1} + d^{n - 1})$
$S_{n + 1} = (c + d)S_n -cdS_{n - 1}$
即
$S_{n + 1} = 2aS_n - (a^2 - b)S_{n - 1}$

看到这个形式就知道接下来就可以构造矩阵

(S n + 1 S n) = (2 a - a 2 + b 1 0) (S n S n - 1)

$\begin{eqnarray*} \ {S_{n + 1} \choose S_n}={2a\ \ -a^2 + b \choose{ 1 \ \ \ \ \ \ 0}}{S_n \choose S_{n - 1}} \end{eqnarray*}$
下面就是矩阵的快速计算了。
注意：结果可能为负的，就是说输出结果+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;
}