Posted at

Project Euler Q66 【ディオファントス方程式】

More than 1 year has passed since last update.

Project Eulerをワンライナーで解いてみる。

間違っていたらコメントください。


問題

次の形式の, $2$次のディオファントス方程式を考えよう:

$x^2 - Dy^2 = 1$

たとえば $D=13$ のとき, $x$ を最小にする解は $649^2 - 13×180^2 = 1$ である.

$D$ が平方数$(square)$のとき, 正整数のなかに解は存在しないと考えられる.

$D = {2, 3, 5, 6, 7}$ に対して $x$ を最小にする解は次のようになる:

$3^2 - 2×2^2 = 1$

$2^2 - 3×1^2 = 1$

$9^2 - 5×4^2 = 1$

$5^2 - 6×2^2 = 1$

$8^2 - 7×3^2 = 1$

したがって, $D ≤ 7$ に対して $x$ を最小にする解を考えると, $D=5$ のとき $x$ は最大である.

$D ≤ 1000$ に対する $x$ を最小にする解で, $x$ が最大になるような $D$ の値を見つけよ.


アプローチ

答えの$D$のとき$x$は$38$桁の数になるようで、ループして計算するやり方だと到底できない。

ペル方程式というものを使用すれば解けるらしいのでやってみる。

(ちなみにペルさんが発見したものではないらしい。)

1.$\sqrt{D}$を連分数展開し、$\sqrt{D}=[a_0;(a_1,a_2,...,a_m)]$とおく。

2.$\sqrt{D}$を元に$[a_0;a_1,a_2,...,a_{m-1}]$を考えると、

これは$\sqrt{D}$の近似分数$\frac{P}{Q}$を表しているが、

実はこの$P$と$Q$が$P^2-DQ^2=±1$を満たす。

3.$m\%2=0$の場合は$P^2-DQ^2=1$であり、$(P,Q)=(x,y)$である。

$m\%2=1$の場合は$P^2-DQ^2=-1$であるがペル方程式の性質により、

$x+y\sqrt{D}=(P+Q\sqrt{D})^2$が成り立つ為、

$x=P^2+Q^2D$、$y=2PQ$となる。

連分数展開についてはこちらを参照(私の記事)。

Project Euler Q64 【奇数周期の平方根】 - Qiita


解答

time seq 1000 |

awk '{print $1,int(sqrt($1))}' |
awk '$2!=sqrt($1)' |
awk '{a=$2;b=a;c=1;for(i=1;i<=999;i++){if(n[$1,a,b,c]!=1){print $1,a;n[$1,a,b,c]=1;c=($1-b^2)/c;a=int((b+$2)/c);b=$2-((b+$2)%c)}else{break}}}' |
awk '{if(k!=$1){printf "\n"$0" "}else{printf $2" "};k=$1}' |
awk -M 'NF{s[-1]=1;b[-1]=0;s[0]=0;b[0]=1;for(i=3;i<=NF-1;i++){s[i-2]=$i*s[i-3]+s[i-4];b[i-2]=$i*b[i-3]+b[i-4]};print $0,$2*b[NF-3]+s[NF-3],b[NF-3]}' |
awk -M 'a=(NF-4)%2==0{print $1,$(NF-1),$NF} !a{print $1,$(NF-1)^2+$NF^2*$1,2*$(NF-1)*$NF}' |
sort -k2,2n |
tail -1 |
awk '$0=$1'
661

real 0m0.624s
user 0m0.528s
sys 0m0.399s


答え合わせ

こちらのサイト様と一致していればOKとした。

http://kingyojima.net/pje/066.html


参考にさせて頂いたサイト様

Project Euler 66 _ ディオファントス方程式 - PEをMathematicaで

ペル方程式 - Wikipedia