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

来源：互联网发布：淘宝昵称修改老说非法编辑：程序博客网时间：2024/06/07 22:31

题意：

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

0 0