← Back to branch summary

~ubuntu-branches/ubuntu/intrepid/ecl/intrepid

« back to all changes in this revision

Viewing changes to src/c/num_rand.d

  • Committer: Bazaar Package Importer
  • Author(s): Albin Tonnerre
  • Date: 2008-06-20 18:00:19 UTC
  • mfrom: (1.1.4 upstream)
  • Revision ID: james.westby@ubuntu.com-20080620180019-7fbz1ln5444vtkkr
Tags: 0.9j-20080306-2ubuntu1
* Enabled unicode support. (Closes: LP #123530)
* Modify Maintainer value to match the DebianMaintainerField specification.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/c/num_rand.d
 
1
/* -*- mode: c; c-basic-offset: 8 -*- */
1
2
/*
2
3
    num_rand.c  -- Random numbers.
3
4
*/
16
17
 
17
18
#include <ecl/ecl.h>
18
19
#include <time.h>
 
20
#include <stdio.h>
 
21
#include <stdlib.h>
 
22
 
 
23
#if 0
 
24
 
 
25
/*
 
26
 * Crappy random number generator
 
27
 */
 
28
 
 
29
cl_object
 
30
init_random_state()
 
31
{
 
32
        return (cl_object)time(0);
 
33
}
 
34
 
 
35
static double
 
36
generate_double(cl_object rs)
 
37
{
 
38
        rs->random.value
 
39
        = rs->random.value
 
40
        + (rs->random.value<<2)
 
41
        + (rs->random.value<<17)
 
42
        + (rs->random.value<<27);
 
43
        rs->random.value = rs->random.value & 0xffffffff;
 
44
        return (double)(rs->random.value>>1) / (4294967296.0/2.0);
 
45
}
 
46
 
 
47
#else
 
48
 
 
49
/*
 
50
 * Mersenne-Twister random number generator
 
51
 */
 
52
 
 
53
/* Period parameters */  
 
54
#define MT_N 624
 
55
#define MT_M 397
 
56
#define MATRIX_A 0x9908b0dfUL   /* constant vector a */
 
57
#define UPPER_MASK 0x80000000UL /* most significant w-r bits */
 
58
#define LOWER_MASK 0x7fffffffUL /* least significant r bits */
 
59
 
 
60
#define ulong unsigned long
 
61
 
 
62
cl_object
 
63
init_random_state()
 
64
{
 
65
        cl_index bytes = sizeof(ulong) * (MT_N + 1);
 
66
        cl_object a = cl_alloc_simple_base_string(bytes);
 
67
        ulong *mt = (ulong*)a->base_string.self;
 
68
        int j;
 
69
#if !defined(_MSC_VER) && !defined(mingw32)
 
70
        FILE *fp = fopen("/dev/urandom","r");
 
71
        if (fp) {
 
72
                fread(mt, sizeof(*mt), MT_N, fp);
 
73
                for (j=0; j < MT_N; j++){
 
74
                        mt[j] &= 0xffffffffUL;
 
75
                }
 
76
                fclose(fp);
 
77
        } else
 
78
#endif  
 
79
        {
 
80
                /* cant get urandom, use crappy source */
 
81
                mt[0] = (rand() + time(0)) & 0xffffffffUL;
 
82
                for (j=1; j < MT_N; j++){
 
83
                        mt[j] = (1812433253UL * (mt[j-1] ^ (mt[j-1] >> 30)) + j);
 
84
                        mt[j] &= 0xffffffffUL;
 
85
                }
 
86
        }
 
87
        mt[MT_N] = MT_N+1;
 
88
        return a;
 
89
}
 
90
 
 
91
ulong
 
92
generate_int32(cl_object state)
 
93
{
 
94
        static ulong mag01[2]={0x0UL, MATRIX_A};
 
95
        ulong y;
 
96
        ulong *mt = (ulong*)state->base_string.self;
 
97
        if (mt[MT_N] >= MT_N){
 
98
                /* refresh data */
 
99
                int kk;
 
100
                for (kk=0; kk < (MT_N - MT_M); kk++) {
 
101
                        y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
 
102
                        mt[kk] = mt[kk + MT_M] ^ (y >> 1) ^ mag01[y & 0x1UL];
 
103
                }
 
104
                for (; kk < (MT_N - 1); kk++) {
 
105
                        y = (mt[kk] & UPPER_MASK) | (mt[kk+1] & LOWER_MASK);
 
106
                        mt[kk] = mt[kk+(MT_M-MT_N)] ^ (y >> 1) ^ mag01[y & 0x1UL];
 
107
                }
 
108
                y = (mt[MT_N-1] & UPPER_MASK) | (mt[0] & LOWER_MASK);
 
109
                mt[MT_N-1] = mt[MT_M-1] ^ (y >> 1) ^ mag01[y & 0x1UL];
 
110
                mt[MT_N] = 0;
 
111
        }
 
112
        /* get random 32 bit num */
 
113
        y = mt[mt[MT_N]++];
 
114
        /* Tempering */
 
115
        y ^= (y >> 11);
 
116
        y ^= (y << 7) & 0x9d2c5680UL;
 
117
        y ^= (y << 15) & 0xefc60000UL;
 
118
        y ^= (y >> 18);
 
119
        return y;
 
120
}
 
121
 
 
122
static double
 
123
generate_double(cl_object state)
 
124
{
 
125
        return generate_int32(state) * (1.0 / 4294967296.0);
 
126
}
 
127
 
 
128
#endif
19
129
 
20
130
static cl_object
21
131
rando(cl_object x, cl_object rs)
22
132
{
23
133
        cl_object z;
24
 
        double d = (double)(rs->random.value>>1) / (4294967296.0/2.0);
 
134
        double d = generate_double(rs->random.value);
25
135
 AGAIN:
26
136
        if (!ecl_plusp(x)) {
27
137
                goto ERROR;
62
172
ecl_make_random_state(cl_object rs)
63
173
{
64
174
        cl_object z = cl_alloc_object(t_random);
65
 
 
66
175
        if (rs == Ct) {
67
 
                z->random.value = time(0);
 
176
                z->random.value = init_random_state();
68
177
        } else {
69
178
                if (Null(rs)) {
70
179
                        rs = ecl_symbol_value(@'*random-state*');
72
181
                if (type_of(rs) != t_random) {
73
182
                        FEwrong_type_argument(@'random-state', rs);
74
183
                }
75
 
                z->random.value = rs->random.value;
 
184
                z->random.value = cl_copy_seq(rs->random.value);
76
185
        }
77
186
        return(z);
78
187
}
79
188
 
80
 
static void
81
 
advance_random_state(cl_object rs)
82
 
{
83
 
        rs->random.value
84
 
        = rs->random.value
85
 
        + (rs->random.value<<2)
86
 
        + (rs->random.value<<17)
87
 
        + (rs->random.value<<27);
88
 
        rs->random.value = rs->random.value & 0xffffffff;
89
 
}
90
 
 
91
 
 
92
189
@(defun random (x &optional (rs ecl_symbol_value(@'*random-state*')))
93
190
@
94
191
        rs = ecl_check_cl_type(@'random', rs, t_random);
95
 
        advance_random_state(rs);
96
192
        @(return rando(x, rs));
97
193
@)
98
194