Project Eulerをシェル芸で解いてみる(Problem 26) #シェル芸

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

しかし、シェル芸[^def]による解答は、私が探した限りでは、 @eban さんの記事しか見つからなかった。
そこで、私自身の勉強も兼ねて、Project Eulerをシェル芸で解いてみることにした。
<br>

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

# 実行環境
* Arch Linux (x86\_64, 4.5.1-1-ARCH)
* bash (4.3.42, POSIX互換モード)
* GNU coreutils (8.25)
* gawk (4.1.3)
<br>

* GNU MPFR (3.1.4)

<!---
* bsd-games (2.17)
* curl (7.47.0)
--->

# 問題文
[原文 ([https://projecteuler.net/problem=26](https://projecteuler.net/problem=26))]
#### Reciprocal cycles

> A unit fraction contains $1$ in the numerator. The decimal representation of the unit fractions with denominators $2$ to $10$ are given:
> 
> $1/2	= 	0.5$
> $1/3	= 	0.(3)$
> $1/4	= 	0.25$
> $1/5	= 	0.2$
> $1/6	= 	0.1(6)$
> $1/7	= 	0.(142857)$
> $1/8	= 	0.125$
> $1/9	= 	0.(1)$
> $1/10	= 	0.1$
> 
> Where $0.1(6)$ means $0.166666...$, and has a 1-digit recurring cycle. It can be seen that $1/7$ has a 6-digit recurring cycle.
> 
> Find the value of $d < 1000$ for which $1/d$ contains the longest recurring cycle in its decimal fraction part.

[日本語訳 ([http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2026](http://odz.sakura.ne.jp/projecteuler/index.php?cmd=read&page=Problem%2026))]
#### 逆数の循環節 その1

> 単位分数とは分子が$1$の分数である. 分母が$2$から$10$の単位分数を10進数で表記すると次のようになる.
> 
> $1/2 = 0.5$
> $1/3 = 0.(3)$
> $1/4 = 0.25$
> $1/5 = 0.2$
> $1/6 = 0.1(6)$
> $1/7 = 0.(142857)$
> $1/8 = 0.125$
> $1/9 = 0.(1)$
> $1/10 = 0.1$
> 
> $0.1(6)$は $0.166666...$ という数字であり, 1桁の循環節を持つ. $1/7$ の循環節は6桁ある.
> 
> $d < 1000$ なる $1/d$ の中で小数部の循環節が最も長くなるような $d$ を求めよ.

正答は、注釈[^ans]に記載する。


# 解法
```sh
1	#!/bin/sh
2	
3	seq 1 1000 \
4	| # 数値の列挙
5	factor \
6	| # 素因数分解
7	tr -d ':' \
8	| # 余分な文字の除去
9	awk '{two=0; five=0; for(i=2;i<=NF;i++){$i==2 ? two++ : ($i==5 ? five++ : 0)}; print $1,two,five}' \
10	| # 各数値における約数2, 5の数の集計
11	awk '{print $1,$1/(2^$2)/(5^$3)}' \
12	| # 各数値を2^aおよび5^bで除算
13	gawk -M '$2==1{print $1,0; next}; {l=1; while(10^l % $2 != 1){l++}; print $1,l}' \
14	| # 循環節の長さの算出
15	sort -k 2,2nr \
16	| # 循環節の長さ順にソート
17	awk 'NR==1 && $0=$1'
18	  # 最大値の出力
```

以下、コマンドライン上での実行用コードである。

```sh
seq 1 1000 | factor | tr -d ':' | awk '{two=0; five=0; for(i=2;i<=NF;i++){$i==2 ? two++ : ($i==5 ? five++ : 0)}; print $1,two,five}' | awk '{print $1,$1/(2^$2)/(5^$3)}' | gawk -M '$2==1{print $1,0; next}; {l=1; while(10^l % $2 != 1){l++}; print $1,l}' | sort -k 2,2nr | awk 'NR==1 && $0=$1'
```

## 解説
まず、解答方針を述べる。

有理小数 $1/n$ を循環小数で表したとき、その循環節の長さ $l$ について、次の性質が成り立つことを利用する。

```math
10^l \equiv 1\left( \bmod \frac{n}{2^a 5^b}\right) \hspace{3em} (aとbは、それぞれ2と5の素因数の個数)
```

よって、整数 $1 \sim 1000$ において、各値を素因数分解し、$2$ と $5$ の素因数の個数を集計した後、上記の式を満たす $l$ を求めれば良い。

- - -

次に、ソースコードの解説を述べる。

まず、設問で定義されている整数 $1 \sim 1000$ を生成する。
この処理は、3行目の`seq`を用いておこなう。

次に、$2$ と $5$ の素因数の個数を求めるために、各整数それぞれについて素因数分解をおこなう。
この処理は、5行目の`factor`でおこなう。
また、以降の処理に支障をきたさないように、7行目の`tr`で余分な文字を除去しておく。

そして、素因数 $2$、$5$ の個数の集計をおこなう。
この処理は、9行目の`awk`を用いておこなう。
`factor`によって、2フィールド目以降に素因数が出力されているため、単純にその中から $2$ と $5$ の個数を集計すれば良い。
ただし、最終的には、設問に該当する整数を求める必要がある。
そのため、後の処理を考慮して、1フィールド目に整数を、2・3フィールドにそれぞれ $2$ と $5$ の個数を出力しておく。

その後、前述した性質 $10^l \equiv 1\left( \bmod \frac{n}{2^a 5^b}\right)$ を利用し、$l$ の値を求める。
この処理は、11行目の`awk`と13行目の`gawk`を用いておこなう。
まず、後の演算を簡略化するために、11行目の`awk`で、$\frac{n}{2^a 5^b}$ の計算をおこなう。
この際、1フィールド目の整数は消さずに残しておく。
次に、$l$ の値を算出する。
ここでは、単純に今までの処理で出力した2フィールド目以降の値を、式に当てはめるだけで良い。
ただし、注意点として、$10^l$ を求める処理過程で、`awk`が扱える数値の上限を超えてしまう場合がある。
そのため、多倍長浮動小数点演算ライブラリの、GNU MPFRを利用する。
したがって、この処理では、`gawk`を`-M`オプションを用いて、GNU MPFRを有効にした上で実行する。
(注：GNU MPFRは、`gawk` Ver 4.1.0以上のみ利用可能です。)
以上の処理によって、1フィールド目に整数、2フィールド目に循環節の長さが出力される。

最後に、循環節の長さが最大の整数を求める。
この処理は、15行目の`sort`および17行目の`awk`を用いておこなう。
具体的には、循環節の長さである2フィールド目を基準にして降順数値ソートをおこない、最後に最大値である1行目の1フィールドを出力する。
<br>

以上が、$1/n (nは1から1000までの整数)$を循環小数で表したとき、循環節の長さが最大になる値を求める方法である。


# 参考にしたもの

1. [Python][ProjectEuler]Project Euler Problem 26 循環小数の循環節 - likr’s labo 〈[http://d.hatena.ne.jp/likr/20100125](http://d.hatena.ne.jp/likr/20100125)〉
1. 『循環小数の秘密』(※Internet Archive) 〈[https://web.archive.org/web/20100106234234/http://web2.incl.ne.jp/yaoki/jyunnkan.htm](https://web.archive.org/web/20100106234234/http://web2.incl.ne.jp/yaoki/jyunnkan.htm)〉

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

* @eban さんが、ご自身の日記にて別解を投稿していらっしゃいます。
    - "Project Euler Problem 26 #シェル芸" [http://jarp.does.notwork.org/diary/201602a.html#201602041](http://jarp.does.notwork.org/diary/201602a.html#201602041)
    - Rubyの正規表現を用いた方法と、剰余を利用した正攻法な方法で解いていらっしゃいます。
<br>
<br>

* (私事がゴタゴタしていたので、すっかり放置してしまったけれど... 第22回シェル芸勉強会が近いので、リハビリも兼ねてマイペースに再開をば。)
    - **第22回シェル芸勉強会は、2016年4月30日(土)開催です！**
    - "jus共催、第4回初心者向け詐称疑惑午前のシェル勉強会/第22回ゴールデンウィークの存在疑惑シェル芸勉強会 - USP友の会 | Doorkeeper" [https://usptomo.doorkeeper.jp/events/41629](https://usptomo.doorkeeper.jp/events/41629)
    - 当日はライブ配信もしてくださるそうです。興味のある方はぜひ！
<br>
<br>

* 初めは、@eban さんのように `bc` で小数部を出力した後、循環節を求めようと考えたのだけれど...
    - フロイドの循環検出法を用いたところ、上手くできなかった...
    - 僕の知識不足だけれど、フロイドの循環検出法は最小周期を求められないそう。
<blockquote class="twitter-tweet" data-conversation="none" data-lang="en"><p lang="ja" dir="ltr"><a href="https://twitter.com/gin_135">@gin_135</a> フロイドもそうだけど、フェルマーの小定理を用いるのはいいけど、あれは最小周期じゃないから。周期の上限を見積もることはできるけど。その約数で考えていけばよい。一般には、離散対数問題になるし。</p>&mdash; tmnghryk (@tmnghryk) <a href="https://twitter.com/tmnghryk/status/724977146892951552">April 26, 2016</a></blockquote>
<script async src="//platform.twitter.com/widgets.js" charset="utf-8"></script>
    - 他にエレガントな解法を思いつかなかったので、ググってカンニングしちゃった(｀；ω；´)
* @eban さんの日記でも言及されているけれど、有理分数 $a/b$ を循環小数で表したとき、循環節の長さは、必ず $b-1$ 以下になる。
    - この性質は、鳩の巣原理(ディリクレの部屋割論法)で証明できるとのこと。
* 解説で、偉そうに「性質が成り立つ」って書いたけれど、証明やどういう原理なのか、実は自分でも分かってなかったりする...
    - これがオイラーの$\phi$関数なのかな？
    - まあ、内容を正確に理解していなくても活用できる点が、定理の利点でもあるのだけれど。
    - それでもやっぱり、きちんと理解していないのは、何だか気持ち悪い...

- - -

[^pe]: About - Project Euler ([https://projecteuler.net](https://projecteuler.net))
[^def]: シェル芸の定義 >> シェル芸 - 上田ブログ ([https://blog.ueda.asia/?page_id=1434](https://blog.ueda.asia/?page_id=1434))
[^ans]: 983