PEをMathematicaで

Project Eulerに挑戦してみよう

Project Euler 233 / 円上の格子点

(0,0), (N,0), (0,N), (N,N) を通る円周上の格子点の個数を f(N) とする。

f(10000) = 36

f(N) = 420 をみたす正の整数 N (≤ 10^11)すべての合計を求めよ。

Problem 233 - Project Euler


円の方程式はこうなります。

(x-n/2)^2+(y-n/2)^2=(n/2)^2

 \therefore 2n^2=(2x-n)^2+(2y-n)^2

これをみたす整数(x, y)が420個ある条件は,2つの平方数の和で 2n^2 を表す方法が420通りあることです。4で割った余りが1の素因数に注目します。

4で割った余りが1の素数を pi,余り3の素数を qi として

 n=2^a \displaystyle\prod_i p_i{}^{b_i} \, \prod_j q_j{}^{c_j}

と素因数分解できたとします。この n に対して 2n^2 を平方数の和として表す方法の数は

 4(2b_1+1)(2b_2+1)\cdot\cdots\cdot(2b_m+1)\ \text{通り}

これが420になることから,n の形が絞り込めます。

 n=xy^2z^3,\, y^3 z^7,\, y^2 z^{10}

このように4で割った余りが1の素数の扱いはさほど難しくないのですが,n の素因数分解における2とqiの扱いは面倒です。

たとえば4で割った余りが1の素数だけに注目して m=5*13^2*17^3 をみつけたとして,これに2や q1=3 などをかけたものも条件をみたします。その和を求めるには,1以上 10^11/m 以下の範囲から4で割った余りが1の素因数をもたない数を探さなければなりません。これは1以上10^9以下の範囲で100以下の素因数をもたない数を探す問題204に近いです。

variee.hatenadiary.com

前に書いたコードを流用してもよかったのですが,組み込み関数の SquaresR を使いました。4で割った余りが1の素数だけに注目して求めた数を m として,(mk)^2 を平方数の和で表す方法が420通りあるような k を求めています。

ループが多かったり最後の処理が適当だったりする割に時間はかからず,約44秒で終わりました。