ループ・再帰・gotoを使わず『フィボナッチ数列』を印字する

以前、『[「ループ・再帰・gotoを使わず1から100までの数値を印字する」Conner Davis 氏の回答の考察 - Qiita](https://qiita.com/xtetsuji/items/19d07c629852876da401)』 がトレンド入りしていました。
今回はその手法を一般化し、様々な数列を表現してみます。

## フィボナッチ数列を印字する
結論を始めに言います。
 **有理数$\frac{10^d}{(10^{2d}-10^{d}-1)}$の小数点部分にフィボナッチ数列が並びます。**
$d$は各項の表示桁数を示します。

具体的な数値を見てみましょう。
$d=2$のとき、

```math
\begin{align}
\frac{10^2}{(10^{4}-10^{2}-1)}&=\frac{100}{9899} \\
&=0.01010203050813213455\dots
\end{align}
```
$d=3$のとき、

```math
\begin{align}
\frac{10^3}{(10^{6}-10^{3}-1)}&=\frac{1000}{998999} \\
&=0.001001002003005008013021034055\dots
\end{align}
```
となります。
本記事ではこの有理数の導出及び一般化の考察を行います。

せっかくなので、元記事に合わせてPerlでコーディングしてみます。
フィボナッチ数列の100項目は21桁、正しく求めるには計算精度2100桁以上が必要になります。
さすがに難しいので、今回は6桁に収まる30項目まで表示してみます。

```pl:eccentric-sequence-davis.pl
use strict;
use warnings;
use feature qw(say);
use Math::BigFloat;

Math::BigFloat->accuracy(7*30-6);
my $x = 10**7 / Math::BigFloat->new(10**14-10**7-1);
my $str = $x->bstr();

say "--- $str";

$str =~ s/^0\.//;
say join " ", map { $_+0 } $str =~ /(\d{7})/g;
```

```text:output
--- 0.000000100000010000002000000300000050000008000001300000210000034000005500000890000144000023300003770000610000098700015970002584000418100067650010946001771100286570046368007502501213930196418031781105142290832040
1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946 17711 28657 46368 75025 121393 196418 317811 514229 832040
```
`832040`は[正しいです。](https://oeis.org/A000045)

また、元記事を辿ると、本来はc++で実行するのが目的だったようです。
c++で多倍長小数を扱うには、[boostのMultiprecisionライブラリ](https://boostjp.github.io/tips/multiprec-float.html)を使います。(実はAtCoderでも使える。)

```cpp:eccentric-sequence-davis.cpp
#include <stdio.h>
#include <iostream>
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <string>

namespace mp = boost::multiprecision;
using cpp_dec_float = mp::number<mp::cpp_dec_float<210-6>>;

int main()
{
    cpp_dec_float x = 1e7 / cpp_dec_float(1e14-1e7-1);
    std::string str = x.str(210, std::ios_base::fixed);
    std::cout << "--- " << str << std::endl;
  
    for (int i = 0; i < 30; i++) {
        std::cout << std::stoi(str.substr(2+i*7, 7)) << ' ';
    }
}
```

```text:output
--- 0.000000100000010000002000000300000050000008000001300000210000034000005500000890000144000023300003770000610000098700015970002584000418100067650010946001771100286570046368007502501213930196418031781105142290832040
1 1 2 3 5 8 13 21 34 55 89 144 233 377 610 987 1597 2584 4181 6765 10946 17711 28657 46368 75025 121393 196418 317811 514229 832040 
```
結果は同じです。

## 数学的な話。
数列$A(n)$の母関数に$10^d-1$を渡せば目的の有理数となります。
元記事ではその発想の流れについて考察されていましたが、母関数を知っている人であれば、割と自然に出てくる発想だと思います。


## 定数を返す数列を作る
**はじめに、$A(n)=1$となるような数列を作ります。**

たとえば$A(n)=1$を表す有理数は、$0.111\dots$や$0.010101\dots$などです。
これらは、$1/9$や$1/99$と分数表記ができます。

**$d$桁ずつの表記を得たい場合、$\frac{1}{10^d-1}$が目的の有理数になりそうです。**
（ここで出てくる10は、十進表記という意味を持ちます。）

$d=4$のとき、$$\frac{1}{10^4-1}=0.000100010001\dots$$となるので、正しそうです。

厳密な証明とは言えませんが、導出を簡単に示しておきます。

目的の数列を$A(n)=1$、対応する有理数を$A_d$とします。このとき$d$桁ずつの表記になるとします。
**新しい数列$B(n)=A(n)-A(n+1)$とする**と、以下のようになります。
(以後簡単のため、$n<0$で$A(n)=0$としています。)

```math
\begin{align}
B(n) &=A(n)-A(n+1) \\
&= [1, 1, 1, 1, \dots] - [0, 1, 1, 1, \dots] \\
&= [1, 0, 0, 0, \dots] \\
&= \left\{
\begin{array}{ll}
1 & (n = 0)\\
0 & (n \neq 0)
\end{array}
\right. \\
&=\delta_{0,n}
\end{align}
```
**ここで$B(n)$に対応する有理数は、$B_d=\frac{1}{10^d}$です。**

$B(n)=A(n)-A(n+1)$という関係を、有理数の領域で表現すると、
$$B_d=A_d-10^{-d}A_d$$
となります。**$10^{-d}$を掛けることで、有理数を小数点下方向に$d$桁分シフトできます。**

```math
\begin{align}
A_d-10^{-d}A_d &= B_d \\
(1-10^{-d})A_d &= B_d \\
A_d &= \frac{B_d}{1-10^{-d}}
\end{align}
```
$B_d=A_d-10^{-d}A_d$なので、

```math
\begin{align}
A_d &= \frac{1}{10^d(1-10^{-d})}\\
 &= \frac{1}{10^d-1}
\end{align}
```

よって目的の有理数$A_d=\frac{1}{10^d-1}$が得られました。

またこれを用いて、**数列$A(n)=a$(定数)に対応する有理数は、$A_d=\frac{a}{10^d-1}$となります。**

## 自然数を返す数列を作る(元ネタ)
元ネタにあるような、**$A(n)=n$となるような数列を作ります。**
**新しい数列$B(n)=A(n)-A(n+1)$とする**と、以下のようになります。

```math
\begin{align}
B(n) &=A(n)-A(n+1) \\
&= [1, 2, 3, 4, \dots] - [0, 1, 2, 3, \dots] \\
&= [1, 1, 1, 1, \dots]
\end{align}
```

このとき$B(n)$は、上で求めた通り$\frac{1}{10^d-1}$と表現できます。
よって目的の、自然数列を表す有理数$A_d$は次のように求まります。

```math
\begin{align}
A_d-10^{-d}A_d &= B_d \\
(1-10^{-d})A_d &= B_d \\
A_d &= \frac{B_d}{1-10^{-d}} \\
&= \frac{1}{(10^d-1)(1-10^{-d})} \\
&= \frac{10^d}{(10^d-1)^2} \\
\end{align}
```
これは[元記事](https://www.quora.com/How-can-I-print-1-to-100-in-C%2B%2B-without-a-loop-goto-or-recursion?ch=10&share=9b09c507&srid=Jqvc)のConner Davis 氏の回答に一致しています。

またこれを用いて、**aの倍数の数列$A(n)=an$に対応する有理数は、$A_d=\frac{10^da}{(10^d-1)^2}$となります。**

## フィボナッチ数列を作る
この記事のタイトル回収です。
**$A(n)=A(n-1) + A(n-2), A(0)=A(1)=1$となるような数列を作ります。**
**新しい数列$B(n)=A(n) - A(n-1)$とする**と、以下のようになります。

```math
\begin{align}
B(n) &= A(n) - A(n-1) \\
&= [1, 1, 2, 3, 5, 8, \dots] - [0, 1, 1, 2, 3, 5,\dots] \\
&= [1, 0, 1, 1, 2, 3, 5, \dots] \\
&= \left\{
\begin{array}{ll}
1 & (n = 0)\\
A(n-3) & (n \neq 0)
\end{array}
\right. \\
&= \delta_{0,n} + A(n-3)
\end{align}
```

このとき$\delta_{0,n}$は、上で求めた通り$\frac{a}{10^d}$と表現できます。
よって目的の、フィボナッチ数列を表す有理数$A_d$は次のように求まります。

```math
\begin{align}
A_d-10^{-d}A_d &= B_d \\
&= \frac{1}{10^d} + 10^{-2d}A_d \\
(1-10^{-d}-10^{-2d})A_d &= \frac{1}{10^d} \\
(10^{2d}-10^{d}-1)A_d &= 10^d \\
A_d &= \frac{10^d}{(10^{2d}-10^{d}-1)} \\
\end{align}
```

具体的な例は、記事の冒頭に挙げています。

## その他の数列を作る
証明が~~めんｄ~~面白かったので、以下の結果の証明は*読者への課題*とします。
(気が向いたら追記します。)

### 累乗
$A(n)=a^n$のとき、
$$A_{a, d}=\frac{1}{10^d-a}$$
例：$n=3, d=2$

```math
\begin{align}
A_{3, 2}&=\frac{1}{10^{2}-3} \\
&=\frac{1}{100-3} \\
&=\frac{1}{97} \\
&=0.01030927\dots
\end{align}
```

### 平方数
$A(n)=n^2$のとき、
$$A_d=\frac{10^d+1}{10^{3d}-3*10^{2d}+3*10^d-1}$$
例：$d=3$

```math
\begin{align}
A_3&=\frac{10^3+1}{10^9-3*10^6+3*10^3-1} \\
&=\frac{1001}{1000000000-3000000+3000-1} \\
&=\frac{1001}{997002999} \\
&=0.000001004009016025036049064081100121\dots
\end{align}
```

### 立方数・n乗数
$A(n)=n^3$のとき、
$$A_d=\frac{10^{2d}+4*10^d+1}{10^{4d}-4*10^{3d}+6*10^{2d}-4*10^{d}+1}$$
例：$d=3$

```math
\begin{align}
A_4&=\frac{10^{8}+4*10^4+1}{10^{16}-4*10^{12}+6*10^{8}-4*10^{4}+1} \\
&=\frac{100040001}{9996000599960001} \\
&=0,000000010008002700640125021603430512072910001331\dots
\end{align}
```
以降$A(n)=n^m$は、二項係数を使って一般化できます。



### カタラン数
$A(n)$がカタラン数のとき、
$$A_d=\frac{1-\sqrt{4-10^{-d}}}{2*10^{-d}}$$
例：$d=3$

```math
\begin{align}
A_3&=\frac{1-\sqrt{4-10^{-3}}}{2*10^{-3}} \\
&=\frac{1-\sqrt{4-0.0001}}{0.0002} \\
&=1.000100020005001400420132042914304\dots
\end{align}
```