Project Euler Problem 66

Problem 66

Consider quadratic Diophantine equations of the form:

x² – Dy² = 1

For example, when D=13, the minimal solution in x is 649² – 13×180² = 1.

It can be assumed that there are no solutions in positive integers when D is square.

By finding minimal solutions in x for D = {2, 3, 5, 6, 7}, we obtain the following:

3² – 2×2² = 1
2² – 3×1² = 1
9² – 5×4² = 1
5² – 6×2² = 1
8² – 7×3² = 1

Hence, by considering minimal solutions in x for D ≤ 7, the largest x is obtained when D=5.

Find the value of D ≤ 1000 in minimal solutions of x for which the largest value of x is obtained.

    这道题事实上是让我们寻找佩尔方程的特解。最naive的思路自然是x从1开始枚举，固定x，再从小到大枚举y，然而随着D的增大，特解中x的增长速度也是惊人的。（D=991时的最小特解y=12055735790331359447442538767，这时x=379516400906811930638014896080）这不仅需要写一个高精度乘法，而且枚举次数也是不现实的。

    好在数学家给我们提供一种更好的方法：若根号N的连分数展开为{a0；<a1,a2,...,an,2a0>}，记p/q=[a0,a1,...an],则这个佩尔方程的特解就是(p,q)(D为偶数),(2p^2+1,2pq)(D为奇数时)。（关于上述符号的意义可以参见Project Euler64&65,定理的证明可以参见二潘的《初等数论》，里面用连分数理论刻画了Pell方程）感觉这道题放在这里别有深意，因为Project Euler Problem 65和Problem 64恰恰是关于连分数计算和二次根式的连分数展开的。

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <malloc.h>

typedef struct {
	int fenzi;
	int fenmu;
	int k;
	int irrationalPart;
}twoReal;  // ÓÐÀíÊý+¶þ´Î¸ùʽ 

int gcd(int a,int b) {
	int da=abs(a),db=abs(b),t=da%db;
	while (t!=0) {
		da = db;
		db = t;
		t = da%db;
	}
	return db;
}

double lengthOfReal(twoReal x )
{
	double ret;
	ret = x.fenzi*1.0/x.fenmu+sqrt(x.irrationalPart)/(x.k);
	return ret;
}

twoReal jieduan(twoReal x) 
{
	twoReal ret;
	int daxiao = (int)lengthOfReal(x);
	ret.fenmu = x.fenmu;
	ret.fenzi = x.fenzi-daxiao*x.fenmu;
	ret.irrationalPart = x.irrationalPart;
	ret.k = x.k;
	return ret;
}

twoReal daoshu(twoReal x) 
{
	twoReal ret;
	int gongyue=1;
	ret.fenzi = x.k*x.k*x.fenmu*x.fenzi;
	ret.fenmu = x.k*x.k*x.fenzi*x.fenzi-x.fenmu*x.fenmu*x.irrationalPart;
	ret.k = (-x.k*x.k*x.fenzi*x.fenzi+x.fenmu*x.fenmu*x.irrationalPart)/(x.k*x.fenmu*x.fenmu);
	ret.irrationalPart = x.irrationalPart;
	gongyue = gcd(ret.fenzi,ret.fenmu);
	ret.fenzi /= gongyue;
	ret.fenmu /= gongyue;
	if (ret.fenmu<0) {
		ret.fenmu*=-1;
		ret.fenzi*=-1;
	}
	return ret;
}

int xiangdeng(twoReal a,twoReal b) {
	int ret=0;
	if (a.fenmu==b.fenmu&&a.fenzi==b.fenzi&&a.irrationalPart==b.irrationalPart&&a.k==b.k) {
		ret = 1;
	}
	return ret;
}

int chazhao(twoReal a,twoReal loop[],int length) 
{
	int ret=0;
	for (int i=0;i<length;i++) {
		if (xiangdeng(a,loop[i])) {
			ret = 1;
			break;
		}
	}
	return ret;
}

void beijia(int fenzi[],int fenmu[],int a) 
{
	int temp[500] = {0};
	for (int i=0;i<500;i++) {
		temp[i] = a*fenmu[i];
	}
	for (int i=0;i<500;i++) {
		if (temp[i]>9&&i<=498) {
			temp[i+1] += temp[i]/10;
			temp[i] %= 10;
		}
		else if (temp[i]>9) {
			printf("yichu\n");
			break;
		}
	}
	for (int i=0;i<500;i++) {
		fenzi[i]+=temp[i];
		if (fenzi[i]>9) {
			if (i==499) {
				printf("yichu");
				break;
			}
			else {
				fenzi[i+1] += fenzi[i]/10;
				fenzi[i] %= 10;
			}
		}
	}
}

void swap(int fenzi[],int fenmu[]) 
{
	int temp[500];
	for (int i=0;i<500;i++) {
		temp[i] = fenzi[i];
		fenzi[i] = fenmu[i];
		fenmu[i] = temp[i];
	}
}

int main () 
{
	twoReal loop[10000],a;
	memset(loop,0,10000*sizeof(twoReal));
	int j=2,length,cnt=0,lianfenshu[500]={0},lengthOfLoop=0,fenzia[500]={1},fenmua[500]={0},jie,mowei,max=0,weishu;
	for (int i=2;i<=1000;i++) {
		if (i==j*j) {
			j++;
			continue;
		}
		
		memset(fenzia,0,500*sizeof(int));
		memset(fenmua,0,500*sizeof(int));
		memset(loop,0,10000*sizeof(twoReal));
		memset(lianfenshu,0,500*sizeof(int));
		fenzia[0]=1;		
		a.fenmu = 1;
		a.fenzi = 0;
		a.irrationalPart = i;
		a.k = 1;
		length = 0;
		lengthOfLoop = 0;
		
		lianfenshu[lengthOfLoop++] = (int)lengthOfReal(a);
		a=jieduan(a);
		loop[length++]=a;
		while (chazhao(a,loop,length-1)==0) {
			a=daoshu(a);
			lianfenshu[lengthOfLoop++] = (int)lengthOfReal(a);
			a=jieduan(a); 
			loop[length++] = a;
		}
		
		for (int k=lengthOfLoop-2;k>=0;k--) {
			swap(fenzia,fenmua);
			beijia(fenzia,fenmua,lianfenshu[k]);
		}
		for (int k=499;k>=0;k--) {
			if (fenzia[k]) {
				weishu=k+1;
				break;
			}
		}
		
		if (!(lengthOfLoop%2)) {
			if (fenzia[weishu-1]>=7) {
				jie=2*weishu+2;
			} else if (fenzia[weishu-1]>=2)  {
				jie=2*weishu+1;
			} else jie=2*weishu;
		} else {
			jie=weishu;
		}   //¼òµ¥¹À¼Æ½âµÄλÊý 
		
		if (jie>=max) {
			printf("%d\t%d\n",jie,i);
			max=jie;
		}
		
//		if (i<=23) printf("%d\t%d\n",i,length-1);
//		if (length%2==0) cnt++;
//		if (length>=10000) printf("yichu");
	}
	printf("%d",cnt);
}

    代码贴下（这道题是用Problem64&65的代码拼起来的，所以可能不那么美观）：