概要

Project Euler¹は、数学的な問題解決能力を要求する、プログラミング問題集・オンラインプログラミングジャッジである。
一般的なオンラインプログラミングジャッジとは異なり、解答するプログラミング言語は、特に限定されていない。
フォームに入力した値だけで正誤を判別するため、適切なアルゴリズムさえ考案できれば、あらゆるプログラミング言語で解答することが可能である。
そのため、Web上では様々な手法による解答が公開されている。

しかし、シェル芸²による解答は、私が探した限りでは、殆ど見つからなかった。
そこで、私自身の勉強も兼ねて、Project Eulerをシェル芸で解いてみることにした。

本記事では、Project EulerのProblem 14を、シェル芸で解いた際の過程を述べる。

実行環境

Arch Linux (x86_64, 4.3.3-2-ARCH)
bash (4.3.42, POSIX互換モード)
GNU coreutils (8.24)
gawk (4.1.3)
nawk (20121220)
mawk (1.3.4)

問題文

[原文 (https://projecteuler.net/problem=14)]

The following iterative sequence is defined for the set of positive integers:

n \rightarrow n/2 \hspace{3em} (n\ \text{is even}) \\
n \rightarrow 3n+1 \hspace{3em} (n\ \text{is odd})

Using the rule above and starting with 13, we generate the following sequence:

13 \rightarrow 40 \rightarrow 20 \rightarrow 10 \rightarrow 5 \rightarrow 16 \rightarrow 8 \rightarrow 4 \rightarrow 2 \rightarrow 1

It can be seen that this sequence (starting at 13 and finishing at 1) contains 10 terms. Although it has not been proved yet (Collatz Problem), it is thought that all starting numbers finish at 1.
Which starting number, under one million, produces the longest chain?

NOTE: Once the chain starts the terms are allowed to go above one million.

[日本語訳 (http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2014)]

正の整数に以下の式で繰り返し生成する数列を定義する.

n \rightarrow n/2 \hspace{3em} (n\ \text{が偶数}) \\
n \rightarrow 3n+1 \hspace{3em} (n\ \text{が奇数})

13からはじめるとこの数列は以下のようになる.

13 \rightarrow 40 \rightarrow 20 \rightarrow 10 \rightarrow 5 \rightarrow 16 \rightarrow 8 \rightarrow 4 \rightarrow 2 \rightarrow 1

13から1まで10個の項になる. この数列はどのような数字からはじめても最終的には 1 になると考えられているが, まだそのことは証明されていない(コラッツ問題)
さて, 100万未満の数字の中でどの数字からはじめれば最長の数列を生成するか.

注意: 数列の途中で100万以上になってもよい

正答は、注釈³に記載する。

解法

prob_014.sh

1	#!/bin/sh
2	
3	seq 1 999999 \
4	| # 100万未満の数の生成
5	gawk 'ORS=""; {n=$0; while(n!=1){print n " "; n=(n%2) ? 3*n + 1 : n/2}; print n "\n"}' \
6	| # コラッツ数列の出力
7	gawk '{print $1, gensub(/[0-9]* /, "s", "g")}' \
8	| # コラッツ数列の長さを'/s+/'に変換
9	awk '{print $1, length($2)}' \
10	| # コラッツ数列の長さの出力
11	sort -k 2n,2 \
12	| # 数列の長さ順にソート
13	tail -n 1 \
14	| # 数列の長さが最大の値の抽出
15	cut -d ' ' -f 1
16	  # 最長のコラッツ数列を生成する数値の出力

(以下、プロンプトでの実行用コード)

seq 1 999999 | gawk 'ORS=""; {n=$0; while(n!=1){print n " "; n=(n%2) ? 3*n + 1 : n/2}; print n "\n"}' | gawk '{print $1, gensub(/[0-9]* /, "s", "g")}' | awk '{print $1, length($2)}' | sort -k 2n,2 | tail -n 1 | cut -d ' ' -f 1

解説

100万未満の自然数について、それぞれその自然数から始まるコラッツ数列を生成し、最長になるものを抽出する。

まず、100万未満の自然数を生成する。
この処理は、3行目のseqを用いておこなう。

次に、各自然数についてコラッツ数列を生成する。
この処理は、5行目のmawkを用いておこなう。
awkのコードは、問題文にあるコラッツ数列の定義、

n \rightarrow n/2 \hspace{3em} (n\ \text{is even}) \\
n \rightarrow 3n+1 \hspace{3em} (n\ \text{is odd})

を、そのまま実装した。
また、この処理は途中で非常に大きい数値を扱うため、扱える数値の上限が大きいgawkを用いた。

次に、コラッツ数列の長さを求める。
この処理は、7~9行目のgawkとawkを用いておこなう。
まず、長さを求める下準備として、7行目のgawkで、コラッツ数列の各項を単一の文字's'に置換する。
そして、9行目のawkで、コラッツ数列の長さを出力している。
7行目で処理系をgawkと明示した理由は、gawkのみに実装されている組込関数、gensub()を利用するためである。

最後に、最長のコラッツ数列を生成する自然数を求める。
この処理は、11~15行目のsort, tail, cutを用いておこなう。
まず、11行目のsortを用いて、自然数を生成するコラッツ数列の長さ順にソートする。
9行目のawkの処理で、各行の2フィールド目は、コラッツ数列の長さとなっている。
そのため、2フィールド目を数値順にソートすれば良い。
次に、13行目のtailを用いて、数列の長さが最大になるものを抽出する。
ソートによって昇順となっているため、最後尾のものが求める自然数である。
最後に、15行目のtrを用いて、出力を整形する。
これは、ソート時に用いた2フィールド目のコラッツ数列長が不要なためである。

以上が、最長のコラッツ数列を生成する100万未満の数字を求める方法である。

雑記

「こんな別解があるよ」や、「こうすればもっと効率化できるぜ」、「こう書けばエレガントですわよ」等の意見がありましたら、気軽に教えていただけますと幸いです🐺
現在の進捗状況 -> https://github.com/gin135/Project_Euler/tree/master/sh

今回の解法は、処理時間がけっこう長い...

$ cat /proc/cpuinfo | grep 'model info' | uniq
model name  : AMD FX(tm)-8320 Eight-Core Processor

$ time seq 1 999999 | gawk 'ORS=""; {n=$0; while(n!=1){print n " "; n=(n%2) ? 3*n + 1 : n/2}; print n "\n"}' | gawk '{print $1, gensub(/[0-9]* /, "s", "g")}' | awk '{print $1, length($2)}' | sort -k 2n,2 | tail -n 1 | cut -d ' ' -f 1
837799
seq 1 999999  0.05s user 0.03s system 0% cpu 1:39.68 total
gawk   80.96s user 0.45s system 80% cpu 1:40.79 total
gawk '{print $1, gensub(/[0-9]* /, "s", "g")}'  100.38s user 0.53s system 100% cpu 1:40.80 total
awk '{print $1, length($2)}'  3.49s user 0.32s system 3% cpu 1:40.80 total
sort -k 2n,2  1.28s user 0.05s system 1% cpu 1:40.97 total
tail -n 1  0.01s user 0.02s system 0% cpu 1:40.97 total
cut -d ' ' -f 1  0.00s user 0.00s system 0% cpu 1:40.97 total

$ time seq 1 999999 | nawk 'ORS=""; {n=$0; while(n!=1){print n " "; n=(n%2) ? 3*n + 1 : n/2}; print n "\n"}' | gawk '{print $1, gensub(/[0-9]* /, "s", "g")}' | awk '{print $1, length($2)}' | sort -k 2n,2 | tail -n 1 | cut -d ' ' -f 1
837799
seq 1 999999  0.05s user 0.01s system 0% cpu 1:56.86 total
nawk   117.68s user 0.33s system 99% cpu 1:58.08 total
gawk '{print $1, gensub(/[0-9]* /, "s", "g")}'  111.05s user 0.43s system 94% cpu 1:58.08 total
awk '{print $1, length($2)}'  3.66s user 0.37s system 3% cpu 1:58.08 total
sort -k 2n,2  1.35s user 0.04s system 1% cpu 1:58.25 total
tail -n 1  0.04s user 0.00s system 0% cpu 1:58.25 total
cut -d ' ' -f 1  0.00s user 0.00s system 0% cpu 1:58.25 total

$ time seq 1 999999 | mawk 'ORS=""; {n=$0; while(n!=1){print n " "; n=(n%2) ? 3*n + 1 : n/2}; print n "\n"}' | gawk '{print $1, gensub(/[0-9]* /, "s", "g")}' | awk '{print $1, length($2)}' | sort -k 2n,2 | tail -n 1 | cut -d ' ' -f 1
854191
seq 1 999999  0.06s user 0.01s system 0% cpu 1:39.61 total
mawk   67.48s user 0.54s system 67% cpu 1:40.69 total
gawk '{print $1, gensub(/[0-9]* /, "s", "g")}'  100.30s user 0.51s system 100% cpu 1:40.70 total
awk '{print $1, length($2)}'  3.52s user 0.30s system 3% cpu 1:40.69 total
sort -k 2n,2  1.24s user 0.06s system 1% cpu 1:40.87 total
tail -n 1  0.02s user 0.01s system 0% cpu 1:40.87 total
cut -d ' ' -f 1  0.00s user 0.00s system 0% cpu 1:40.87 total

悪あがきでmawkを使ってみたけれど、元々のアルゴリズム自体が遅かった...orz
でも、こう見るとmawkって、AWK処理系の中ではかなり速いんだな。
(mawkだと、途中で演算結果が扱える数値の上限を超えてしまうので、正しい結果が出ない／(^o^)＼)

About - Project Euler (https://projecteuler.net) ↩
シェル芸の定義 >> シェル芸 - 上田ブログ (https://blog.ueda.asia/?page_id=1434) ↩
837799 ↩