HDU - 3292（佩尔方程）

原创于 2020-02-04 15:44:02 发布 · 275 阅读

0 ·

CC 4.0 BY-SA版权

数学专栏收录该内容

20 篇文章

订阅专栏

本文深入探讨了佩尔方程x^2 - ny^2 = 1的求解方法，特别是当n为非平方数时，如何找到最小的正整数解(x1, y1)，并利用矩阵快速幂算法高效计算第k小的解。通过暴力搜索最小特解，再运用矩阵运算，文章提供了一种解决此类问题的有效途径。

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

题意的意思就是给你个N，K;
求第K小的X,满足 $X^2-N*Y^2=1$ （X,Y要为正整数）；
显然这是一个佩尔方程。
定义：若一个不定方程具有这样的形式： $x^2-ny^2=1$
则称此二元二次不定方程为佩尔方程.
在整数域解中
（1） $n$ 为完全平方数时
则原式化为 $x2−(ny)2=1x^2-(\sqrt{n}y)^2=1$
即 $(x−ny)⋅(x+ny)=1(x-\sqrt{n}y)\cdot(x+\sqrt{n}y)=1$
显然
要使等式成立，只有| $(x−ny)∣=∣(x+ny)∣=1(x-\sqrt{n}y)|=|(x+\sqrt{n}y)|=1$
即 $X$ 只有解 $x = 1, - 1, y = 0$
(2) $n$ 为非平方数
若有两组解（ $x_1,y_1$ ）( $x_2,y_2$ )
则 $（x12−ny12）⋅(x22−ny22)=1（x_1^2-ny_1^2）\cdot(x_2^2-ny_2^2)=1$
即 $x_1^2x_2^2+n^2y_1^2y_2^2)-n(x_1^2y_2^2+x_2^2y_1^2)=1$
左右边加上 $nx_1x_2y_1y_2$ 等式仍成立
所以 $x_1x_2+ny_1y_2)^2-n(x_1,y_2+y_1x_2)^2=1$
所以
$x_3=x_1x_2+ny_1y_2$
$y_3=x_1y_2+y_1x_2$
回到本题
$X, Y$ 都必须是正整数，所以我们要暴力求出一组最小特组 $x_1,y_1)$

  ll y=1,x;
  while(1)
  {
      x=sqrt(n*y*y+1);
      if(x*x-n*y*y==1)break;
      y++;
  }

则 $a n s w e r =$ $(x1,ny1y1,x1)k−1{x_1,ny_1\choose y_1,x_1}^{k-1}$ $(x1y1){x_1\choose y_1}$
由矩阵快速幂求即可

#include<bits/stdc++.h>
#define ll long long
#define endl '\n'
#define rep(i,l,r) for(int i=l;i<=r;i++)
#define per(i,r,l) for(int i=r;i>=l;i--)
const int MX=4e2+7;
const int mod=8191;
using namespace std;
int p[MX],k[MX];
const int m=2;
ll qpow(ll a,ll b,ll MOD=mod){for(ll ans=1;;a=a*a%MOD,b>>=1){if(b&1)ans=ans*a%MOD;if(!b)return ans;}}
ll inv(ll a,ll MOD=mod){return qpow(a,MOD-2,MOD);}
ll __gcd(ll a,ll b){return a*b/__gcd(a,b);}
ll a[m][m],b[m][m];
void mul(ll a[m][m],ll b[m][m])
{
    ll c[m][m]={};
    for(int i=0;i<m;i++)
    {
        for(int j=0;j<m;j++)
        {
            for(int k=0;k<m;k++)
                c[i][j]=(c[i][j]+b[i][k]*a[k][j])%mod;
        }
    }
    memcpy(a,c,sizeof c);
}

int main()
{
  ios::sync_with_stdio(0),cin.tie(0);
  ll n,k;
  while(cin>>n>>k){
    int m=sqrt(n+0.5);
  if(m*m==n||k<=0){
    cout<<"No answers can meet such conditions"<<endl;
    continue;
  }
  ll y=1,x;
  while(1)
  {
      x=sqrt(n*y*y+1);
      if(x*x-n*y*y==1)break;
      y++;
  }
  a[0][0]=x;
  a[0][1]=n*y%mod;
  a[1][0]=y;
  a[1][1]=x;
  b[0][0]=x;
  b[1][0]=y;
  k--;
  while(k)
  {
      if(k&1)
        mul(b,a);
      mul(a,a);
      k>>=1;
  }
  cout<<b[0][0]<<endl;
  }
}