SSE4.2: RDTSC と CRC32 の２命令だけで乱数らしきものを作る

# 値の生成方法

生成処理は、Intel IA32 の SSE4.2 に対応したプロセッサ命令で

```asm:生成プログラム
       RDTSC
       CRC32   EAX, EDX
```

二命令だけです。 EAX の値が乱数のようなものになるだろうという期待をします。

## RDTSC と CRC32 の関係

CPU の CRC32 は元となる crc32 計算

<dl><dd>eax = crc32(edx)</dd></dl>

とは少し違い

<dl><dd>eax = crc32(eax ^ edx)</dd></dl>

という計算をします。（レジスタのビット反転を含むので詳しくは仕様書を参照してください）

元の y = CRC32(x) では、複数の x に対して、y が１つ（重複）となることがありません。符号なし32ビット整数における全単射になるでしょうか。([過去記事](https://qiita.com/ikiuo/items/492a133d05edc55e4543))

RDTSC は EDX:EAX にプロセッサのタイムスタンプ（分解能は CPU クロック) を 64 ビットで返します。EDX は 2<sup>32</sup> クロック置きに更新されます。4GHz 動作で約 1 秒間隔です。

これを生成プログラムのように CRC32 に渡すと、RDTSC[EDX:EAX] と CRC32 の関係は

| EDX | EAX | CRC32 |
|--:|--:|--:|
| <code>00000000</code> | <code>00000000</code> | <code>00000000</code> |
| <code>00000000</code> | <code>00000001</code> | <code>dd45aab8</code> |
| <code>00000000</code> | <code>00000002</code> | <code>bf672381</code> |
| <code>00000000</code> | <code>00000003</code> | <code>62228939</code> |
| ... | ... | ... |
| <code>00000000</code> | <code>fffffffe</code> | <code>6add1e80</code> |
| <code>00000000</code> | <code>ffffffff</code> | <code>b798b438</code> |
| <code>00000001</code> | <code>00000000</code> | <code>dd45aab8</code> |
| <code>00000001</code> | <code>00000001</code> | <code>00000000</code> |
| <code>00000001</code> | <code>00000002</code> | <code>62228939</code> |
| <code>00000001</code> | <code>00000003</code> | <code>bf672381</code> |
| ... | ... | ... |
| <code>fffffffe</code> | <code>fffffffe</code> | <code>00000000</code> |
| <code>fffffffe</code> | <code>ffffffff</code> | <code>dd45aab8</code> |
| <code>ffffffff</code> | <code>00000000</code> | <code>b798b438</code> |
| <code>ffffffff</code> | <code>00000001</code> | <code>6add1e80</code> |
| ... | ... | ... |
| <code>ffffffff</code> | <code>fffffffe</code> | <code>dd45aab8</code> |
| <code>ffffffff</code> | <code>ffffffff</code> | <code>00000000</code> |

<details><summary>CRC32命令の計算</summary><div>

```python:sse42crc32.py
def crc32(eax, edx):
    eax ^= edx
    for _ in range(32):
        f = eax & 1
        eax >>= 1
        if not f:
            continue
        eax ^= 0x82f63b78
    return eax
```

</div></details>

となるので、<code> TSC = (EDX << 32) | EAX </code> に対する CRC32[TSC] を乱数列（質に疑問符あり）と捉えます。周期は 2<sup>64</sup> CPU サイクルで、4GHz 動作だと百数十年になります。

生成プログラムは、呼び出された時点の RDTSC（時間とともに増える）から CRC32[TSC] で引き出すことになります。結果、乱数列臭いものが出てくることは予想できますが、実行状況によっては乱数列とよべないものが並ぶことも考えられます。CPU の動作に詳しい人なら、狙ったタイミングで RDTSC を実行する難しさはご存じでしょう。よって、乱数列とよべないものが並ぶことは滅多にないと考えられますが、確率がゼロでないことは実際に発生しうることを意味します。


# Python モジュール spcrand

:::note info
RDTSC を使う方法なので疑似乱数ではありません。
:::

:::note warn
RDTSC はシステムによって無効になっている場合があります。
:::

どんな数列になるかは実行してみなければ分からないので Python モジュールにしてみました。


メソッドは１つです

```text:
spcrand() -> int:  符号なし整数(32ビット)で生成プログラム一回分
spcrand(count: int) -> bytearray:  生成プログラム count 回分で count * 4 バイトの数列
```

**動作確認している環境**

OS
　　macOS 12.2.1

C コンパイラ (Xcode 13.2.1)

　　Apple clang version 13.0.0 (clang-1300.0.29.30)
　　Target: x86_64-apple-darwin21.3.0

Python : /usr/bin/python3

　　Python 3.8.9 (default, Oct 26 2021, 07:25:54)
　　[Clang 13.0.0 (clang-1300.0.29.30)] on darwin


```c:spcrand.c
#ifndef PY_SSIZE_T_CLEAN
#define PY_SSIZE_T_CLEAN
#endif /* PY_SSIZE_T_CLEAN */

#include <Python.h>

#if (defined(i386) ||       \
     defined(__i386) ||     \
     defined(__i386__) ||   \
     defined(__i486__) ||   \
     defined(__i586__) ||   \
     defined(__i686__) ||   \
     defined(__amd64)   ||  \
     defined(__amd64__) ||  \
     defined(__x86_64) ||   \
     defined(__x86_64__) || \
     0)
#define ISA_X86  1
#else
#define ISA_X86  0
#endif

/*
 * Method: spcrand
 */

#if ISA_X86
#define HAS_CHECK_SSE42  1

typedef struct cpuid_result {
    unsigned int eax; /* reg:00 */
    unsigned int ecx; /* reg:01 */
    unsigned int edx; /* reg:10 */
    unsigned int ebx; /* reg:11 */
} cpuid_result;

inline static int
check_sse42()
{
    cpuid_result res;

    res.eax = 1;
    res.ecx = 0;
    res.edx = 0;
    res.ebx = 0;

    __asm__ volatile (
        "mov %0,%%eax;"
        "mov %1,%%ecx;"
        "xor %%edx,%%edx;"
        "xor %%ebx,%%ebx;"
        "cpuid;"
        "mov %%eax,%0;"
        "mov %%ecx,%1;"
        "mov %%edx,%2;"
        "mov %%ebx,%3;"
        : "+m" (res.eax),
          "+m" (res.ecx),
          "+m" (res.edx),
          "+m" (res.ebx)
        : /**/
        : "eax", "ecx", "edx", "ebx");

    return (res.ecx >> 20) & 1;
}

inline static uint32_t ezrand()
{
#if defined(__GNUC__)

    uint32_t eax, edx;

    __asm__ volatile (
        "rdtsc\n"
        "crc32l\t%%edx, %%eax"
        : "=a" (eax), "=d" (edx));

    return eax;

#else  /* !__GNUC__ */
#error "unknown compiler"
#endif
}

#else  /* !ISA_X86 */
#error "unknown processor"
#endif

static PyObject *
method_spcrand(PyObject *unused, PyObject *const *args, Py_ssize_t size)
{
    PyByteArrayObject *ba;
    PyObject *r = NULL;
    uint32_t *buffer;
    Py_ssize_t count;
    Py_ssize_t i;

    (void) unused;

    if (size == 0)
        return PyLong_FromUnsignedLong(ezrand());
    if (size != 1) {
        PyErr_SetString(PyExc_TypeError, "spcrand expected at 1 argument");
        return NULL;
    }
    if (!PyLong_CheckExact(args[0])) {
        PyErr_SetObject(PyExc_TypeError, args[0]);
        return NULL;
    }

    count = PyLong_AsSsize_t(args[0]);
    if (count < 0) {
        if (!PyErr_Occurred())
            PyErr_SetObject(PyExc_TypeError, args[0]);
        return NULL;
    }

    if (!(r = PyByteArray_FromStringAndSize("", 0)))
        return NULL;
    if (PyByteArray_Resize(r, count * sizeof(uint32_t)) < 0) {
        Py_XDECREF(r);
        return NULL;
    }
    ba = (PyByteArrayObject *) r;
    buffer = (uint32_t *) ba->ob_start;
    for (i = 0; i < count; i++)
        buffer[i] = ezrand();
    return r;
}

/*
 * Module: spcrand
 */

static PyMethodDef spcrand_methods[] = {
    {"spcrand", (PyCFunction) method_spcrand, METH_FASTCALL, NULL},
    {NULL, NULL, 0, NULL}, /* end */
};

static PyModuleDef spcrand_def = {
    PyModuleDef_HEAD_INIT,
    .m_name = "spcrand",
    .m_doc = NULL,
    .m_size = -1,
    .m_methods = spcrand_methods,
};

PyMODINIT_FUNC
PyInit_spcrand(void)
{
#if HAS_CHECK_SSE42
    if (check_sse42() < 0)
        return NULL;
#endif /* HAS_CHECK_SSE42 */
    return PyModule_Create(&spcrand_def);
}
```

```python:setup.py
#!/usr/bin/env python3

import os
from distutils.core import setup, Extension

def getenv(name, defval=None):
    if name in os.environ:
        return os.environ[name]
    return defval

DEBUG = getenv('DEBUG') in ('true', 'yes')

MAJOR_VERSION = 0
MINOR_VERSION = 1
DEBUG_VERSION = 0
VERSION = '%d.%d.%d' % (MAJOR_VERSION, MINOR_VERSION, DEBUG_VERSION)

DEFINE_MACROS = [
    ('MAJOR_VERSION', MAJOR_VERSION),
    ('MINOR_VERSION', MINOR_VERSION),
    ('DEBUG_VERSION', DEBUG_VERSION),
]
UNDEF_MACROS = []

EXTRA_COMPILE_ARGS = [
    '-W',
    '-Wall',
    '-Wno-invalid-offsetof',
    '-Wno-deprecated-declarations',
]

if DEBUG:
    DEFINE_MACROS.append(('DEBUG', 1))
    UNDEF_MACROS.append('NDEBUG')
    EXTRA_COMPILE_ARGS.append('-O0')
    pass

setup(name='spcrand',
      version=VERSION,
      description='',
      ext_modules=[Extension(
          name='spcrand',
          define_macros=DEFINE_MACROS,
          undef_macros=UNDEF_MACROS,
          extra_compile_args=EXTRA_COMPILE_ARGS,
        sources=['spcrand.c'])])
```

<sub>setuptools を使えという警告が出るのはわかっているのですが…</sub>

# ビルドとテスト

```console:
$ env ARCHFLAGS="-arch x86_64" python3 setup.py build
running build
running build_ext
building 'spcrand' extension
clang -Wno-unused-result -Wsign-compare -Wunreachable-code -fno-common -dynamic -DNDEBUG -g -fwrapv -O3 -Wall -iwithsysroot/System/Library/Frameworks/System.framework/PrivateHeaders -iwithsysroot/Applications/Xcode.app/Contents/Developer/Library/Frameworks/Python3.framework/Versions/3.8/Headers -Werror=implicit-function-declaration -arch x86_64 -DMAJOR_VERSION=0 -DMINOR_VERSION=1 -DDEBUG_VERSION=0 -I/Applications/Xcode.app/Contents/Developer/Library/Frameworks/Python3.framework/Versions/3.8/include/python3.8 -c spcrand.c -o build/temp.macosx-10.14-x86_64-3.8/spcrand.o -W -Wall -Wno-invalid-offsetof -Wno-deprecated-declarations
creating build/lib.macosx-10.14-x86_64-3.8
clang -bundle -undefined dynamic_lookup -Wl,-headerpad,0x1000 -arch x86_64 build/temp.macosx-10.14-x86_64-3.8/spcrand.o -o build/lib.macosx-10.14-x86_64-3.8/spcrand.cpython-38-darwin.so
$
$ (cd build/lib.macosx-10.14-x86_64-3.8; python3)
Python 3.8.9 (default, Feb 18 2022, 07:45:34) 
[Clang 13.1.6 (clang-1316.0.21.2)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> from spcrand import spcrand
>>> hex(spcrand())
'0x1ff4677b'
>>> hex(spcrand())
'0xe7806714'
>>> hex(spcrand())
'0xf996f420'
>>> r = spcrand(4096)
>>> [r.count(i) for i in range(256)]
[72, 69, 60, 75, 61, 62, 60, 57, 59, 64, 66, 51, 58, 68, 58, 68, 63, 84, 68, 73, 61, 63, 65, 60, 66, 70, 68, 66, 55, 68, 78, 63, 60, 71, 64, 64, 49, 73, 60, 61, 79, 59, 62, 71, 52, 69, 63, 67, 68, 45, 51, 69, 57, 73, 70, 45, 72, 61, 58, 75, 55, 62, 72, 62, 55, 60, 70, 73, 63, 63, 47, 63, 61, 58, 56, 72, 69, 62, 72, 75, 68, 57, 51, 66, 59, 57, 69, 57, 67, 59, 59, 68, 68, 55, 57, 64, 82, 62, 63, 70, 59, 61, 76, 60, 62, 46, 64, 64, 72, 66, 64, 57, 69, 59, 60, 67, 64, 67, 54, 68, 68, 62, 71, 68, 67, 60, 58, 53, 78, 65, 67, 65, 57, 62, 56, 48, 56, 53, 60, 59, 69, 72, 58, 61, 51, 64, 66, 66, 65, 62, 63, 53, 59, 58, 55, 65, 73, 62, 63, 72, 64, 70, 61, 84, 54, 71, 72, 76, 67, 63, 67, 55, 59, 73, 62, 61, 56, 68, 64, 63, 72, 67, 77, 63, 66, 69, 54, 60, 64, 73, 59, 67, 67, 76, 58, 74, 66, 62, 64, 61, 57, 70, 69, 61, 64, 75, 79, 77, 59, 76, 67, 51, 67, 66, 58, 64, 57, 53, 58, 63, 68, 65, 74, 52, 69, 72, 64, 59, 65, 63, 69, 59, 61, 51, 65, 88, 57, 71, 63, 57, 62, 62, 60, 61, 67, 73, 61, 66, 69, 65, 75, 65, 75, 51, 62, 74]
>>> min([r.count(i) for i in range(256)])
45
>>> max([r.count(i) for i in range(256)])
88
>>> len(r)
16384
>>> def dump(l):
...     for s in range(0, len(l), 16):
...         print(' '.join(['%02x' % d for d in l[s:s+16]]))
... 
>>> dump(r[:256])
a2 f8 99 3c a5 4f fc ea 06 b8 db 59 b5 19 3e b1
9d fa 3c d5 59 4d 38 2f 1f 9d 7d 8c a3 0f 7e da
e6 9b 1d 9f 4f 5b 78 44 27 28 bd b3 68 8b 9c 9e
5f 0d ba 1f e6 6c d9 a6 28 ec 9f 34 1f 6a b9 b5
b6 aa dc 6e de d9 19 99 10 59 5f 0b 0f 3c 7b ee
92 3e 1e 52 2e ac 1d 04 8d 5b 3a b7 c2 f8 1b 9a
aa 8b de 6d b5 ee fa 88 09 7c f9 de 66 9d d5 37
4e 7e d7 53 79 f8 f1 d2 36 5b d0 ff f7 e8 70 d3
c0 6e 56 52 8f cd 77 7f f4 ac 56 35 48 3e 55 63
0d aa 36 26 f5 90 b4 68 49 02 b7 3e 35 1f f6 19
1d fc f4 7d 0d 5d f2 1f c3 dd b4 8d f4 5b 92 0c
0c 61 10 42 b0 f3 13 14 4d cd 35 8c e4 0d 50 57
79 0f 35 eb c5 9d 36 bd 5b db 75 e7 78 33 d7 b6
c4 a1 d4 e0 5e df d1 31 90 5f 97 a3 97 08 20 10
09 4e 63 4a b5 dc 60 1c 96 34 c2 4d 2a a6 c1 1b
31 fb a3 75 7e 58 82 58 38 88 c7 fb 0f 0e e1 7a
>>> dump(r[256:256+256])
c1 8e a7 e8 04 05 41 4f b8 97 42 19 c4 8a 03 3e
fc c8 07 38 c1 79 63 d1 84 ed 00 94 38 7f 03 c2
7c d7 82 da 41 66 e6 33 bc 58 c0 ab 7b ab e2 b7
b5 2b a4 25 82 ad 82 a4 aa 4e 80 c0 89 a6 22 91
35 34 21 c7 2a 51 05 22 7c 39 0b a9 c0 ab 08 ff
e2 7f 48 f3 5e ed 4b a5 fc 26 8e 4b c1 97 ea a2
5f d1 a9 f8 fb 5a ee 26 d3 b9 ec 42 e4 3f ca c3
2a bf 8c 51 d2 85 0e 1f e5 03 28 9e aa a0 09 b3
17 f9 2c 81 ab 6b 2f d7 6f dc 2b 2d d3 4e 28 7b
aa 57 cd 8a 1c f2 8c b4 d2 72 ca 26 90 9a c9 0e
5e 1a 8f 9c 69 9c a9 1d 10 85 4c ec 5f 26 6d c1
68 a0 4b 40 c1 60 2e 9b 48 92 af b6 f4 00 ac e0
eb 65 88 05 3b bc 08 2f f5 3c 4e bd d7 e8 0e b1
6b 7a 0d e7 cf f1 4a 39 51 b7 09 63 ed 25 0a 35
4f ee cf db e6 2e aa 00 d1 a8 8c 81 9e 0b ad ac
a6 49 a9 aa 29 65 ca f6 66 c6 eb db b6 1f 6b f1
```

このときのバイト列の分布としては 16384/256=64, min=64-19, max=64+24 の範囲に収まっているようです。


<details><summary>参考：count 回ループ部アセンブラ出力（抜粋）</summary><div>


```C:C
    for (i = 0; i < count; i++)
        buffer[i] = ezrand();
```

```asm:コンパイル結果（抜粋）
LBB1_14:
        rdtsc
        crc32l  %edx, %eax
        movl    %eax, (%rcx,%rsi,4)
        addq    $1, %rsi
        cmpq    %rsi, %r15
        jne     LBB1_14
```

</div></details>


# 使用例

「地上波アナログ時代の砂嵐風画面」の GIF 動画を作るプログラム

```python:
#!/usr/bin/env python3

from mkgif import File as GIF
from spcrand import spcrand

width = 320
height = 240
pixels = width * height

gif = GIF()
gif.append_logical_screen_descriptor(
    width, height, global_color=tuple((n, n, n) for n in range(256)))
gif.append_animaition_loop()
for n in range(25):
    print('Frame:', n)
    gif.append_graphic_control_extension(2)
    gif.append_image_descriptr(spcrand() & 0xff for _ in range((pixels)))
    pass
gif.append_trailer()

import sys
open(sys.argv[1], 'wb').write(gif)
```

mkgif は記事「[GIF 画像データを作る](https://qiita.com/ikiuo/items/7057c7b0401b4a6d22cd)」にあります。

![sandstorm-320x240.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/226677/66730282-d15d-2a90-ed6a-a2fef768399b.gif)

乱数は利用目的に適合していれば、いわゆる乱数評価方法において不適格なものでも問題ありません。画面上での演出など、簡単な乱数生成で十分な場合が多々あります。

---

:::note warn
タイトルからして、この記事は「ポンコツ」という意味を含ませています。
よって「この乱数アルゴリズムを使え」みたいなコメントは
その通り過ぎて鬱陶しいだけなので一切不要です。
:::

---

:::note
検証する時は、RDTSC+CRC32 自身に数列を決定する能力がないので、「プログラム・実行環境・実行時」がセットで決定されることに注意しましょう。実行環境には物理的要因(主に機器温度：温度上昇を抑えるためにパフォーマンスを低下させる仕組みが入っている場合)も考えられます。
:::

検証記事は見つけたら追加していきます。

https://qiita.com/taqu/items/d613983657d0731566df