2021-12-27

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]

2021-12-25

Project Euler 169 / ある数を2のべき乗の和で表せる方法の数の探査

dp mathematica project euler python

問題概略
2 進法表示の桁数字に注目
実装
- python で解く
- mathematica で解く

問題概略

整数 $n$ を 2 のべき乗の和で表すことを考える。ただし各数は高々 2 回しか使ってはいけないものとする。
この表し方の数を $f(n)$ とする。 $f(0)=1$ と定義する。
例として $n=10$ を考える。
\begin{align*}
10=1+1+8=1+1+4+4=1+1+2+2+4=2+4+4=2+8
\end{align*}
5 通りの異なる表し方があるので $f(10)=5$ である。
$f(10^{25})$ を求めよ。
https://projecteuler.net/problem=169

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

drive.google.com

2 進法表示の桁数字に注目

$n$ の 2 進法表示を考えます。普通は桁数字に 0 と 1 しか使えませんが，この問題では 2 も使えます。

これはたとえば $10$ という部分があったら 1 を右にずらして $10\to 02$ と変更できるということです。右に 1 があるときはこういう操作はできません。

$n$ の偶奇で場合分けして $2^0$ の位の桁数字に注目します。

$n$ が奇数のとき $n=2m+1$ の 2 進法表示の右端は 1 で，ここは変更できません。他の桁数字の表し方は $m$ の表し方と同じだけあります。

$f(2m+1)=f(m)$

$n$ が偶数のとき $n=2m$ の 2 進法表示の右端は 0 です。ここを変更しない表し方が $f(m)$ 通りあります。

この他に $2^1$ の位の桁数字を減らすかわりに $2^0$ の位の桁数字を 2 に変更する表し方もあります。この部分は正直ちょっとわかりにくく，少し実験が必要でした。

$2=2^1=2^0+2^0$ より $f(2)=2$
$4=2^2=2^1+2^1=2^1+2^0+2^0$ より $f(4)=3$

偶数だからといって $2^1$ の位の桁数字は 1 とは限りませんが，この問題の設定では上の位の数字を右にずらしていって $2^1$ の位の桁数字を 1 または 2 にすることができます。これをもう 1 回右にずらすわけです。

この表し方は $(2m-2)\div 2=m-1$ の表し方と同じ $f(m-1)$ 通りあります。

$f(2m)=f(m)+f(m-1)$

以上を $f(n)=\cdots$ の形に直して dp すれば解けます。

\begin{align*}
f(n)=\begin{cases}
f\left(\dfrac{n-1}{2}\right)& (\text{$n$が奇数のとき})\\[11pt]
f\left(\dfrac{n}{2}\right)+f\left(\dfrac{n}{2}-1\right)& (\text{$n$が偶数のとき})
\end{cases}
\end{align*}

実装

python で解く

$10^{25}$ はかなり大きいので python では LRU キャッシュを効かせないと固まったようになります。キャッシュの効果は絶大で，効かせたら瞬時に答えが出ました。

    import time
    from functools import lru_cache
    
    @lru_cache(maxsize=None)
    def dp(n):
        if n <= 1:
            return 1
        if n % 2 == 1:
            return dp((n - 1) // 2)
        else:
            return dp(n // 2) + dp(n // 2 - 1)
    
    if __name__ == '__main__':
        start = time.time()
        print(dp(10 ** 25))
        elapsed_time = time.time() - start
        print("elapsed_time: {0} [sec]".format(elapsed_time))

mathematica で解く

mathematica でも解きました。こちらもメモ化して高速化しています。計算時間は 0.0007 秒でした。

    f[0] = 1; f[1] = 1;
    f[n_] := f[n] = If[OddQ[n], f[(n - 1)/2], f[n/2] + f[n/2 - 1]];
    ans = f[10^25]

2021-12-24

Project Euler 304 / Primonacci

mathematica project euler

問題概略
どう解くか
- 素数は NextPrime で処理
- フィボナッチ数は MatrixPowerMod で処理
実装
- 並列処理
- 前の項を使い回す

問題概略

$n$ より大きい最小の素数を $\verb|next_prime|(n)$ であらわす。
これを用いて $a_n$ を次のように定義する。

$a_1=\verb|next_prime|(10^{14})$

$a_n=\verb|next_prime| (a_{n-1}) \quad (n>1)$

また，次の式をみたすフィボナッチ数 $f(n)$ を用いて $b_n=f(a_n)$ と定義する。

$f(0)=0,\, f(1)=1$

$f(n)=f(n-1)+f(n-2) \quad (n>1)$

$\displaystyle\sum_{n=1}^{10^5} b_n=\sum_{n=1}^{10^5} f(a_n)$ を $\bmod\, 1234567891011$ で求めよ。
https://projecteuler.net/problem=304

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

drive.google.com

どう解くか

この問題のポイントは「大きな素数を求めること」「大きなフィボナッチ数を求めること」です。これらをうまく処理できる言語にとっては，この問題は「やるだけ」です。

mathematica はこの条件をみたしていて，問題文をほぼそのままコードに落とすだけで解けました。

素数は NextPrime で処理

大きな素数には組み込み関数の NextPrime が使えます。例として $a_{100000}$ を求めてみました。

    In[]:= a[n_] := NextPrime[10^14, n];
        a[10^5]
    
    Out[]= 100000003235443

フィボナッチ数は MatrixPowerMod で処理

フィボナッチ数 $f_n$ は行列を使って求めます。
まず漸化式 $f_{n+2}=f_{n+1}+f_n$ を行列のかけ算に直します。

$\begin{pmatrix} f_{n+2}\\f_{n+1} \end{pmatrix} =\begin{pmatrix} f_{n+1}+f_n\\f_{n+1} \end{pmatrix} =\begin{pmatrix} 1&1\\1&0 \end{pmatrix} \begin{pmatrix} f_{n+1}\\f_{n} \end{pmatrix}$

$f_n$ は次のベクトルの第 1 成分です。

$\begin{pmatrix} f_{n}\\f_{n-1} \end{pmatrix} =\begin{pmatrix} 1&1\\1&0 \end{pmatrix}^{n-1} \begin{pmatrix} f_{1}\\f_{0} \end{pmatrix} =\begin{pmatrix} 1&1\\1&0 \end{pmatrix}^{n-1} \begin{pmatrix} 1\\0 \end{pmatrix}$

この計算を組み込み関数の MatrixPowerMod で処理すると，同時に $\bmod$ の計算もできて楽です。この関数は現行のマニュアルには載っていませんが，何の問題もなく使えます。StackExchange を参考にしました。

mathematica.stackexchange.com

$b_n$ の $\bmod$ をとったものを $c_n$ とおいて，次のように書きました。

    modnum = 1234567891011;
    mat = {{1, 1}, {1, 0}}; 
    c[n_] := Mod[First@#, modnum] &@
        (Algebra`MatrixPowerMod[mat, a@n - 1, modnum] . {1, 0});

実装

並列処理

上にあげたコードを使って計算したらちょっと遅かったので並列処理しました。計算時間は約 200 秒です。

    nmax = 10^5;
    a[n_] := NextPrime[10^14, n];
    modnum = 1234567891011;
    mat = {{1, 1}, {1, 0}}; 
    c[n_] := Mod[First@#, modnum] &@
        (Algebra`MatrixPowerMod[mat, a@n - 1, modnum] . {1, 0});
    ans = Mod[Total[ParallelMap[c, Range@nmax]], modnum]

前の項を使い回す

もう少し速くしたかったので， $b_n$ を毎回一から計算するかわりに $b_{n-1}$ を再利用する方法を試してみました。

\begin{align*}
\begin{pmatrix}
b_n\\ \square
\end{pmatrix}
=\begin{pmatrix}
1&1\\1&0
\end{pmatrix}^{a_n-1}
\begin{pmatrix}
1\\0
\end{pmatrix}\quad \therefore
\begin{pmatrix}
b_n\\ \square
\end{pmatrix}
=\begin{pmatrix}
1&1\\1&0
\end{pmatrix}^{a_n-a_{n-1}}
\begin{pmatrix}
b_{n-1}\\ \square
\end{pmatrix}
\end{align*}

計算回数はかなり減りますが，この方法だと並列化できないのでかなり遅くなります。

$n=10^3$ までの和をとる時間を比較しました。元のコードが約 4 秒で，新しく書いたコードは約 12 秒です。

結論として，mathematica では「できるだけ組み込み関数を使って，並列化する」のが最適解だと思います。

    nmax = 10^3;
    a[n_] := a[n] = NextPrime[10^14, n];
    modnum = 1234567891011;
    mat = {{1, 1}, {1, 0}};
    v[1] = Algebra`MatrixPowerMod[mat, a@1 - 1, modnum] . {1, 0};
    v[n_] := v[n] = 
      Mod[#, modnum] &@
        (Algebra`MatrixPowerMod[mat, a[n] - a[n - 1], modnum] . v[n - 1]);
    ans = First@Mod[Total[v /@ Range@nmax], modnum]

砂場で遊ぼう

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

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

問題概略

n の約数の個数の条件

n が奇数のとき

n が偶数のとき

n の素因数分解

Project Euler 169 / ある数を2のべき乗の和で表せる方法の数の探査

問題概略

2 進法表示の桁数字に注目

実装

python で解く

mathematica で解く

Project Euler 304 / Primonacci

問題概略

どう解くか

素数は NextPrime で処理

フィボナッチ数は MatrixPowerMod で処理

実装

並列処理

前の項を使い回す