连分数法解佩尔方程特解

连分数法解佩尔方程特解

一、佩尔方程的形式:

二、关于佩尔方程的特解:

         特解是指佩尔方程的最小整数解,容易发现当x最小的时候y也同样达到最小。在一般情况下,佩尔方程的特解是通过暴利枚举方法得到的,本文将介绍如何用应用连分数法求特解。

         本文将不涉及证明,只介绍方法。

三、连分数法:

         一个实数的简单连分数表示,是指将一个实数用以下方法表示:


可以把连分数简记为:

有理数的连分数有两种表示形式:

 

所有无限连分数都是无理数,而所有无理数都可以用一种精确的方式表示成无限连分数,可以用这种方法逼近,无理数的值。

四、关于一个非完全平方数的平方根的连分数表示:

可以证明:一个非完全平方数的平方根的连分数是以周期呈现的。

比如:


简写为:


在之后就会循环出现1,2,4,2,1,8

我们不妨这样记这种连分数的形式:


显然循环节的长度是6

并且还有个重要的特点:这个循环节一定是从开始的,且最后一个数一定是2倍的

五、求解佩尔方程的最小特解:

我们将写成连分数的形式:

并且我们记:

(关于计算p,q:只要按照连分数的展开形式,迭代计算即可)

其中如果记循环节长度为s

那么有如下结论:

1、如果s为偶数时。最小特解为:

2、如果s为奇数时,最小特解为:

六、计算

我们希望得到准确的连分数展开,那么关键在于不用浮点型计算。接下来以为例,解释如何计算的连分数。

我们记当前展开为,那么首先





按照这种方式,我们计算出了的连分数:

然后可以计算出来:

由于循环节长度6是偶数,那么佩尔方程的最小特解是:

之后我们参照上面的例子,来设计计算连分数的算法:

我们记:

那么显然有:

之后我们可以得到:

可以证明,这里一定是大于0的,这个实际上就是下次的

继续推导有:

可以证明,分母是可以被整除的。那么上式就可以写成:

那么容易得到新的b,c是:

还有,结果很大1000以内好多结果都超long long了。。。要改成大数才行。。。

七、关于如何解佩尔方程:

         这个请参考AekdyCoin牛的空间:http://hi.baidu.com/aekdycoin/item/a45f7c37850e5b9db80c03d1

八、代码

#include<cstdio>
#include<cstring>
#include<cmath>
#include<cstdlib>
using namespace std;
typedef long long ll;
ll a[20000];
bool pell_minimum_solution(ll n,ll &x0,ll &y0){
	ll m=(ll)sqrt((double)n);
	double sq=sqrt(n);
	int i=0;
	if(m*m==n)return false;//当n是完全平方数则佩尔方程无解
	a[i++]=m;
	ll b=m,c=1;
	double tmp;
	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]);
	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;}
	else{x0=2*p*p+1;y0=2*p*q;}
	return true;
}

int main(){
	ll n,x,y;
	while(~scanf("%lld",&n)){
		if(pell_minimum_solution(n,x,y)){
			printf("%lld^2-%lld*%lld^2=1\t",x,n,y);
			printf("%lld-%lld=1\n",x*x,n*y*y);
		}
	}

一个java打表的代码:
import java.io.*;
import java.math.*;
import java.util.*;

public class main {
	static long [] a = new long [1000];
	static BigInteger x;
	static BigInteger y;
	public static void main(String [] args) throws FileNotFoundException{
		FileReader fin = new FileReader ("data.in");
		File fout = new File("data.out");
		PrintStream pw = new PrintStream(fout);
		Scanner cin = new Scanner(fin);
		System.setOut(pw);
		while(cin.hasNext()){
			long n=cin.nextLong();
			if(pell_solution(n)){
				
				System.out.print("\""+x+"\",");
			}else{
				System.out.print("\"no solution\",");
			}
		}
	}
	public static boolean pell_solution(long D){
		double sq=Math.sqrt((double)D);
		long m=(long) Math.floor(sq);
		int i=0;
		if(m*m==D)return false;
		a[i++]=m;
		long b=m,c=1;
		double tmp;
		do{
			c=(D-b*b)/c;
			tmp=(sq+b)/c;
			a[i++]=(long)(Math.floor(tmp));
			b=a[i-1]*c-b;
		}while(a[i-1]!=2*a[0]);
		BigInteger p=new BigInteger("1");
		BigInteger q=new BigInteger("0");
		BigInteger t;
		for(int j=i-2;j>=0;j--){
			t=p;
			p=q.add(p.multiply(BigInteger.valueOf(a[j])));
			q=t;
		}
		if((i-1)%2==0){
			x=p;y=q;
		}else{
			x=BigInteger.valueOf(2).multiply(p).multiply(p).add(BigInteger.ONE);
			y=BigInteger.valueOf(2).multiply(p).multiply(q);
		}
		return true;
	}
}


  • 3
    点赞
  • 17
    收藏
    觉得还不错? 一键收藏
  • 4
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 4
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值