1. 程式人生 > >連分數與佩爾方程特解(最小整數解)

連分數與佩爾方程特解(最小整數解)

在這裡插入圖片描述
在這裡插入圖片描述
在這裡插入圖片描述
在這裡插入圖片描述

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll a[20000];
bool pell_minimum_solution(ll n,ll &x0,ll &y0){//求PELL方程最小整數解
	ll m=(ll)sqrt((double)n);//M是N的平方根向下取整
	if(m*m==n)return false;//當n是完全平方數則佩爾方程無解(不討論正負一,零)
	//下面是把N用連分數形式存,B,C,TMP(即AI)的遞推見解釋以22為例參考
	int i=0;    //連分數的數位
	a[i++]=m;   //A0位整數部分m=4
	ll b=m,c=1; //B=4即整數部分位,C=1即求RN時的分母
	double sq=sqrt(n);//SQ是N的高精度根,相當於r0
	double tmp;//tp在下面的迴圈就是rn
	do{
		c=(n-b*b)/c;
		tmp=(sq+b)/c;
		a[i++]=(ll)(floor(tmp));
		b=a[i-1]*c-b;
		//printf("%lld %lld %lld\n",a[i-1],b,c);
	}while(a[i-1]!=2*a[0]);//當有一位等於整數兩倍就結束
    //下面就是要把連分數形式化成分子分母的形式,求PQ兩個值
	ll p=1,q=0;
	for(int j=i-2;j>=0;j--){
		ll t=p;
		p=q+p*a[j];
		q=t;
		//printf("a[%d]=%lld %lld %lld\n",j,a[j],p,q);
	}
	if((i-1)%2==0){x0=p;y0=q;}//如果I是奇數,X0與Y0都是0
	else{x0=2*p*p+1;y0=2*p*q;}//如果I是偶數,X0是兩倍P方+1,y0是兩倍PQ
	return true;
}
int main(){
	ll n,x,y;
	while(~scanf("%lld",&n)){
		if(pell_minimum_solution(n,x,y)){//輸入N求X,Y的值,有解就輸出
            printf("x=%lld  y=%lld\t",x,y);//x,y
			printf("%lld^2-%lld*%lld^2=1\t",x,n,y);//x^2-n*y^2=1
			printf("%lld-%lld=1\n",x*x,n*y*y);//(x^2)-(n*y^2)=1
		}
	}
    return 0;
}