シェルスクリプトで乱数生成 (Xorshift32) を実装する

## はじめに

シェルスクリプトで乱数を扱う場合、`$RANDOM` 変数を使うか、`/dev/random` （または `/dev/urandom`） を使うか、awk の `rand` 関数を使うのが一般的だと思います。しかし `$RANDOM` 変数は POSIX で規定されたシェルの機能ではないので全てのシェルで使えるわけではありません。 `/dev/random` (`/dev/urandom`) は全ての環境で実装されているとは限らずシード値による再現可能な乱数が欲しい場合には使えません。awk を含む外部コマンドは呼び出しが遅いです。といった問題があります。そこで POSIX 準拠のシェルでも簡単に実装可能な Xorshift32 を実装してみました。Xorshiftについては「[Google Chromeが採用した、擬似乱数生成アルゴリズム「xorshift」の数理](https://blog.visvirial.com/articles/575) 」などを参照して下さい。

## C 言語版

計算する値は乱数なので出力結果を見ても正しいかどうかぱっと分からないので、まず C 言語版を実装しました。実装が簡単すぎるせいか手頃なライブラリが見つからなかったので Wikipedia の [Xorshift](https://en.wikipedia.org/wiki/Xorshift#Example_implementation) 等を参考にして作成しています。

```c
#include <stdio.h>
#include <stdint.h>
#include <inttypes.h>

struct xorshift32_state {
  uint32_t a;
};

uint32_t xorshift32(struct xorshift32_state *state)
{
  uint32_t x = state->a;
  x ^= x << 13;
  x ^= x >> 17;
  x ^= x << 5;
  return state->a = x;
}

int main() {
  struct xorshift32_state state = { 2463534242 };
  for (int i = 0; i < 10; i++) {
    printf("%"PRIu32"\n", xorshift32(&state));
  }
  return 0;
}
```

出力

```
723471715
2497366906
2064144800
2008045182
3532304609
374114282
1350636274
691148861
746858951
2653896249
```

注意点として `2463534242` という値はシード値ですので、実際には 32 bit の適当な値を初期値として使います。ただし 0 以外です。計算してみるとわかると思いますが 0 だと何度やっても 0 です。よって取得できる乱数の値も 1 から 4294967295 となります。シード値には時刻や PID を使用すると良いでしょう。詳細は「[シード値について](#シード値について)」を参照してください。

出力数を 10 個から 10000 個 に変更して `sha1sum` した結果（とついでに計算速度）です。これを使って比較するとしましょう。

```console
$ time ./xorshift32 | sha1sum
947d263b6ed277bb90991cad5bdf4e0f0681ea41  -

real    0m0.012s
user    0m0.000s
sys     0m0.031s
```

## シェルスクリプト版

さて C 言語版で正しい結果が分かったので安心してシェルスクリプト版を実装できます。

### 1. 基本

```r
xorshift32() {
  # 計算途中の値はsetで位置パラメータ$1に再代入して無駄にグローバル変数を使わないようにしています
  # 速度の点でいえば変数を使った方が速いです
  set -- $(( $1 ^ (($1 << 13) & 0xFFFFFFFF) ))
  set -- $(( $1 ^ ($1 >> 17) ))
  set -- $(( $1 ^ (($1 << 5) & 0xFFFFFFFF) ))
  xorshift32=$1 # 計算結果をechoではなく変数で返しているのは遅いコマンド置換を避けるためです
}

state=2463534242
for i in $(seq 10000); do
  xorshift32 "$state"
  state=$xorshift32
  echo "$state"
done
```

```console
$ time sh ./xorshift32.sh | sha1sum
947d263b6ed277bb90991cad5bdf4e0f0681ea41  -

real    0m0.107s
user    0m0.078s
sys     0m0.094s
```

C 言語版と同じ出力になりました。使用するのは XOR と ビットシフトだけなので本質的な処理は同じです。違う所は数値を 32 bitの範囲に収めるために `0xFFFFFFFF` でビットマスクをしていることぐらいです。

### 2. 符号付き32bit対策

1.のコードでほとんどのシェルで動作するのですが `mksh` は例外です。`mksh` は OS が 64 bit であっても、（POSIX に厳密に従って）数値を符号付き 32 bit として扱うので 2の31乗-1 を超える値はマイナス値になってしまいます。またテストしてないのですが（したはずなのですが結果を忘れました ^^;） 32bit OSだと他のシェルも同様の結果になると思います。

そこでそれらに対するワークアラウンドを追加します。まず（マイナス値を）右シフトした場合のビットマスク（`& 0x7FFF`）を追加しています。そして計算の最後でマイナスになっている場合は桁を２つに分けて計算し、文字列として結合しています。（ただし一般的にはランダム値をそのまま使うことは少なく何かしらの計算を行ってから使うはずで、結局その時にマイナス値として扱われるので、マイナスをプラスにする処理はいらない気もしなくはないです。）

```r
xorshift32() {
  set -- $(( $1 ^ (($1 << 13) & 0xFFFFFFFF) ))
  set -- $(( $1 ^ (($1 >> 17) & 0x7FFF) ))
  set -- $(( $1 ^ (($1 << 5) & 0xFFFFFFFF) ))

  if [ "$1" -lt 0 ]; then
    set -- $((429496729 + $1 / 10)) $((6 + $1 % 10))
    [ "$2" -lt 0 ] && set -- $(($1 - 1)) $(($2 + 10))
    set -- "${1#0}$2"
  fi
  xorshift32=$1
}
```

### 3. pdksh対策

通常は符号付き32bit対策までやれば十分でしょう。これは古すぎて今では使われてるとは思えない pdksh 対策です。pdksh は `0xFFFF` という `0x` で始まる16進数表記を認識できません。そのため 10 進数に置き換えます。（もしくは pdksh だけ `16#FFFF` という表記を使います。）

```r
xorshift32() {
  set -- $(( $1 ^ (($1 << 13) & 4294967295) ))
  set -- $(( $1 ^ (($1 >> 17) & 32767) ))
  set -- $(( $1 ^ (($1 << 5) & 4294967295) ))

  if [ "$1" -lt 0 ]; then
    set -- $((429496729 + $1 / 10)) $((6 + $1 % 10))
    [ "$2" -lt 0 ] && set -- $(($1 - 1)) $(($2 + 10))
    set -- "${1#0}$2"
  fi
  xorshift32=$1
}
```

以上で完成です。

## シード値について


シード値とは（疑似）乱数を計算するときの初期値です。同じ値を指定すると呼び出すたびに同じ乱数値が生成されます。意図的に同じ乱数を生成する場合を除き異なる値を指定する必要があります。例えば [awk](https://pubs.opengroup.org/onlinepubs/9699919799/utilities/awk.html) でシード値を指定する `srand` 関数は引数を省略した場合に現在の時刻を使用します。


> srand([expr])
>   Set the seed value for rand to expr or use the time of day if expr is omitted. The previous seed value shall be returned.

シェルスクリプトであれば `date +%s` （UNIXタイム）を使うのが手っ取り早いでしょう。ただし `%s` は POSIX では規定されていないため（殆どの環境で対応していますが）対応してない可能性があります。POSIX に完全に準拠したい場合は「[シェルスクリプトでUNIXTIMEの取得と日時との相互変換（POSIX準拠dateコマンドのみ使用）](https://qiita.com/ko1nksm/items/e61728f863c1b2bc829d)」を参照してください。もっともシード値として使用する場合は UNIXタイム である必要はなく、`date +%m%d%H%M%S` （月日時分秒）でも十分だと思います。（注 `%m` は常に 2 桁 で 10 未満の場合に頭に 0 がつくので、計算する際は 8 進数と解釈されないように注意してください。`${value#0}` とすることで頭の 0 は取り除くことが出来ます。）

シード値に秒を使用した場合、プログラムを 1 秒間に何回も起動するとシード値が同じになる可能性があります。これを避けたい場合は PID （`$$` 変数）を組み合わせる（何らかの計算を行う）のが良いと思います。Xorshift32 の場合、シード値は 32 bit なのでこれを超えないように注意してください。

PID はプロセスごとに値が異なるため、シード値に「秒 + PID」を使用した場合に同一マシンで短時間で同じ値になることはまずないでしょう。しかし仮想マシンやコンテナを同時に何百台も起動したら多数の同一のシード値になる可能性は否定できません。（えぇ、考えすぎだと思いますよ？）

そういう場合を懸念する場合に、秒と PID が同じであっても異なる値を返す方法です。仕組みとしては `sleep 0` を実行する間に何回足し算を行えるか？を利用しており、環境の僅かな違いによる誤差で"それなりに"異なる値になることが期待できます。

```sh
#!/bin/sh
set -eu
seed() {
  (
    trap '' INT
    (
      i=0
      trap 'echo $i; exit 0' TERM
      while :; do i=$((i+1)); done
    ) &
    sleep 0
    kill -TERM $!
    wait
  )
}
seed
```

まあ、ここまでシード値にこだわる必要はないと思います。正直 POSIX の範囲でやれることを思いついてしまったから書いただけなので。