第22章 共用体(union)类型

22.1 伪随机数生成程序

我们可以通过不同的方法生成一个介于0~1之间的随机浮点数。最简单的做法是:使用Mersenne twister(马特赛特旋转演算法)之类的PRNG[1]生成32位DWORD值,把这个值转换为单精度float之后再除以RAND_MAX(本例中是0xFFFFFFFF)。这样就可以得到介于0~1之间的单精度浮点数。

但是除法的运算速度非常慢。而且从效率的角度考虑,随机函数同样应当尽量少用FPU的操作指令。那么,我们是否可以彻底脱离除法运算实现随机函数呢?

单精度浮点数(float)型数据由符号位、有效数字和指数三部分构成。这种数据结构决定,只要随机填充有效数字位就可以生成随机的单精度数!

依据有关规范,随机浮点数的指数部分不可以是0。那么我们不妨把指数部分直接设置为01111111(即十进制的1),用随机数字填充有效数字部分,再把符号位设置为零(表示正数)。瞧!像模像样!这就可以生成一个介于1~2之间的随机数。再把它减去1就可以得到一个标准的随机函数。

本例利用了线性同余随机数生成随机数的算法[2]。它使用UNIX格式的系统时间作为PRNG的随机种子,继而生成32位数字。

在采用这种算法时,我们可采用共用体类型(union)的数据表示单精度浮点数(float)。这个方法可利用C/C++的数据结构,按照一种与读取方式不同的数据类型存储浮点数据的各个组成部分。本例中,我们创建一个union型变量,然后把它当作float型或uint32_t型数据进行读取。换句话说,这是一种hack,而且还是比较深度的hack。

第20章的例子介绍过整数型PRNG的有关程序。因而本节不再复述其数据格式。

include <stdio.h>
#include <stdint.h>
#include <time.h>

// integer PRNG definitions, data and routines:

//constants from the Numerical Recipes book
const uint32_t RNG_a=1664525;
const uint32_t RNG_c=1013904223;
uint32_t RNG_state; // global variable

void my_srand(uint32_t i)
{
         RNG_state=i;
};
uint32_t my_rand()
{
         RNG_state=RNG_state*RNG_a+RNG_c;
         return RNG_state;
};

// FPU PRNG definitions and routines:

union uint32_t_float
{
         uint32_t i;
         float f; 
};

float float_rand()
{
         union uint32_t_float tmp;
         tmp.i=my_rand() & 0x007fffff | 0x3F800000;
         return tmp.f-1;
};

// test

int main() 
{
         my_srand(time(NULL)); // PRNG initialization

         for (int i=0; i<100; i++)
                  printf ("%f\n", float_rand());

         return 0;
};
22.1.1 x86

使用MSVC 2010(启用选项/Ox)编译上述程序,可得到如下所示的代码。

指令清单22.1 Optimizing MSVC 2010

$SG4232  DB     '%f', 0aH, 00H

__real@3ff0000000000000 DQ 03ff0000000000000r ;1

tv140 = -4                
_tmp$ = -4                
?float_rand@@YAMXZ PROC
        push    ecx
        call    ?my_rand@@YAIXZ
; EAX=pseudorandom value
        and     eax, 8388607        ; 007fffffH
        or      eax, 1065353216     ; 3f800000H
; EAX=pseudorandom value & 0x007fffff | 0x3f800000
; store it into local stack:
        mov     DWORD PTR _tmp$[esp+4], eax
; reload it as float point number:
        fld     DWORD PTR _tmp$[esp+4]
; subtract 1.0:
        fsub    QWORD PTR __real@3ff0000000000000
; store value we got into local stack and reload it:
        fstp    DWORD PTR tv130[esp+4] ; \  these instructions are redundant
        fld     DWORD PTR tv130[esp+4] ; /
        pop     ecx
        ret     0
?float_rand@@YAMXZ ENDP

_main   PROC
        push    esi
        xor     eax, eax
        call    _time
        push    eax
        call    ?my_srand@@YAXI@Z
        add     esp, 4
        mov     esi, 100
$LL3@main:
        call    ?float_rand@@YAMXZ
        sub     esp, 8
        fstp    QWORD PTR [esp]
        push    OFFSET $SG4238
        call    _printf
        add     esp, 12
        dec     esi
        jne     SHORT $LL3@main
        xor     eax, eax
        pop     esi
        ret     0
_main   ENDP

在使用C++编译器编译之后,可执行程序中的函数名称彻底走样。这是编译器对函数的命名规范所决定的。本书51.1.1节将详细介绍这种规范。

如果使用MSVC 2012编译器进行编译,那么可执行程序则会使用面向FPU的SIMD指令。本书27.5节会对SMID指令进行单独介绍。

22.1.2 MIPS

指令清单22.2 Optimizing GCC 4.4.5

float_rand:

var_10          = -0x10
var_4           = -4
                lui       $gp, (__gnu_local_gp >> 16)
                addiu     $sp, -0x20
                la        $gp, (__gnu_local_gp & 0xFFFF)
                sw        $ra, 0x20+var_4($sp)
                sw        $gp, 0x20+var_10($sp)
; call my_rand():
                jal       my_rand
                or        $at, $zero ; branch delay slot, NOP
; $v0=32-bit pseudorandom value
                li        $v1, 0x7FFFFF
; $v1=0x7FFFFF
                and       $v1, $v0, $v1
; $v1=pseudorandom value & 0x7FFFFF
                lui       $a0, 0x3F80
; $a0=0x3F800000
                or        $v1, $a0
; $v1=pseudorandom value & 0x7FFFFF | 0x3F800000
; matter of the following instruction is still hard to get:
                lui       $v0, ($LC0 >> 16)
; load 1.0 into $f0:
                lwc1      $f0, $LC0
; move value from $v1 to coprocessor 1 (into register $f2)
; it behaves like bitwise copy, no conversion done:
                mtc1      $v1, $f2
                lw        $ra, 0x20+var_4($sp)
; subtract 1.0. leave result in $f0:
                sub.s     $f0, $f2, $f0
                jr        $ra
                addiu     $sp, 0x20 ; branch delay slot
main:

var_18          = -0x18
var_10          = -0x10
var_C           = -0xC
var_8           = -8
var_4           = -4
                lui       $gp, (__gnu_local_gp >> 16)
                addiu     $sp, -0x28
                la        $gp, (__gnu_local_gp & 0xFFFF)
                sw        $ra, 0x28+var_4($sp)
                sw        $s2, 0x28+var_8($sp)
                sw        $s1, 0x28+var_C($sp)
                sw        $s0, 0x28+var_10($sp)
                sw        $gp, 0x28+var_18($sp)
                lw        $t9, (time & 0xFFFF)($gp)
                or        $at, $zero ; load delay slot, NOP
                jalr      $t9
                move      $a0, $zero ; branch delay slot
                lui       $s2, ($LC1 >> 16) # "%f\n"
                move      $a0, $v0
                la        $s2, ($LC1 & 0xFFFF) # "%f\n"
                move      $s0, $zero
                jal       my_srand
                li        $s1, 0x64 # 'd' ; branch delay slot
loc_104:
                jal       float_rand
                addiu     $s0, 1
                lw        $gp, 0x28+var_18($sp)
; convert value we got from float_rand() to double type (printf() need it):
                cvt.d.s   $f2, $f0
                lw        $t9, (printf & 0xFFFF)($gp)
                mfc1      $a3, $f2
                mfc1      $a2, $f3
                jalr      $t9
                move      $a0, $s2
                bne       $s0, $s1, loc_104
                move      $v0, $zero
                lw        $ra, 0x28+var_4($sp)
                lw        $s2, 0x28+var_8($sp)
                lw        $s1, 0x28+var_C($sp)
                lw        $s0, 0x28+var_10($sp)
                jr        $ra
                addiu     $sp, 0x28 ; branch delay slot

$LC1:           .ascii "%f\n"<0>
$LC0:           .float 1.0

上述程序同样出现了无实际意义的LUI指令。17.5.6节解释过这种情况了,本节不再对其进行复述。

22.1.3 ARM (ARM mode)

指令清单22.3 Optimizing GCC 4.6.3 (IDA)

float_rand
                STMFD   SP!, {R3,LR}
                BL      my_rand
; R0=pseudorandom value
                FLDS    S0, =1.0
; S0=1.0
                BIC     R3, R0, #0xFF000000
                BIC     R3, R3, #0x800000
                ORR     R3, R3, #0x3F800000
; R3=pseudorandom value & 0x007fffff | 0x3f800000
; copy from R3 to FPU (register S15).
; it behaves like bitwise copy, no conversion done:
                FMSR    S15, R3
; subtract 1.0 and leave result in S0:
                FSUBS   S0, S15, S0
                LDMFD   SP!, {R3,PC}
flt_5C          DCFS 1.0

main
                STMFD   SP!, {R4,LR}
                MOV     R0, #0
                BL      time
                BL      my_srand
                MOV     R4, #0x64 ; 'd'
loc_78
                BL      float_rand
; S0=pseudorandom value
                LDR     R0, =aF         ; "%f"
; convert float type value into double type value (printf() will need it):
                FCVTDS  D7, S0
; bitwise copy from D7 into R2/R3 pair of registers (for printf()):
                FMRRD   R2, R3, D7
                BL      printf
                SUBS    R4, R4, #1
                BNE     loc_78
                MOV     R0, R4
                LDMFD   SP!, {R4,PC}

aF              DCB "%f",0xA,0

通过objdump工具对上述程序进行分析,可看到objdump所显示的FPU指令和IDA所显示的指令并不一致。这大概是因为IDA与binutils的研发人员使用的是不同的指令手册。不管原因如何,同一个FPU指令会有不同对应名称的情况确实存在。

指令清单22.4 Optimizing GCC 4.6.3 (objdump)

00000038 <float_rand>:
  38:   e92d4008        push    {r3, lr}
  3c:   ebfffffe        bl      10 <my_rand>
  40:   ed9f0a05        vldr    s0, [pc, #20]        ; 5c <float_rand+0x24>
  44:   e3c034ff        bic     r3, r0, #-16777216   ; 0xff000000
  48:   e3c33502        bic     r3, r3, #8388608     ; 0x800000
  4c:   e38335fe        orr     r3, r3, #1065353216  ; 0x3f800000
  50:   ee073a90        vmov    s15, r3
  54:   ee370ac0        vsub.f32        s0, s15, s0
  58:   e8bd8008        pop     {r3, pc}
  5c:   3f800000        svccc   0x00800000

00000000 <main>:
   0:   e92d4010        push    {r4, lr}
   4:   e3a00000        mov     r0, #0
   8:   ebfffffe        bl      0 <time>
   c:   ebfffffe        bl      0 <main>
   10:  e3a04064        mov     r4, #100        ; 0x64
   14:  ebfffffe        bl      38 <main+0x38>
   18:  e59f0018        ldr     r0, [pc, #24]   ; 38 <main+0x38>
   1c:  eeb77ac0        vcvt.f64.f32    d7, s0
   20:  ec532b17        vmov    r2, r3, d7
   24:  ebfffffe        bl      0 <printf>
   28:  e2544001        subs    r4, r4, #1
   2c:  1afffff8        bne     14 <main+0x14>
   30:  e1a00004        mov     r0, r4
   34:  e8bd8010        pop     {r4, pc}
   38:  00000000        andeq   r0, r0, r0

在float_rand()函数里地址5c处的指令和main()函数里地址38处的指令都是随机噪音。

22.2 计算机器精度

单精度浮点数的“机器精度/machine epsilon”指的是相对误差的上限,即FPU操作的最小值。因此,浮点数的数据位越多,误差越小、精度越高。故而单精度float型数据的最高精度为223=1.19e7,而双精度double型的最高精度为252=2.22e16

所以计算某一数值的机器精度是可能的。

#include <stdio.h>
#include <stdint.h>

union uint_float
{
        uint32_t i;
        float f; 
};

float calculate_machine_epsilon(float start)
{
        union uint_float v;
        v.f=start;
        v.i++;
        return v.f-start;
}; 
void main() 
{
        printf ("%g\n", calculate_machine_epsilon(1.0));
};

上述程序从IEEE 754格式的浮点数中提取小数部分,把这部分当作整数处理并且给它加上1。运算的中间值是“输入值+机器精度”。通过测算输入值(单精度型数值)与中间值的差值来测算具体float型数据的机器精度。

程序使用union型数据结构解析IEEE 754格式的float型浮点数,并利用这种数据结构把float数据的小数部分提取为整数。“1”实际上加到浮点数小数部分中去了。当然,这可能造成溢出,可能把最高位的进位加到浮点数的指数部分。

22.2.1 x86

指令清单22.5 Optimizing MSVC 2010

tv130 = 8
_v$ = 8
_start$ = 8
_calculate_machine_epsilon PROC
        fld      DWORD PTR _start$[esp-4]
        fst      DWORD PTR _v$[esp-4]      ;冗余代码
        inc      DWORD PTR _v$[esp-4]
        fsubr    DWORD PTR _v$[esp-4]
        fstp     DWORD PTR tv130[esp-4]    ;冗余代码
        fld      DWORD PTR tv130[esp-4]    ;冗余代码
        ret      0
_calculate_machine_epsilon ENDP

上述程序的第二个FST指令是冗余指令:没有理由把输入值在同一个地址存储两次。这是编译器的问题:它决定把变量v的存储地址和传递参数所用的栈内地址分配成同一个地址。

接下来的INC指令把输入值的小数部分当作整数数据处理。处理结果的中间值再被当作IEEE 754型数据传递给FPU。FSUBR指令计算差值,最后返回值被存储到ST0之中。

程序中最后两条FSTP/FLD指令没有实际作用。编译器没能进行相应的优化。

22.2.2 ARM64

把数据类型扩展为64位数据的程序如下所示。

#include <stdio.h>
#include <stdint.h>

typedef union
{
        uint64_t i;
        double  d;
} uint_double;

double calculate_machine_epsilon(double start)
{
        uint_double v;
        v.d=start;
        v.i++;
        return v.d-start;
}
void main() 
{
        printf ("%g\n", calculate_machine_epsilon(1.0));
};

ARM64平台的指令不能直接在FPU的D-字头寄存器里进行加法运算。因此输入值首先从D0寄存器复制到GPR,在那里进行运算并把结果复制给D1寄存器,从而在FPU里进行减法运算。

指令清单22.6 Optimizing GCC 4.9 ARM64

calculate_machine_epsilon:
        fmov     x0, d0      ; load input value of double type into X0
        add      x0, x0, 1   ; X0++
        fmov     d1, x0      ; move it to FPU register
        fsub     d0, d1, d0  ; subtract
        ret

可见上述x64程序使用到了SIMD指令。有关介绍请参见本书27.4节。

22.2.3 MIPS

面向MIPS平台的编译器使用MTC1(Move To Coprocessor 1)指令把GPR的数据传递给FPU寄存器。

指令清单22.7 Optimizing GCC 4.4.5 (IDA)

calculate_machine_epsilon:
                mfc1    $v0, $f12
                or      $at, $zero ; NOP
                addiu   $v1, $v0, 1
                mtc1    $v1, $f2
                jr      $ra
                sub.s   $f0, $f2, $f12 ; branch delay slot
22.2.4 本章小结

本章程序所使用的技巧,在实际程序中出现的几率不大。但是这些技巧可充分演示IEEE 754型数据和C/C++ UNIONS型数据结构的特点。

22.3 快速平方根计算

浮点数解释为整数的另一个著名的算法,是快速平方根计算。

指令清单22.8 取自http://go.yurichev.com/17364的源代码

/* Assumes that float is in the IEEE 754 single precision floating point format
 * and that int is 32 bits. */
float sqrt_approx(float z)
{
    int val_int = *(int*)&z; /* Same bits, but as an int */
    /*
     * To justify the following code, prove that
     *
     * ((((val_int / 2^m) - b) / 2) + b) * 2^m = ((val_int - 2^m) / 2) + ((b + 1) / 2) * 2^m)
     *
     *where
     *
     *b = exponent bias
     * m = number of mantissa bits
     *
     *
     .
     */

    val_int -= 1 << 23; /* Subtract 2^m. */
    val_int >>= 1; /* Divide by 2. */
    val_int += 1 << 29; /* Add ((b + 1) / 2) * 2^m. */

    return *(float*)&val_int; /* Interpret again as float */
}

作为一个练习,你可以尝试编译这个函数,以理解它是如何工作的。还有一个著名的快速计算\frac{1}{\sqrt{x}} 的算法。据说这一算法因为用于QuakeⅢArena中而变得很流行。

该算法的描述和源代码参见http://go.yurichev.com/17360


[1] Pseudorandom number generator,伪随机数生成函数。

[2] 参考了网络文章http://xor0110.wordpress.com/2010/09/24/how-to-generate-floating-point-random-numbers-efficiently。