我们可以通过不同的方法生成一个介于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;
};
使用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.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.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处的指令都是随机噪音。
单精度浮点数的“机器精度/machine epsilon”指的是相对误差的上限,即FPU操作的最小值。因此,浮点数的数据位越多,误差越小、精度越高。故而单精度float型数据的最高精度为2−23=1.19e−7,而双精度double型的最高精度为2−52=2.22e−16。
所以计算某一数值的机器精度是可能的。
#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.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指令没有实际作用。编译器没能进行相应的优化。
把数据类型扩展为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节。
面向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
本章程序所使用的技巧,在实际程序中出现的几率不大。但是这些技巧可充分演示IEEE 754型数据和C/C++ UNIONS型数据结构的特点。
浮点数解释为整数的另一个著名的算法,是快速平方根计算。
指令清单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 */
}
作为一个练习,你可以尝试编译这个函数,以理解它是如何工作的。还有一个著名的快速计算 的算法。据说这一算法因为用于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。