[HDU3292] 求第K大数

---

题目描述

　　战国时期齐国君主非常喜欢管类乐器，尤其喜欢竽，他有自己的一支奏竽的乐队，其中乐队的人数是一个正方形数（每行每列都有X个人）。每天下午他都会让自己的乐队排成正方形为自己演奏。
　　国王不知道，其中有一个叫南郭的人，其实根本不会吹竽，也混在其中。国王去世后他的儿子继承父位，新国王不喜欢父亲的欣赏方式，他打算把父亲的竽乐队分成若干几个正方形数，然后分别为他演奏，南郭先生听说后非常害怕，于是偷偷逃走了。结果国王重排乐队之后，发现剩下的人正好分组Y^2人的N组。

---

输入格式

输入有多组数据，每组占一行包含两个正整数N(2<=N<29),K(K<10^9)

---

输出格式

输出第K大的解满足条件的X mod 8191的值，如果没有这样的X存在，则输出“No Answers”。

---

样例数据

样例输入

2 999888
3 1000001
4 8373

样例输出

7181
600
No Answers

---

题目分析

x^2-n*y^2=1：佩尔方程
n为完全平方数时方程无解
暴力求出最小整数解（当然也可以用麻烦的连分数解法）然后矩阵快速幂

$$\left[ \begin{matrix} x_k \\ y_k \end{matrix} \right] = \left[ \begin{matrix} x_1 & dy_1 \\ y_1 & x1 \end{matrix} \right]^{k-1} \left[ \begin{matrix} x_1 \\ y_1 \end{matrix} \right]$$

---

源代码

#include<algorithm>
#include<iostream>
#include<iomanip>
#include<cstring>
#include<cstdlib>
#include<vector>
#include<cstdio>
#include<cmath>
#include<queue>
using namespace std;
typedef long long LL;
inline const int Get_Int() {
    int num=0,bj=1;
    char x=getchar();
    while(x<'0'||x>'9') {
        if(x=='-')bj=-1;
        x=getchar();
    }
    while(x>='0'&&x<='9') {
        num=num*10+x-'0';
        x=getchar();
    }
    return num*bj;
}
const int maxn=5,mod=8191;
struct Matrix {
    LL n,m,a[maxn][maxn];
    Matrix(LL n,LL m) {
        this->n=n;
        this->m=m;
        memset(a,0,sizeof(a));
    }
    Matrix(LL n,LL m,char E) { //单位矩阵
        this->n=n;
        this->m=m;
        memset(a,0,sizeof(a));
        for(int i=1; i<=n; i++)a[i][i]=1;
    }
    void init(LL n,LL m) {
        this->n=n;
        this->m=m;
        memset(a,0,sizeof(a));
    }
    LL* operator [] (const LL x) {
        return a[x];
    }
    Matrix operator = (Matrix b) {
        init(b.n,b.m);
        for(int i=1; i<=b.n; i++)
            for(int j=1; j<=b.m; j++)
                a[i][j]=b[i][j];
    }
    Matrix operator + (Matrix& b) {
        for(int i=1; i<=b.n; i++)
            for(int j=1; j<=b.m; j++)
                a[i][j]+=b[i][j];
    }
    void operator += (Matrix& b) {
        *this=*this+b;
    }
    Matrix operator - (Matrix& b) {
        for(int i=1; i<=b.n; i++)
            for(int j=1; j<=b.m; j++)
                a[i][j]-=b[i][j];
    }
    void operator -= (Matrix& b) {
        *this=*this-b;
    }
    Matrix operator * (LL k) {
        for(int i=1; i<=n; i++)
            for(int j=1; j<=m; j++)
                a[i][j]*=k;
    }
    void operator *= (LL k) {
        *this=*this*k;
    }
    Matrix operator * (Matrix& b) {
        Matrix c(n,b.m);
        for(int i=1; i<=n; i++)
            for(int j=1; j<=b.m; j++)
                for(int k=1; k<=m; k++)
                    c[i][j]=(c[i][j]+a[i][k]*b[k][j])%mod;
        return c;
    }
    void operator *= (Matrix& b) {
        *this=*this*b;
    }
    Matrix operator ^ (LL b) {
        Matrix ans(n,m,'e'),a=*this;
        while(b>0) {
            if(b&1)ans=ans*a;
            a*=a;
            b>>=1;
        }
        return ans;
    }
    void Output() {
        for(int i=1; i<=n; i++) {
            for(int j=1; j<=m; j++)
                printf("%lld ",a[i][j]);
            putchar('\n');
        }
    }
};
LL x,y;
void Get_Answer(LL d) { //找到佩尔方程最小特解 
    y=1;
    while(true) {
        x=sqrt(d*y*y+1);
        if(x*x-d*y*y==1)break;
        y++;
    }
} 
//x^2-n*y^2=1
LL n,k;
int main() {
    while(scanf("%lld%lld",&n,&k)!=EOF) {
        int t=sqrt(n);
        if(t*t==n) {
            puts("No Answers");
            continue;
        }
        Get_Answer(n);
        Matrix a(2,2);
        a.a[1][1]=x%mod;
        a.a[1][2]=y*n%mod;
        a.a[2][1]=y%mod;
        a.a[2][2]=x%mod;
        a=a^(k-1);
        x=(a.a[1][1]*x%mod+a.a[1][2]*y%mod)%mod;
        printf("%lld\n",x);
    }
    return 0;
}

---