ペル方程式

ペル方程式とは，ディオファントス方程式の一種で，

X**2 = d * Y**2 + 1

の形で表される．任意のdに対し，(X,Y)=(1,0)は自明な解であるが，
dが完全平方ではない自然数であれば，必ず自明でない自然数の組の解(X,Y)が存在する．
解は，可算個存在するが，そのうち，Xの値がもっとも小さいものを基本解といい，
dの平方根の連分数近似から，求める事ができる．
(式から直ちに，X/Yは，dの平方根の近似値であることがわかる)

****************************************************************************
pell.sas
Pell方程式
2008.7.2 翔
****************************************************************************;

options nocenter;

/*
連分数近似を使って，ペル方程式
 X**2 = d * Y**2 + 1 (X,Y自然数，d完全平方でない自然数)
を解く
連分数近似
与えられた数値nから整数部を取り出し，
残った小数部の逆数から再び整数部を取り出し．．を繰り返す
得られた整数の列をn0,n1,n2,..とすると
n=n0+1/(n1+1/(n2+..))
と連分数が得られる．
これを適当なところで打ち切り，nの近似値とする
各段階での近似分数の分子分母(Xn/Yn)は，
X(-2)=0,X(-1)=1,X(n+2)=n2X(n+1)+X(n)
Y(-2)=1,Y(-1)=0,Y(n+2)=n2Y(n+1)+Y(n)
*/


data pell;
  d=97; /*d=61 桁落ちのため，解が求まらない*/
  n=sqrt(d);

  x0=0; /* x(-2) */
  x1=1; /* x(-1) */
  y0=1; /* y(-2) */
  y1=0; /* y(-1) */

  do i=0 to 100;

  n2=int(n);
  x2=n2*x1+x0;
  y2=n2*y1+y0;

  solve=x2**2-d*y2**2; /* x**2 - d*y**2 =1 ? */

  put  d=  "X=" x2 "Y=" y2 solve=;
  if solve=1 then leave;

  n=1/(n-n2);
  x0=x1;
  x1=x2;
  y0=y1;
  y1=y2;

  end;
run;
/*
d=97 X=9 Y=1 solve=-16
d=97 X=10 Y=1 solve=3
d=97 X=59 Y=6 solve=-11
d=97 X=69 Y=7 solve=8
d=97 X=128 Y=13 solve=-9
d=97 X=197 Y=20 solve=9
d=97 X=325 Y=33 solve=-8
d=97 X=522 Y=53 solve=11
d=97 X=847 Y=86 solve=-3
d=97 X=4757 Y=483 solve=16
d=97 X=5604 Y=569 solve=-1
d=97 X=105629 Y=10725 solve=16
d=97 X=111233 Y=11294 solve=-3
d=97 X=661794 Y=67195 solve=11
d=97 X=773027 Y=78489 solve=-8
d=97 X=1434821 Y=145684 solve=9
d=97 X=2207848 Y=224173 solve=-9
d=97 X=3642669 Y=369857 solve=8
d=97 X=5850517 Y=594030 solve=-11
d=97 X=9493186 Y=963887 solve=3
d=97 X=53316447 Y=5413465 solve=-16
d=97 X=62809633 Y=6377352 solve=1 <-基本解
 */