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