rand(3)

assiss

昨天写了个小程序，用rand(3)模拟随机过程，结果概率不对。
后来换random(3)，结果才正确。

今天一看man 3 rand，

These interfaces are obsoleted by random(3).

系统是BSD的同学可以试一下下面的程序（我看了FreeBSd的rand.c，和NetBSD的算法一样），用rand和random分别看一下结果：

1. 
   
2. #include <stdio.h>
   
3. #include <time.h>
   
4. 
   
5. #define SEAT_COUNT 100
   
6. 
   
7. int seats[SEAT_COUNT];
   
8. 
   
9. void clear_seats() {
   
10.         int i;
    
11.         for(i = 0; i < SEAT_COUNT; i++) seats[i] = -1;
    
12. }
    
13. 
    
14. int get_rand_seat() {
    
15.         int i;
    
16.         int ss[SEAT_COUNT];
    
17.         int count = 0;
    
18.         for(i = 0; i < SEAT_COUNT; i++) {
    
19.                 if(seats[i] == -1) ss[count++] = i;
    
20.         }
    
21.         return ss[rand() % count];
    
22. }
    
23. 
    
24. int main() {
    
25.         int i;
    
26.         int j;
    
27.         int right = 0, wrong = 0;
    
28. 
    
29.         srand(time(NULL));
    
30. 
    
31.         for(j = 0; j < 100000; j++) {
    
32.                 clear_seats();
    
33.                 seats[rand() % SEAT_COUNT] = 0; //被第一位抢占
    
34.                 for(i = 1; i < SEAT_COUNT; i++) {
    
35.                         if(seats[i] == -1) { //为空
    
36.                                 seats[i] = i;
    
37.                         }
    
38.                         else { //被抢
    
39.                                 seats[get_rand_seat()] = i;
    
40.                         }
    
41.                 }
    
42.                 if(seats[SEAT_COUNT - 1] == SEAT_COUNT - 1) right++;
    
43.                 else wrong++;
    
44.         }
    
45.         printf("正确的次数: %d\n", right);
    
46.         printf("错误的次数: %d\n", wrong);
    
47.         printf("概率: %f\n", (float) right / 100000.0);
    
48.         return 0;
    
49. }
    

复制代码

gvim

光换rand -> random就可以吗？
还没仔细看代码，不过和你的感觉相反，换成random之后，概率就固定为某个值了。

assiss

原帖由 gvim 于 2006-8-12 10:17 发表
光换rand -> random就可以吗？
还没仔细看代码，不过和你的感觉相反，换成random之后，概率就固定为某个值了。

当然对应的srand(3)要换成srandom(3)，呵呵。

assiss

1. 
   
2. /*
   
3. /usr/src/lib/libc/stdlib/rand.c
   
4. */
   
5. 
   
6. #include <sys/types.h>
   
7. #include <stdlib.h>
   
8. 
   
9. static u_long next = 1;
   
10. 
    
11. int
    
12. rand()
    
13. {
    
14.         /* LINTED integer overflow */
    
15.         return (int)((next = next * 1103515245 + 12345) % ((u_long)RAND_MAX + 1));
    
16. /*
    
17. 错误似乎出在这里，如果不+1,结果就正确。
    
18. 另外很奇怪为什么要用%，这造成rand速度极慢。
    
19. */
    
20. }
    
21. 
    
22. void
    
23. srand(seed)
    
24.         u_int seed;
    
25. {
    
26.         next = seed;
    
27. }
    

复制代码

gvim

原帖由 assiss 于 2006-8-12 10:21 发表

当然对应的srand(3)要换成srandom(3)，呵呵。

再看看

gvim

sys/lib/libkern/arch/i386/random.S

/* Here is a very good random number generator.  This implementation is
* based on ``Two Fast Implementations of the "Minimal Standard" Random
* Number Generator'', David G. Carta, Communications of the ACM, Jan 1990,
* Vol 33 No 1.  Do NOT modify this code unless you have a very thorough
* understanding of the algorithm.  It's trickier than you think.  If
* you do change it, make sure that its 10,000'th invocation returns
* 1043618065.
*
* Here is easier-to-decipher pseudocode:
*
* p = (16807*seed)<30:0>       # e.g., the low 31 bits of the product
* q = (16807*seed)<62:31>      # e.g., the high 31 bits starting at bit 32
* if (p + q < 2^31)
*      seed = p + q
* else
*      seed = ((p + q) & (2^31 - 1)) + 1
* return (seed);
*
* The result is in (0,2^31), e.g., it's always positive.
*/
#include <machine/asm.h>

        .data
randseed:
        .long   1
        .text
ENTRY(random)
        movl    $16807,%eax
        imull   randseed
        shld    $1,%eax,%edx
        andl    $0x7fffffff,%eax
        addl    %edx,%eax
        js      1f
        movl    %eax,randseed
        ret
1:
        subl    $0x7fffffff,%eax
        movl    %eax,randseed
        ret

数学不够，看不懂......
>>>It's trickier than you think

assiss

原帖由 gvim 于 2006-8-12 10:40 发表
sys/lib/libkern/arch/i386/random.S

数学不够，看不懂......
>>>It's trickier than you think

我昨天为这个问题，写了大约10个程序，分别用数组、链表、Hash等数据结构，用C/C++/PYTHON等语言，用LIBC/GLIB等库，

结果表明：
PYTHON的几个程序都正确。
用GLIB中的rand函数的程序都正确。
用libc中的rand函数的程序结果都错。

这也可能是NetBSD手册里建议使用random(3)代替rand(3)的原因吧。