IchigoJamで階乗を求める

# 何をするか

IchigoJam BASIC を用い、階乗 (1からその数までの正の整数を1個ずつ掛けた値) を求める。
なるべく大きな数に対応する。

# 設計

## 数の表現

変数1個では表現できない大きな数を表現するには、数を数桁ずつに区切って格納するのが基本である。

16進数の桁で区切ると、コンピュータ上の数の表現と相性が良く、格納の効率は良くなる。
しかし、10進数で出力したい場合は10進数に変換しなければならず、処理が面倒になりそうである。

そこで、今回は10進数の桁で区切ることにした。
配列領域だけでなくキャラクターパターンの領域も活用して数を格納するため、
変数(2バイト)ではなく1バイト単位で処理を行うことにし、1バイトに10進数の2桁を格納することにした。

## 掛け算の処理

階乗の値は、小さい値の階乗でもどんどん大きくなっていくことが知られている。
従って、入力としては小さい値(大きくても変数1個分)のみをサポートしておけばよさそうである。
よって、今回の掛け算は「多倍長整数×普通の整数」ができればよい。

この処理を行うためには、掛けられる数の単位(今回はバイト)にそれぞれ掛ける数を掛け、
繰り上がりの処理を行うとよい。
すなわち、前の単位からの繰り上がりを$c_0$、今の単位(掛ける前)を$d_0$、掛ける数を$m$、
今の単位(掛けた結果)を$d$、次の単位への繰り上がりを$c$とすると、今回は1単位に10進数2桁を格納するので、
$p = d_0 \times m + c_0$として、$d = p \mod 100$、$c = \lfloor p / 100 \rfloor$となる。

ただし、$d_0$は最大99になる可能性があるので、$m \ge 331$になると、
$d_0 \times m$が32,767を超えてオーバーフローしてしまう可能性が出てくる。
(実際は繰り上がりもあるので、さらに小さい$m$でオーバーフローする可能性がある)

そこで、$m = m_1 \times 100 + m_0 (0 \le m_0 \le 99)$ と分割する。
すると、$p = d_0 (m_1 \times 100 + m_0) + c_0 = d_0 \times m_1 \times 100 + d_0 \times m_0 + c_0$ となる。
これを$d$と$c$を求める式に代入すると、
$(d_0 \times m_1 \times 100) \mod 100 = 0$、$(d_0 \times m_1 \times 100) / 100 = d_0 \times m_1$ なので、
$d =(d_0 \times m_0 + c_0) \mod 100$、$c = d_0 \times m_1 + \lfloor (d_0 \times m_0 + c_0) / 100 \rfloor$ となる。

$d_0 \times m_0 + c_0$ がオーバーフローしないためには、$d_0$も$m_0$も99以下なので、
$99 \times 99 + c_0 \le 32767$ すなわち $c_0 \le 22966$ であればよい。
一番下の単位については、その前の単位からの繰り上がりは無いため、$c_0=0$となり、条件を満たす。
それ以外の単位については、$c \le 99 \times m_1 + \lfloor 32767 / 100 \rfloor \le 22966$
すなわち $m_1 \le 228$ であれば条件を満たす。
言い換えれば、この計算方法であれば $m \le 22899$ においてオーバーフローせずに計算ができることになる。

## 対応範囲の設定

計算結果を格納する領域さえ十分あれば22,899までを掛ける処理がオーバーフローせずにできることがわかったが、
実際は計算結果を格納する領域に制限がありそうである。
今回は、キャラクターパターン領域と変数・配列領域のうち、
計算に使う変数の領域を除いた部分を計算結果の格納に用いる。
実装してみたところ変数を5個使ったので、
512バイトから2バイトの変数5個分の10バイトを除いた502バイトが計算結果の格納に使えることになる。
すなわち、1バイトに2桁を格納するので、1004桁までの計算結果を扱うことができる。

そこで、以下の関数を用いて、Pythonのインタラクティブモードで階乗の桁数を求めてみた。
`x`が階乗を求める関数、`xx`がその桁数を求める関数である。

```py
def x(n):
  r = 1
  for i in range(1, n + 1):
    r *= i
  return r

def xx(n):
  return len(str(x(n)))
```

その結果、`451!`の桁数は1003桁、`452!`の桁数は1006桁と求まった。
よって、1004桁に収まるのは`451!`までである。

ちなみに、`22899!`は89,894桁らしい。

# 実装

```
10 ' ｶｲｼﾞｮｳ ｦ ﾓﾄﾒﾙ
20 INPUT"0 ｲｼﾞｮｳ 451 ｲｶ ﾉ ｶｽﾞ : ",Z
30 IFZ<0OR451<Z?"ｿﾉ ｶｽﾞ ﾊ ｼﾞｮｳｹﾝ ﾆ ｱﾜﾅｲ":GOTO20
40 POKE#700,1:FORY=#701TO#8F5:POKEY,0:NEXT
50 IFZ=0GOTO90
60 FORY=1TOZ
70 X=0:FORW=#700TO#8F5:V=Y%100*PEEK(W)+X:X=Y/100*PEEK(W)+V/100:POKEW,V%100:NEXT
80 NEXT
90 ?Z;"! = ";
100 X=0:FORY=#8F5TO#700STEP-1
110 W=PEEK(Y)
120 IFX=0GOTO150
130 IFW<10?"0";
140 ?W;:GOTO160
150 IFW>0?W;:X=1
160 NEXT
170 ?CHR$(10);
```

# 実行結果

以下のように、入力した数の階乗を求めることができた。

![319の階乗を求めた結果](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/99656/5463894a-356b-1765-a700-de5b64a678bc.png)

今回の出力の仕様で32×24の1画面に収まるのは、332の階乗までのようである。

![332の階乗を求めた結果](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/99656/99e57f45-e2b2-0755-f1c8-a3c66c0dafe0.png)


プログラムに

```
35 CLT
175 ?TICK()
```

を追加して1回測定した実行時間は、以下のようになった。
(`-`は時間がかかりすぎることが予想されるため、測定を行わなかったことを表す)

IchigoKamuy のファームウェアは 1.4.1 (`VER()=14114`)、
IchigoJam R のファームウェアは 1.5b (`VER()=15001`) である。

|階乗を求める値|IchigoKamuy<br>ビデオ出力有|IchigoKamuy<br>ビデオ出力無|IchigoJam R<br>ビデオ出力有|IchigoJam R<br>ビデオ出力無|
|--:|--:|--:|--:|--:|
|10|2059|1351|250|219|
|20|3868|2539|470|412|
|50|9296|6101|1130|992|
|100|18344|12040|2230|1957|
|150|27394|17978|3329|2922|
|200|-|23919|4428|3887|
|250|-|29861|5528|4853|
|300|-|-|6628|5818|
|350|-|-|7728|6783|
|400|-|-|8828|7749|
|451|-|-|9950|8733|

ビデオ出力有のとき、IchigoJam R では`451!`が3分弱で求まったのに対し、
IchigoKamuy ではその約3分の1の計算を行う`150!`を求めるのに約7分半かかっている。
IchigoJam R の圧倒的な速さがうかがえる。

なお、パソコン上のPythonでは、`451!`が一瞬(1秒未満)で求まった。

# おわりに

* IchigoJamはjig.jpの登録商標です。
* IchigoKamuyは株式会社syushuの登録商標です。