【Excel not VBA】Excel で xoshiro128**を書いてみた

どうも、3 月 5 日を勝手に「世界ビット演算デー」にしたいたぬきです。なんで 3 月 5 日なのかって? $3_{10} = 11_2, 5_{10} = 101_2$ なので並べるだけで真理値表になるからです。ちなみに実際は「軍縮不拡散に関する啓発のための国際デー」らしいです。めちゃシリアスじゃん。

本記事は [misskey.dev ユーザー Advent Calendar 2023](https://adventar.org/calendars/9688) の 13 日目の記事でもあります。Qiitadon 爆破事件のあと dev に流れ着いた身なので、まあ Qiita に投げ込んでも問題ないやろ。

misskey.dev にまつわる記事と思いきや、

> misskey.dev のユーザーならなんでもどうぞ

とのことなので、好き勝手書きます。私が好きなものといえばビット演算と Excel です。

## 初期空間を指定できる乱数関数がほしい

Excel には乱数を得る関数が色々ありますが、シード値を設定することができません。モンテカルロ法のシミュレーションで Excel を使うとすると、出た乱数を全部「値としてコピー」しなきゃいけないわけです。複数パターンの検証とか大変ですよね。

じゃあ自分で実装すればいいじゃない。というわけで xorshift 派生アルゴリズムの xoshiro ファミリーを実装していきます。Excel の数値型は(通貨を除いて)おおむね binary64 なので、32 bits 整数範囲は正確に表現できるはずです。よって xoshiro128**を実装します。

VBA ではなく LAMBDA で実装するため、マクロ禁止環境でも利用できます。

## ユーティリティ関数ども

まず先に実装の際に便利ないくつかの関数を定義します。

### AS32BITS

32 bits 整数範囲になるようマスクする関数です。binary64 では 53 bits までなら正確に表現できますが、今回の要件の場合シフト命令とかで正確に表現されると困るので、この関数を使います。

```vb:AS32BITS
= LAMBDA(var_n, BITAND(var_n, 4294967295))
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/241027/e22a24c3-bd08-4d16-8875-896ec8a97203.png)

### BITLSHIFT32

`BITLSHIFT` 後 `AS32BITS` してくれる関数です。それだけなのですが、毎回関数を二つ書くのは面倒なので(あと書き忘れとか怖い)この関数を使います。

```vb:BITLSHIFT32
= LAMBDA(var_n, var_width, AS32BITS(BITLSHIFT(var_n, var_width)))
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/241027/1eea3f3a-9e1c-5572-012e-b3c0a96a9c1d.png)

### BITLROTATE32

`BITLSHIFT32` してはみ出たビットを右につけたす、いわゆるローテート演算を行う関数です。ローテートを利用するのが xoshiro ファミリーの特徴です。

```vb:BITLROTATE32
= LAMBDA(var_n, var_width, BITOR(BITLSHIFT32(var_n, var_width), BITRSHIFT(var_n, 32 - var_width)))
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/241027/fc0c1c95-bcfd-5143-d6cb-b405b430a7ab.png)

### FLATINDEX

Excel に存在する `INDEX` 関数は行と列の添字をそれぞれ指定しますが、これを Z 字にたどることを想定して一つの添字で参照するための関数です。

<table>
  <tr>
    <td>(1, 1)</td><td>(1, 2)</td><td>(1, 3)</td>
  </tr>
  <tr>
    <td>(2, 1)</td><td>(2, 2)</td><td>(2, 3)</td>
  </tr>
  <tr>
    <td>(3, 1)</td><td>(3, 2)</td><td>(3, 3)</td>
  </tr>
</table>

これが公式の `INDEX` ですが、これを

<table>
  <tr>
    <td>1</td><td>2</td><td>3</td>
  </tr>
  <tr>
    <td>4</td><td>5</td><td>6</td>
  </tr>
  <tr>
    <td>7</td><td>8</td><td>9</td>
  </tr>
</table>

のように参照できるようにするのが `FLATINDEX` です。

`SEQUENCE` 関数と併用することでこのあとに定義する `RESHAPE` 関数の実装を容易にします。

```vb:FLATINDEX
= LAMBDA(var_array, var_i, LET(var_cols_num, COLUMNS(var_array), var_row_i, QUOTIENT(var_i - 1, var_cols_num) + 1, var_col_i, MOD(var_i - 1, var_cols_num) + 1, INDEX(var_array, var_row_i, var_col_i)))
```

これだと処理を追いづらいので整形すると

```vb:FLATINDEX
= LAMBDA(var_array, var_i,
    LET(
        var_cols_num, COLUMNS(var_array),
        var_row_i,    QUOTIENT(var_i - 1, var_cols_num) + 1,
        var_col_i,    MOD(var_i - 1, var_cols_num) + 1,

        INDEX(var_array, var_row_i, var_col_i)
    )
)
```

といった感じです。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/241027/692432e1-4a0b-7362-fe1b-a6e3a9700299.png)

### RESHAPE

配列の形状を要素数が同じ別の形状に変更します。とりあえず 1 列でバーッと乱数を出したあと所望の形状に変更するために使います。

```vb:RESHAPE
= LAMBDA(var_array, var_rows_num, var_cols_num, FLATINDEX(var_array, SEQUENCE(var_rows_num, var_cols_num)))
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/241027/594d7663-c6f3-37e1-d488-7e383b94e5c9.png)

## 状態遷移関数

いよいよ xoshiro128**を実装していきます。といってもリファレンス実装を Excel 数式で書き下していくだけです。ただ、Excel 数式はプログラミング風にいうと宣言的なので、変数の書き換えとかはできません。よって一時変数をふたつ増やして実装しています。

必要な行ぶんの状態を一気に出していきます。繰り返しは使えないので再帰で実装します。1 行だけの場合は与えられた空間をそのまま吐きます。2 行以上の場合は再帰で 1 行前までの状態を吐き、その下に新しい行を追加します。

ところでこの「行を追加する」、`HSTACK` / `VSTACK` 関数ができてほんとに楽に書けるようになりました。配列リテラルは中身もリテラルじゃないといけない(数式を書けない)問題などがあり、かつては [`SEQUENCE` 関数の結果を `IFS` 関数で書き換えるみたいなハック](https://qiita.com/yosatonet/items/47a21439702e6485f054)が必要でした。

```vb:XOSHIRO128STATE
= LAMBDA(var_states_num, var_state_1, var_state_2, var_state_3, var_state_4, IF(var_states_num = 1, HSTACK(AS32BITS(var_state_1), AS32BITS(var_state_2), AS32BITS(var_state_3), AS32BITS(var_state_4)), LET(var_old, XOSHIRO128STATE(var_states_num - 1, var_state_1, var_state_2, var_state_3, var_state_4), var_old_state_1, INDEX(var_old, var_states_num - 1, 1), var_old_state_2, INDEX(var_old, var_states_num - 1, 2), var_old_state_3, INDEX(var_old, var_states_num - 1, 3), var_old_state_4, INDEX(var_old, var_states_num - 1, 4), var_t, BITLSHIFT32(var_old_state_2, 9), var_u, BITXOR(var_old_state_3, var_old_state_1), var_v, BITXOR(var_old_state_4, var_old_state_2), var_new_state_1, BITXOR(var_old_state_1, var_v), var_new_state_2, BITXOR(var_old_state_2, var_u), var_new_state_3, BITXOR(var_u, var_t), var_new_state_4, BITLROTATE32(var_v, 11), var_new, HSTACK(var_new_state_1, var_new_state_2, var_new_state_3, var_new_state_4), VSTACK(var_old, var_new))))
```

これも整形すると、

```vb:XOSHIRO128STATE
= LAMBDA(var_states_num, var_state_1, var_state_2, var_state_3, var_state_4,
    IF(var_states_num = 1,
        HSTACK(AS32BITS(var_state_1), AS32BITS(var_state_2), AS32BITS(var_state_3), AS32BITS(var_state_4)),
        LET(
            var_old,         XOSHIRO128STATE(var_states_num - 1, var_state_1, var_state_2, var_state_3, var_state_4),
            var_old_state_1, INDEX(var_old, var_states_num - 1, 1),
            var_old_state_2, INDEX(var_old, var_states_num - 1, 2),
            var_old_state_3, INDEX(var_old, var_states_num - 1, 3),
            var_old_state_4, INDEX(var_old, var_states_num - 1, 4),
            var_t,           BITLSHIFT32(var_old_state_2, 9),
            var_u,           BITXOR(var_old_state_3, var_old_state_1),
            var_v,           BITXOR(var_old_state_4, var_old_state_2),
            var_new_state_1, BITXOR(var_old_state_1, var_v),
            var_new_state_2, BITXOR(var_old_state_2, var_u),
            var_new_state_3, BITXOR(var_u, var_t),
            var_new_state_4, BITLROTATE32(var_v, 11),
            var_new,         HSTACK(var_new_state_1, var_new_state_2, var_new_state_3, var_new_state_4),

            VSTACK(var_old, var_new)
        )
    )
)
```

となります。いやだいぶ長いな…

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/241027/b3b5899f-fc23-c63b-9a93-ae7920e4bc61.png)

## 出力関数

次にこれらの内部空間から乱数の出力を得る関数を実装します。xoshiro128**では内部空間のうち 2 番目(1 オリジンで)のみを使います。

```vb:XOSHIRO128STARSTAROUTPUT
= LAMBDA(var_states, LET(var_states_num, ROWS(var_states), var_state_2_s, MAP(SEQUENCE(var_states_num), LAMBDA(row_i, INDEX(var_states, row_i, 2))), AS32BITS(BITLROTATE32(AS32BITS(var_state_2_s * 5), 7) * 9)))
```

整形すると、

```vb:XOSHIRO128STARSTAROUTPUT
= LAMBDA(var_states,
    LET(
        var_states_num, ROWS(var_states),
        var_state_2_s,  MAP(SEQUENCE(var_states_num), LAMBDA(row_i, INDEX(var_states, row_i, 2))),

        AS32BITS(BITLROTATE32(AS32BITS(var_state_2_s * 5), 7) * 9)
    )
)
```

となります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/241027/e4209d82-36a9-99ce-b29e-2dcd9992a2b3.png)

なお、`OFFSET` 関数を使わず `SEQUENCE` 関数の戻り値に `MAP` 関数をかけるというややまだるっこしいことをしていますが、`OFFSET` 関数は、`LET` 関数で束縛された式などのような「実際にシート上にあるわけではない配列」に対してエラーを返す場合があるなど、やや動作が非直観的なため、ここでは使用していません。

## UI

最後に、ユーザーが実際に使用するための関数を実装します。

```vb:XOSHIRO128STARSTAR
= LAMBDA(var_rows_num, var_cols_num, var_seed, LET(var_rands_num, var_rows_num * var_cols_num, var_seed_is_array, ROWS(var_seed) * COLUMNS(var_seed) > 1, var_seed_1, IF(var_seed_is_array, FLATINDEX(var_seed, 1), var_seed), var_seed_2, IF(var_seed_is_array, FLATINDEX(var_seed, 2), MOD(48271 * var_seed_1, 2147483647)), var_seed_3, IF(var_seed_is_array, FLATINDEX(var_seed, 3), MOD(48271 * var_seed_2, 2147483647)), var_seed_4, IF(var_seed_is_array, FLATINDEX(var_seed, 4), MOD(48271 * var_seed_3, 2147483647)), var_states, XOSHIRO128STATE(var_rands_num, var_seed_1, var_seed_2, var_seed_3, var_seed_4), var_rands, XOSHIRO128STARSTAROUTPUT(var_states), RESHAPE(var_rands, var_rows_num, var_cols_num)))
```

整形します。

```vb:XOSHIRO128STARSTAR
= LAMBDA(var_rows_num, var_cols_num, var_seed,
    LET(
        var_rands_num,     var_rows_num * var_cols_num,
        var_seed_is_array, ROWS(var_seed) * COLUMNS(var_seed) > 1,
        var_seed_1,        IF(var_seed_is_array, FLATINDEX(var_seed, 1), var_seed),
        var_seed_2,        IF(var_seed_is_array, FLATINDEX(var_seed, 2), MOD(48271 * var_seed_1, 2147483647)),
        var_seed_3,        IF(var_seed_is_array, FLATINDEX(var_seed, 3), MOD(48271 * var_seed_2, 2147483647)),
        var_seed_4,        IF(var_seed_is_array, FLATINDEX(var_seed, 4), MOD(48271 * var_seed_3, 2147483647)),
        var_states,        XOSHIRO128STATE(var_rands_num, var_seed_1, var_seed_2, var_seed_3, var_seed_4),
        var_rands,         XOSHIRO128STARSTAROUTPUT(var_states),

        RESHAPE(var_rands, var_rows_num, var_cols_num)
    )
)
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/241027/fefc1bf3-a298-5bd1-cf8d-96167548e587.png)

要求量の乱数を生成して `RESHAPE` します。シードはひとつの値でも 4 個(以上)の値を持つ配列でも構いません。ひとつの値を与えた場合、残りの空間は線形合同法で埋めることにしました(正直微妙ではある)。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/241027/56e1724f-d8d4-5e89-81de-bfbc5ae3fb42.png)

リファレンス実装の結果と比較。問題なさそうです。

## おまけ

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/241027/2fc43b58-9d2d-ae46-086f-2c6e2415fd11.png)

ほぎゃーーーーーっ!!

はい、Excel の再帰は割と貧弱なのでそんなにたくさんの乱数は作れません[^recursive_limit]。だめじゃん。

[^recursive_limit]: 今回の関数の場合 861 個が限界でした。