PARI/GPと「天に向かって続く数」

[数学 Advent Calendar 2016](http://qiita.com/advent-calendar/2016/math) 2日目の記事です。

# 天に向かって続く数とは

$$
\dots 11000557423423230896109004106619977392256259918212890625
$$

のように無限に大きな方向に伸びる列のことで、10進展開を以下のような集合として定義することができます。

$$
S_{10} = \bigg\\{ \sum_{n=0}^{\infty} a_n10^n \ \bigg| \ 0 \le a_n \lt 10, a_n \in \mathbb{Z} \bigg\\}
$$

より詳細な定義や数学的バックグラウンドについては、こちらの書籍を参照して下さい。

[天に向かって続く数 / 加藤文元 中井保行 著 - 日本評論社](https://www.nippyo.co.jp/shop/book/7214.html)

書籍の冒頭で触れられていますが、この集合の中で上記のような特徴のある列が存在します。下 $n$ 桁を2乗してみましょう。

$$
\begin{align}
 90625^2 &= 8212890625 \equiv 90625 \pmod{10^5} \\\\
 8212890625^2 &= 67451572418212890625 \equiv 8212890625 \pmod{10^{10}}
\end{align}
$$

このように $x^2 = x$ の解となりそうな $0, 1$ 以外の列が得られます。
今回はこのような列を生成する方法を PARI という代数計算用の言語で見ていきましょう。

PARIとそのインタラクティブシェル `gp` のインストールは[公式サイト](http://pari.math.u-bordeaux.fr/)を参照してください。
macOSならHomebrewで一発です。

```
$ brew install pari
```

## 中国剰余定理による剰余の計算

今回の計算では使うものは中国剰余定理のみです。
PARIでは `chinese(m, n)` として定義されており、 `m,n` は剰余数 `Mod(k, p)` つまり、$k \pmod{p}$ です。
有名な「3で割ると2余り、5で割ると3余り、7で割ると2余る数は何か」という問題は次のように計算できます。

```
? chinese(chinese(Mod(2,3), Mod(3,5)), Mod(2,7))
%1 = Mod(23, 105)
```

2進展開で $m$, 5進展開で$n$ であるような10進展開の列は、中国剰余定理により導くことができます。
(一意性などの証明はここでは飛ばします）

逆に、10進展開の数は2進展開、5進展開の数へそれぞれ射影が可能です。

$$
f : S_{10} \to S_5 : x \mapsto \sum_{n=0}^{\infty} b_n5^n \\\\
b_n = ((x \bmod 5^{n+1}) - (x \bmod 5^n)) / 5^n
$$

10進展開の数は $10^{n+1}$ 以上の項は全て $5^{n+1}$ で割り切れることに注意すると、 $b_n$ はたかだか有限個の項で計算できます。

さてこれを使って、10進展開で $x^2 = x$ の解になるようなものを見つけていきます。

## 計算

2進展開、5進展開では $x = 0,1$ が $x^2 = x$ の自明な解であり、10進展開で $0 = \dots 00000$ や $1 = \dots 00001$ ならば2進展開、5進展開でも0や1であるのは自明です。
それならば2つがねじれた状態（2進で1, 5進で0）のような状態であれば $x^2 = x$ の解でありながら非自明な列になるのでは？という発想が生まれます。実際に計算してみましょう。

```
? chinese(Mod(1,2^10),Mod(0,5^10))
%2 = Mod(8212890625, 10000000000)
```

先程の値が計算されました。逆のねじれはどうでしょうか？

```
? chinese(Mod(0,2^10),Mod(1,5^10))
%3 = Mod(1787109376, 10000000000)
```

実はこの列も $x^2 = x$ の解です。
次のように表記しておきましょう。

$$
(1_{(2)}, 0_{(5)}) = \dots 8212890625 \\\\
(0_{(2)}, 1_{(5)}) = \dots 1787109376
$$

もっと長い列を確認してみたい場合は指数を増やすと良いです。

```
? chinese(Mod(0,2^500),Mod(1,5^500))
%17 = Mod(33676827594484243897647456005000654391916198809258469939943944255181290307214900224081949924583571472291837988649753193941836723828323234739062471943155785513806039500165527193278093329582759905765533802187573309212464055383301491935363862833615950970780658118090418340475522138153859087121701561568296751826571113427262336853480895011970552039185326239496042803106285328198624380944537003185235736096046992680891830197061490109937833490419136188999442576576769103890995893380022607743740081787109376, 100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000)
```

剰余の値だけを見たい場合は `lift` を使います。

```
? lift(chinese(Mod(0,2^500),Mod(1,5^500)))
%18 = 33676827594484243897647456005000654391916198809258469939943944255181290307214900224081949924583571472291837988649753193941836723828323234739062471943155785513806039500165527193278093329582759905765533802187573309212464055383301491935363862833615950970780658118090418340475522138153859087121701561568296751826571113427262336853480895011970552039185326239496042803106285328198624380944537003185235736096046992680891830197061490109937833490419136188999442576576769103890995893380022607743740081787109376
```

かなり巨大な数も計算できてしまいます。

続いて、$x^3 = x$ には $x = -1, 0, 1$ という解が存在します。 $0,1$ の非自明な組み合わせは先程試しましたので、その他の組み合わせを計算してみましょう。

```
? chinese(Mod(0,2^10),Mod(-1,5^10))
%21 = Mod(8212890624, 10000000000)
? chinese(Mod(0,2^10),Mod(-1,5^10))^2
%22 = Mod(1787109376, 10000000000)
? chinese(Mod(0,2^10),Mod(-1,5^10))^3
%23 = Mod(8212890624, 10000000000)
```

$x^2 \ne x$ かつ $x^3 = x$ となる列が生成されました。ところで2乗した時の数、どこかで見ませんでした？

$$
(0_{(2)}, -1_{(5)}) = \dots 8212890624 \\\\
(0_{(2)}, -1_{(5)})^2 = (0_{(2)}, 1_{(5)}) = \dots 1787109376
$$

ちゃんと対応が取れていますね。だとすると次のような計算ができてしまうのでは…？

$$
(1_{(2)}, -1_{(5)})^2 = (1_{(2)}, 1_{(5)}) = \dots 0000000001 = 1_{(10)}
$$

```
? chinese(Mod(1,2^10),Mod(-1,5^10))
%24 = Mod(6425781249, 10000000000)
? chinese(Mod(1,2^10),Mod(-1,5^10))^2
%25 = Mod(1, 10000000000)
```

正しいようです。 $x^2 = 1$ の非自明な解になっています。

## まとめ

このような無限列の法則の発見と計算を手元ですぐに試せるのはなかなか楽しいです。
PARIは他にも代数計算のための数多くの機能を備えていますので、気になった方は触れてみてはいかがでしょうか？