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的代码拼起来的，所以可能不那么美观）：