hdu 3292 No more tricks, Mr Nanguo（矩阵快速幂解佩尔方程）

题意：

给N和K，求方程：

X^ 2 - N * Y ^ 2 = 1。的第K大的X值，若没有这样的X值存在，则输出balbalbalabalabla...

解析：

佩尔方程的求解。

若N为完全平方数则无解。

若有解，先暴力出佩尔方程的最小特解x1，y1，然后根据迭代公式递推可得矩阵：

            | Xk |    -    | X1    D * y1 |  ^ k - 1      | X1 |

            | Yk |    -    | Y1         X1 |                 | Y1 |

所以由快速幂直接可得第k大解。

代码：

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <algorithm>
#include <cstring>
#include <cmath>
#include <stack>
#include <vector>
#include <queue>
#include <map>
#include <climits>
#include <cassert>
#define LL long long

using namespace std;
const int inf = 0x3f3f3f3f;
const double eps = 1e-8;
const double pi = 4 * atan(1.0);
const double ee = exp(1.0);

const int maxn = 4;
const int mod = 8191;

struct Matrax
{
    LL num[maxn][maxn];
} d, per;

LL n;
LL x, y;
LL D, K;

void searchXY()
{
    y = 1;
    while (1)
    {
        x = (LL) sqrt(D * y * y + 1.0);
        if (x * x - D * y * y == 1)
            break;
        y++;
    }
}

void init()
{
    d.num[0][0] = x % mod;
    d.num[0][1] = D * y % mod;
    d.num[1][0] = y % mod;
    d.num[1][1] = x % mod;
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            per.num[i][j] = (i == j);
        }
    }
}

Matrax add(Matrax a, Matrax b)
{
    Matrax c;
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            c.num[i][j] = (a.num[i][j] + b.num[i][j]) % mod;
        }
    }
    return c;
}

Matrax multi(Matrax a, Matrax b)
{
    Matrax c;
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            c.num[i][j] = 0;
            for (int k = 0; k < n; k++)
            {
                c.num[i][j] += a.num[i][k] * b.num[k][j];
            }
            c.num[i][j] %= mod;
        }
    }
    return c;
}

Matrax power(int k)
{
    Matrax c, t;
    Matrax res = per;
    t = d;
    while (k)
    {
        if (k & 1)
        {
            res = multi(res, t);
            k--;
        }
        else
        {
            k >>= 1;
            t = multi(t, t);
        }
    }
    return res;
}

Matrax sum(int k)
{
    if (k == 1)
        return d;
    Matrax b, t;
    t = sum(k >> 1);
    if (k & 1)
    {
        b = power((k >> 1) + 1);
        t = add(t, multi(t, b));
        t = add(t, b);
    }
    else
    {
        b = power(k >> 1);
        t = add(t, multi(t, b));
    }
    return t;
}

int main()
{
#ifdef LOCAL
    freopen("in.txt", "r", stdin);
#endif // LOCAL
    while (scanf("%I64d%I64d", &D, &K) == 2)
    {
        LL t = sqrt((double)D);
        if (t * t == D)
        {
            printf("No answers can meet such conditions\n");
        }
        else
        {
            n = 2;
            searchXY();
            init();
            d = power(K - 1);
            x = (d.num[0][0] * x % mod + d.num[0][1] * y % mod) % mod;
            printf("%I64d\n", x);
        }
    }
    return 0;
}