第20章 线性同余法与伪随机函数

“线性同余法”大概是生成随机数的最简方法。虽然现在的随机函数基本不采用这种技术了[1],但是它原理简单(只涉及乘法、加法和与运算),仍然值得研究。

#include <stdint.h>

// constants from the Numerical Recipes book
#define RNG_a 1664525
#define RNG_c 1013904223

static uint32_t rand_state;

void my_srand (uint32_t init)
{
        rand_state=init;
}

int my_rand ()
{
        rand_state=rand_state*RNG_a;
        rand_state=rand_state+RNG_c;
        return rand_state & 0x7fff;
}

上述程序包含两个函数。第一个函数用于初始化内部状态,第二个函数生成伪随机数字。

这个程序中的两个常量来自于参考书目[Pre+07]。本文直接使用C/C++的#define语句,把它们定义为两个宏。C/C++对宏和常量的处理方法有所区别。C/C++编译器的预处理程序会把宏直接替换为相应的值。所以宏并不像变量那样占用内存。相对而言,编译器把常量当作只读变量处理。因为指针的本质是内存地址,所以C/C++的指针可以指向常量,但是无法指向宏。

函数最后使用了“&”/与操作符。根据C语言有关标准,my_rand()函数的返回值应当介于0~32767之间。如果您有意生成32位的伪随机数,那么就不必在此处进行与运算。

20.1 x86

指令清单20.1 Optimizing MSVC 2013

_BSS     SEGMENT
_rand_state DD  01H DUP (?)
_BSS     ENDS

_init$ = 8
_srand  PROC
        mov     eax, DWORD PTR _init$[esp-4]
        mov     DWORD PTR _rand_state, eax
        ret     0
_srand  ENDP

_TEXT   SEGMENT
_rand   PROC
        imul    eax, DWORD PTR _rand_state, 1664525
        add     eax, 1013904223         ; 3c6ef35fH
        mov     DWORD PTR _rand_state, eax
        and     eax, 32767              ; 00007fffH
        ret     0
_rand   ENDP

_TEXT   ENDS

编译器把源代码中的宏直接替换为宏所绑定的常量。这个例子证明。编译器不会给宏单独分配内存。my_srand()函数直接把输入值传递给内部的rand_state变量。

此后,my_rand()函数接收了这个值,以此计算rand_state。在调整了它的宽度之后,把返回值保留在EAX寄存器里。

更为详尽的计算方法可参见如下所示的非优化编译而生成的程序。

指令清单20.2 Non-optimizing MSVC 2013

_BSS    SEGMENT
_rand_state DD  01H DUP (?)
_BSS    ENDS

_init$ = 8
_srand  PROC
        push    ebp
        mov     ebp, esp
        mov     eax, DWORD PTR _init$[ebp]
        mov     DWORD PTR _rand_state, eax
        pop     ebp
        ret     0
_srand  ENDP

_TEXT   SEGMENT
_rand   PROC
        push    ebp
        mov     ebp, esp
        imul    eax, DWORD PTR _rand_state, 1664525
        mov     DWORD PTR _rand_state, eax
        mov     ecx, DWORD PTR _rand_state
        add     ecx, 1013904223         ; 3c6ef35fH
        mov     DWORD PTR _rand_state, ecx
        mov     eax, DWORD PTR _rand_state
        and     eax, 32767              ;00007 fffH
        pop     ebp
        ret     0
_rand   ENDP

_TEXT   ENDS

20.2 x64

x64的程序与x86的程序几乎相同。因为函数的返回值属于int型数据,所以它没有使用64位寄存器,而是使用了32位寄存器的助记符。但是,my_srand()函数从ECX寄存器获取的所需参数,没有通过栈读取参数,这构成了64位程序的显著特征。

指令清单20.3 Optimizing MSVC 2013 x64

_BSS     SEGMENT
rand_state DD   01H DUP (?)
_BSS     ENDS

init$ = 8
my_srand PROC
; ECX = input argument
        mov      DWORD PTR rand_state, ecx
        ret      0
my_srand ENDP

_TEXT   SEGMENT
my_rand PROC
        imul     eax, DWORD PTR rand_state, 1664525     ; 0019660dH
        add      eax, 1013904223                        ; 3c6ef35fH
        mov      DWORD PTR rand_state, eax
        and      eax, 32767                             ; 00007fffH
        ret      0
my_rand ENDP

_TEXT ENDS

GCC生成的程序与此类似。

20.3 32位ARM

指令清单20.4 Optimizing Keil 6/2013 (ARM mode)

my_srand PROC
        LDR      r1,|L0.52|  ; load pointer to rand_state
        STR      r0,[r1,#0]  ; save rand_state
        BX       lr
        ENDP

my_rand PROC
        LDR      r0,|L0.52|  ; load pointer to rand_state
        LDR      r2,|L0.56|  ; load RNG_a
        LDR      r1,[r0,#0]  ; load rand_state
        MUL      r1,r2,r1
        LDR      r2,|L0.60|  ; load RNG_c
        ADD      r1,r1,r2
        STR      r1,[r0,#0]  ; save rand_state
; AND with 0x7FFF:
        LSL      r0,r1,#17
        LSR      r0,r0,#17
        BX       lr
        ENDP

|L0.52|
        DCD      ||.data||
|L0.56|
        DCD      0x0019660d
|L0.60|
        DCD      0x3c6ef35f

        AREA ||.data||, DATA, ALIGN=2

rand_state
        DCD      0x00000000

ARM模式的单条指令本来就是32位opcode,所以单条指令不可能容纳得下32位常量。因此,Keil单独开辟了一些空间来存储常量,再通过额外的指令读取它们。

值得关注的是,常量0x7fff也无法封装在单条指令之中。Keil把rand_state左移17位之后再右移17位,这相当于C/C++程序中的“(rand_state<<17)>>17”语句。虽然看上去这是画蛇添足,但是它清除了寄存器的高17位,保留了低15位数据,与源代码的功能相符。

Keil优化编译而生成的程序与此雷同,本文不再单独介绍。

20.4 MIPS

指令清单20.5 Optimizing GCC 4.4.5 (IDA)

my_srand:
; store $a0 to rand_state:
                lui     $v0, (rand_state >> 16)
                jr      $ra
                sw      $a0, rand_state
my_rand:
; load rand_state to $v0:
                lui     $v1, (rand_state >> 16)
                lw      $v0, rand_state
                or      $at, $zero  ; load delay slot
; multiplicate rand_state in $v0 by 1664525 (RNG_a):
                sll     $a1, $v0, 2
                sll     $a0, $v0, 4
                addu    $a0, $a1, $a0
                sll     $a1, $a0, 6
                subu    $a0, $a1, $a0
                addu    $a0, $v0
                sll     $a1, $a0, 5
                addu    $a0, $a1
                sll     $a0, 3
                addu    $v0, $a0, $v0
                sll     $a0, $v0, 2
                addu    $v0, $a0
; add 1013904223 (RNG_c)
; the LI instruction is coalesced by IDA from LUI and ORI
                li      $a0, 0x3C6EF35F
                addu    $v0, $a0
; store to rand_state:
                sw      $v0, (rand_state & 0xFFFF)($v1)
                jr      $ra
                andi    $v0, 0x7FFF ; branch delay slot

在上述程序中,常量只有0x3C6EF35F (即1013904223)。另一个值为1664525的常量去哪了?

编译器通过位移和加法运算实现了“乘以1664525”的运算。我们假设编译器自行生成了以下函数:

#define RNG_a 1664525

int f (int a) 
{
         return a*RNG_a;
}

那么,上述程序的汇编指令应当如下所示。

指令清单20.6 Optimizing GCC 4.4.5 (IDA)

f:
                sll     $v1, $a0, 2
                sll     $v0, $a0, 4
                addu    $v0, $v1, $v0
                sll     $v1, $v0, 6
                subu    $v0, $v1, $v0
                addu    $v0, $a0
                sll     $v1, $v0, 5
                addu    $v0, $v1
                sll     $v0, 3
                addu    $a0, $v0, $a0
                sll     $v0, $a0, 2
                jr      $ra
                addu    $v0, $a0, $v0 ; branch delay slot

我们确实可以在可执行程序中找到对应的指令!

MIPS的重新定位

本节重点介绍寄存器与内存交换数据的具体方法。在展现MIPS程序的细节方面,IDA略有不足(“智能”整理得太严重)。因此本文使用objdump程序,分别观测到汇编层面的指令清单及重定位表(relocation list)。

指令清单20.7 Optimizing GCC 4.4.5 (objdump)

# objdump -D rand_O3.o
...

00000000 <my_srand>:
   0:   3c020000    lui     v0,0x0
   4:   03e00008    jr      ra
   8:   ac440000    sw      a0,0(v0)

0000000c <my_rand>:
   c:   3c030000    lui     v1,0x0
  10:   8c620000    lw      v0,0(v1)
  14:   00200825    move    at,at
  18:   00022880    sll     a1,v0,0x2
  1c:   00022100    sll     a0,v0,0x4
  20:   00a42021    addu    a0,a1,a0
  24:   00042980    sll     a1,a0,0x6
  28:   00a42023    subu    a0,a1,a0
  2c:   00822021    addu    a0,a0,v0
  30:   00042940    sll     a1,a0,0x5
  34:   00852021    addu    a0,a0,a1
  38:   000420c0    sll     a0,a0,0x3
  3c:   00821021    addu    v0,a0,v0
  40:   00022080    sll     a0,v0,0x2
  44:   00441021    addu    v0,v0,a0
  48:   3c043c6e    lui     a0,0x3c6e
  4c:   3484f35f    ori     a0,a0,0xf35f
  50:   00441021    addu    v0,v0,a0
  54:   ac620000    sw      v0,0(v1)
  58:   03e00008    jr      ra
  5c:   30427fff    andi    v0,v0,0x7fff

...
# objdump -r rand_O3.o
...

RELOCATION RECORDS FOR [.text]: 
OFFSET   TYPE              VALUE
00000000 R_MIPS_HI16       .bss
00000008 R_MIPS_LO16       .bss
0000000c R_MIPS_HI16       .bss
00000010 R_MIPS_LO16       .bss
00000054 R_MIPS_LO16       .bss
...

函数my_srand()两次出现了relocation现象。第一次出现在地址0处,其类型为R_MIPS_HI16。第二处出现在地址8处,类型为R_MIPS_L016。这意味着.bss段的起始地址要写到地址为0、8(分别为16位的高、低地址)的指令中去。

位于.bss段起始位置的值,是变量rand_state。

在前几条LUI和SW的指令中,某些操作符是零。因为编译器在编译阶段(complier)还不能确定这些值,所以此时把它空了出来。编译器将会在链接阶段(linker)处理这些数据,把地址的高地址位传递给LUI指令的相应操作符中,并把低地址位传递给SW指令。SW指令会对地址的低地址位和$V0寄存器的高地址位进行加权求和。

有了上述概念之后我们就不难理解my_rand()函数中的重定位:重定位标记R_MIPS_H16告诉linker程序“要把.bss段地址的高地址位传递给LUI指令”。所以变量rand_state地址的高地址位将会存储在$V1寄存器。地址为0x10的LW指令再把变量地址的低半部分读出来,然后加到$V0寄存器中。地址为0x54的SW指令再次进行这种求和,把新的数值传递给全局变量rand_state。

在读取MIPS程序的重定位信息之后,IDA程序会把数据与指令进行智能匹配。因此我们无法在IDA中看到这些原始的汇编指令。不过,无论IDA是否显示它们,重定位信息确实存在。

20.5 本例的线程安全改进版

请参见本书65.1节查询本例的线程安全(thread-safe)改进版。


[1] 多数随机函数都采用梅森旋转算法(Mersenne twister)。