佩尔方程(超详细推导+例题讲解) 每日一遍，算法再见！

佩尔方程

第一类佩尔方程

定义：形如$x^2-d{y}^2=1$(d>1,且d不是完全平方数)
要求第一类佩尔方程的解都是正整数解,也即$(x,y),x>0,y>0$

为什么要求d不是完全平方数呢？
我们假设d是完全平方数
则$x^2-d{y}^2=1$等价于$(x+\sqrt{d}\ast y)(x-\sqrt{d}\ast y)=1$
所以它的解只有可能是(1,0),(-1,0)，不是正整数解。

为什么d>1呢?因为如果d<=1，解也不会是正整数解，自己证明试试看。

显然佩尔方程$x^2-d{y}^2=1$有无数多个解，那么如何求最小的正整数解呢

我们把式子化简一下就变成了，$x=\sqrt{1+dy2}$,我们从y=1开始枚举y，如果$\sqrt{1+dy2}$是整数的话，此时的x,y就是最小的正整数，也就得到了最小正整数解

那么如何得到所有解呢，也就是说我们最小正整数解是(x1,y1)，那么如何得到(x2,y2),(x3,y3),…？

这里给出通解的迭代公式
$x_n=x_{n-1}{x}_1+d{y}_{n-1}{y}_1$
$y_n=x_{n-1}{y}_1+y_{n-1}{x}_1$
这个通解公式怎么推导的呢？我们要知其然，而知其所以然

推导：我们令(x1,y1),(x2,y2)为方程的两个特解
$x_1^2-d{y}_1^2=1$----①
$x_2^2-d{y}_2^2=1$----②
两式想乘得到
$(x_1^2-d{y}_1^2)\ast (x_2^2-d{y}_2^2)=1$
展开得到
$x_1^2{x}_2^2-d{x}_1^2{y}_2^2-d{y}_1^2{x}_2^2+d^2{y}_1^2{y}_2^2=1$
-> $[x_1^2{x}_2^2+d^2{y}_1^2{y}_2^2]-d[d{x}_1^2{y}_2^2+y_1^2{x}_2^2]=1$
-> $[x_1^2{x}_2^2+d^2{y}_1^2{y}_2^2+2d{x}_1{x}_2{y}_1{y}_2]-d[d{x}_1^2{y}_2^2+y_1^2{x}_2^2+2x_1{x}_2{y}_1{y}_2]=1$
-> $(x_1{x}_2-d{y}_1{y}_2{)}^2-d(x_1{y}_2+x_2{y}_1{)}^2=1$
所以$x_3=x_2{x}_1+d{y}_2{y}_1,y_3=x_2{y}_1+x_1{y}_2$
我们把下标换为n就得到
$x_n=x_{n-1}{x}_1+d{y}_{n-1}{y}_1$
$y_n=x_{n-1}{y}_1+y_{n-1}{x}_1$

这样的话我们只要知道第一个特解$(x1,y1)$，和第n-1个特解$(x_{n-1},y_{n-1})$,就能够求出第n个解，但是很明现如果暴力求第n个解的话，时间复杂度是O(n),如果n很大的话就会超时，我们这时候可以用矩阵快速幂logn的复杂度来就快速求取，第n个解。

步骤：

step1:先通过暴力找到最小正整数解，也就是第一组解(x1,y1)，也有一种快速找最小正整数解的方法，叫做连分数法，感兴趣的小伙伴自行百度。

代码：

#include<bits/stdc++.h>
typedef long long ll;
void solve(ll &x,ll &y,ll d)
{
     y=1;
     for(y=1;;y++)
     {
      tempc = ceil(sqrt(d*y*y+1));
      tempf = floor(sqrt(d*y*y+1));
     if(tempc==tempf) break;//只有整数的ceil和floor相同，也就是向上取整和向下取整
     }
	 x=tempc;
	 return ; 
}

step2:根据迭代公式找到所有我们要找的解，也就是
$x_n=x_{n-1}{x}_1+d{y}_{n-1}{y}_1$
$y_n=x_{n-1}{y}_1+y_{n-1}{x}_1$

方法1:暴力迭代法
代码：

#include<bits/stdc++.h>
typedef long long ll;
using namespace std;
void solve(ll &x,ll &y,ll d)
{
     ll tempc,tempf; 
     for(y=1;;y++)
     {
      tempc = ceil(sqrt(d*y*y+1));
      tempf = floor(sqrt(d*y*y+1));
     if(tempc==tempf) break;
     }
	 x=tempc;
	 return ; 
}
int main()
{
	ll x1,y1,d;
	cin>>d;//输入佩恩方程的d 
	solve(x1,y1,d);
	ll xxi=x1,yyi=y1;
	ll xi,yi;
	int n; 
	for(int i=2;i<=n;i++)//这里可以用于暴力求第n个解输入n就行 
	{
	   xi = xxi*x1+d*yyi*y1;
	   yi = xxi*y1+yyi*x1;
	   xxi = xi;
	   yyi = yi; 
	}
	cout<<xi<<' '<<yi<<endl;
	return 0;
}

方法2:矩阵快速幂迭代法
如果我们需要用矩阵快速幂来求的话，我们就需要把通解转化为矩阵的表达形式。

$\left[\begin{array}{c}{x}_n\\ y_n\end{array}\right]$=${\left[\begin{array}{cc}{x}_1& d{y}_1\\ y_1& x_1\end{array}\right]}^{n-1}$$\left[\begin{array}{c}{x}_1\\ y_1\end{array}\right]$
这样我们就可以用矩阵快速幂来求取方程的第n个解，时间复杂度是O(logn),矩阵快速幂就类似于普通的快速幂，只不过我们的底数是一个矩阵，而不是一个常数，如果不了解矩阵快速幂的，请出门右转(矩阵快速幂)。

代码：

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll Maxn=10;
/*矩阵乘法*/ 
void Matrix_ksm(ll a[][Maxn],ll b[][Maxn],ll size)
{

	ll temp[Maxn][Maxn]={0};//注意一定要初始化为0 
	for(ll i=1;i<=size;i++)
	 for(ll j=1;j<=size;j++)
	  for(ll k=1;k<=size;k++)
	  {
	  	  temp[i][j]+=a[i][k]*b[k][j];
	  }
	  for(ll i=1;i<=size;i++)
	   for(ll j=1;j<=size;j++)
	   a[i][j]=temp[i][j];
	   
	  return ;
}
void solve(ll &x,ll &y,ll d)
{
     ll tempc,tempf; 
     for(y=1;;y++)
     {
      tempc = ceil(sqrt(d*y*y+1));
      tempf = floor(sqrt(d*y*y+1));
     if(tempc==tempf) break;
     }
	 x=tempc;
	 return ; 
}
int main()
{
	ll x1,y1;
	ll n,d;
	printf("请输入d:\n");
	cin>>d;
	solve(x1,y1,d);//我们先求第一个解，也就是最小正整数解 
//	cout<<x1<<' '<<y1<<'\n';
	
	printf("请输入要求的第n个矩阵,n=:\n");
	cin>>n;
	ll ans[Maxn][Maxn]={0};
    for(ll i=1;i<=Maxn;i++) ans[i][i]=1;
    ll dishu[Maxn][Maxn];
    dishu[1][1]=x1;
    dishu[1][2]=d*y1;
    dishu[2][1]=y1;
    dishu[2][2]=x1;
    ll zhishu=n-1;
    /*矩阵快速幂*/ 
    while(zhishu)
    {
    	if(zhishu&1)
    	{
    		zhishu--;
    		Matrix_ksm(ans,dishu,2);
		}
		else 
		{
			zhishu>>=1;
			Matrix_ksm(dishu,dishu,2);
		}
	}
    ll xn = ans[1][1]*x1+ans[1][2]*y1;
	ll yn = ans[2][1]*x1+ans[2][2]*y1;
	cout<<xn<<' '<<yn<<endl;
	return 0; 
}

第一类佩尔方程例题讲解

有句熟话说得好“光说不练假把式”，下面我们就来分析几道例题

Street Numbers
题目大意:找出n，m使得$1+2+3+....+(n-1)=(n+1)+(n+2)....+m$,输出前十个解(n,m)。每个数字占是个字符，向右对齐，且每一排升序排列。

分析：我们不难发现，等式的左右都是一个等差数列，那么我们处理一下就可以得到这样一个等式，$(2m+1{)}^2-8n^2=1$,我们令$x=(2m+1),y=n,d=8$,这不就是一个裸的佩尔方程吗？我们只需要找到x，y就能找到m，n。

AC代码：

#include<iostream>
#include<cmath>
#include<cstdio>
using namespace std;
typedef long long ll;
/*找佩尔方程的第一个解，也就是正整数解*/
void solve(ll &x,ll &y,ll d)
{
     ll tempc,tempf; 
     for(y=1;;y++)
     {
      tempc = ceil(sqrt(d*y*y+1));
      tempf = floor(sqrt(d*y*y+1));
     if(tempc==tempf) break;
     }
	 x=tempc;
	 return ; 
}
int main()
{
	ll x1,y1,d=8;
	solve(x1,y1,d);
	ll xxi=x1,yyi=y1;
	ll xi,yi;
	/*通过迭代公式依次找到前十个解，注意这道题答案，第一个解不算在里面
	，所以前十个解是从第2到第11个解*/
	for(int i=2;i<=11;i++)//这里可以用于暴力求第n个解输入n就行 
	{
	   xi = xxi*x1+d*yyi*y1;
	   yi = xxi*y1+yyi*x1;
	   xxi = xi;
	   yyi = yi; 
	   printf("%10lld%10lld\n",yyi,(xxi-1)/2);/*注意题目要求每一个数字占
	   十个字符，每一排升序排列，所以我们先输出n，
	   再输出m，因为(n<m),*/
	}
	return 0;
} 

我们再来一道板子题No more tricks, Mr Nanguo

题目大意：就是一开始乐队有n个人，可以站成一个边长为x的正方形，现在走了一个人，把剩余的n-1个人分成d(这个d题目会给出)个边长为y的正方形站列，显然n有无穷多个解，对应的x，y也有无穷多组解，现在让你求第k(k题目会给出)小的正整数解x

这道题自己尝试着写一下，我就不写题解了，不过提醒一点，运算的时候注意取模，不然数据太大long long都可能超,由于k会比较大，所以不建议用暴力迭代法会超时，用矩阵快速幂来解决。

AC代码：

#include<iostream>
#include<cmath>
#include<cstdio>
using namespace std;
typedef long long ll;
const int mod=8191;
const ll Maxn=10;
/*矩阵乘法*/ 
void Matrix_ksm(ll a[][Maxn],ll b[][Maxn],ll size)
{

	ll temp[Maxn][Maxn]={0};//注意一定要初始化为0 
	for(ll i=1;i<=size;i++)
	 for(ll j=1;j<=size;j++)
	  for(ll k=1;k<=size;k++)
	  {
	  	/*在矩阵运算的时候，我们把每一个可能超范围的地方都取模*/
	  	  temp[i][j]=(temp[i][j]%mod+a[i][k]%mod*b[k][j]%mod)%mod;
	  }
	  for(ll i=1;i<=size;i++)
	   for(ll j=1;j<=size;j++)
	   a[i][j]=temp[i][j];
	   
	  return ;
}
void solve(ll &x,ll &y,ll d)
{
     ll tempc,tempf; 
     for(y=1;;y++)
     {
      tempc = ceil(sqrt(d*y*y+1));
      tempf = floor(sqrt(d*y*y+1));
     if(tempc==tempf) break;
     }
	 x=tempc;
	 return ; 
}
int main()
{
	ll x1,y1;
	ll n,d;
	while(scanf("%lld %lld",&d,&n)!=EOF)
	{
	if(d<=1||ceil(sqrt(d))==floor(sqrt(d)))
	{
		cout<<"No answers can meet such conditions"<<'\n';
		continue;
	} 
	solve(x1,y1,d);//我们先求第一个解，也就是最小正整数解 
	ll ans[Maxn][Maxn]={0};
    for(ll i=1;i<=Maxn;i++) ans[i][i]=1;
    ll dishu[Maxn][Maxn];
    dishu[1][1]=x1;
    dishu[1][2]=d*y1;
    dishu[2][1]=y1;
    dishu[2][2]=x1;
    ll zhishu=n-1;
    /*矩阵快速幂*/ 
    while(zhishu)
    {
    	if(zhishu&1)
    	{
    		zhishu--;
    		Matrix_ksm(ans,dishu,2);
		}
		else 
		{
			zhishu>>=1;
			Matrix_ksm(dishu,dishu,2);
		}
	}
    ll xn = ans[1][1]%mod*x1%mod+ans[1][2]%mod*y1%mod;
	ll yn = ans[2][1]%mod*x1%mod+ans[2][2]%mod*y1%mod;
	cout<<xn%mod<<'\n';
     }
	return 0; 
}

第二类佩尔方程

定义:形如$x^2-d{y}^2=k$的方程叫做第二类佩恩方程
第二类佩尔方程可以看成是第一类佩尔的拓展，也相对复杂一点

解第二类佩恩方程，需要用到第一类方程的解。

推导：
我们假设$x^2-d{y}^2=k-----①$的一个特解为$p，q$
设$x^2-d{y}^2=1-----②$的正整数解为$x1,y1$，①×②式就可以得到下面式子。
->$(p^2-d{q}^2)\ast (x_1^2-d{y}_1^2)=k$
->$(p{x}_1\pm dq{y}_1{)}^2-d(p{y}_1\pm q{x}_1{)}^2=k$
所以通解为,$x=p{x}_1\pm dq{y}_1$,$y=p{y}_1\pm q{x}_1$
我们假设②式中d=2，那么$x^2-2y^2=1$，解得最小正整数解是$x_1=3,y_1=2$,那么①式的解就是$x_n=3x_{n-1}\pm 4y_{n-1}$,$y_n=2x_{n-1}\pm 3y_{n-1}$。

遇到第二类佩尔方程也不用害怕，只需要和第一类佩尔方程联系起来求出通解即可。我们更具第二类佩尔方程的通解，
$x=p{x}_1+dq{y}_1$,$y=p{y}_1+q{x}_1$，可以转化为矩阵乘法的形式
$\left[\begin{array}{c}{x}_n\\ y_n\end{array}\right]$=${\left[\begin{array}{cc}p& dq\\ q& p\end{array}\right]}^{n-1}$$\left[\begin{array}{c}{x}_1\\ y_1\end{array}\right]$,这样就方便用矩阵快速幂求第k个解

请关注我看更多数论分支知识体系讲解，也别忘点个赞额！