Project Euler 176 / 隣辺を共有する直角三角形

問題概略
n の約数の個数の条件
- n が奇数のとき
- n が偶数のとき
n の素因数分解

問題概略

辺の長さが $(9,12,15)$ , $(12,16,20)$ , $5,12,13)$ , $(12,35,37)$ の 4 つの直角三角形はすべて隣辺（直角に隣接する辺）として長さ 12 の辺をもっている。
長さ 12 の隣辺をもつ整数辺の直角三角形はこれら 4 つだけである。
47547 個の整数辺の直角三角形の隣辺になる最小の整数を求めよ。
https://projecteuler.net/problem=176

解説の pdf も作りました。きれいなレイアウトで読みたい方はこちらをどうぞ。

drive.google.com

n の約数の個数の条件

直角三角形の3辺を $x$ , $y$ , $n$ とおきます。 $y$ が斜辺です。 $x^2+n^2=y^2$ から

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

$n$ に対応する三角形の個数を $f(n)$ とおきます。

$f(n)$ は上の方程式の自然数解の個数で，要するにこの問題で聞いていることは「 $n^2$ を積の形に直す方法は何通りありますか？」です。

$x$ , $y$ の条件は「 $y+x>y-x>0$ 」かつ「 $y\pm x$ の偶奇が一致する」です。
以下， $k$ の約数の個数を $\verb|div_count|(k)$ であらわします。

n が奇数のとき

$n^2$ は奇数なので $y\pm x$ は両方とも奇数です。

$n^2$ の約数のうち $y+x=y-x=n$ の 1 個を除外してから $y+x>y-x$ を考えて 2 で割ります。

$f(n)=\frac{1}{2}\left\{\verb|div_count|(n^2)-1\right\}$

$f(n)=47547$ のとき

$\verb|div_count|(n^2)=47547\cdot 2+1=5\cdot 7\cdot 11 \cdot 13 \cdot 19\mbox{　……(1)}$

n が偶数のとき

$n^2$ は偶数なので $y\pm x$ の偶奇が一致するとき $y\pm x$ は両方とも偶数です。

$y+x=2a$ , $y-x=2b$ とおくと $ab=\frac{n^2}{4}$ になります。 $n$ が奇数の場合と同様に考えると

$f(n)=\frac{1}{2}\left\{\verb|div_count|\left(\frac{n^2}{4}\right)-1\right\}$

$f(n)=47547$ のとき

$\verb|div_count|\left(\frac{n^2}{4}\right)=47547\cdot 2+1=5\cdot 7\cdot 11 \cdot 13 \cdot 19\mbox{　……(2)}$

n の素因数分解

$p_1,\, p_2,\, \cdots$ を素数として $k_1,\, k_2,\, \cdots$ をその指数とします。

$n=p_1{}^{k_1}p_2{}^{k_2}\cdots$ のとき $n^2=p_1{}^{2k_1}p_2{}^{2k_2}\cdots$ の約数の個数は $(2k_1+1)(2k_2+1)\cdots$ 個です。

(1)(2)から次のことがわかります。

$n$ に使う素数は 5 個
$5=2\cdot 2+1$ などより指数は 2, 3, 5, 6, 9

最小の $n$ を求めたいので「素数は 2 または 3 からはじまる連続する 5 個」「小さな素数の指数は大きい」もわかります。
求める $n$ の候補は次の 2 つです。

$n$ が奇数のとき $n=3^9\cdot 5^6\cdot 7^5\cdot 11^3\cdot 13^2$
$n$ が偶数のとき $n=2\cdot(2^9\cdot 3^6\cdot 5^5\cdot 7^3\cdot 11^2)$

最小値は $n$ が偶数のときの値です。

桁数の大きな電卓アプリならそのまま計算できそうな値ですが，project euler の精神にのっとって mathematica で解きました。計算時間は約 0.00005 秒です。

$n$ の指数のリスト powerList を作る
powerList の長さをもとに素数のリスト primesEven, primesOdd を作る
$n$ の候補 2 つのうち小さい方が答え

    powerList = (First@# - 1)/2 & /@ Reverse@FactorInteger[47547*2 + 1];
    len = Length@powerList;
    primesEven = Prime@Range@len;
    primesOdd = Prime@Range[2, len + 1];
    numOdd = Inner[Power, primesOdd, powerList, Times];
    numEven = 2*Inner[Power, primesEven, powerList, Times];
    ans = Min[numOdd, numEven]

砂場で遊ぼう

c++/python/mathematica などの練習帳