|
7
|
1 // f-rand.cc -*- C++ -*- |
|
1
|
2 /* |
|
|
3 |
|
497
|
4 Copyright (C) 1993, 1994 John W. Eaton |
|
1
|
5 |
|
|
6 This file is part of Octave. |
|
|
7 |
|
|
8 Octave is free software; you can redistribute it and/or modify it |
|
|
9 under the terms of the GNU General Public License as published by the |
|
|
10 Free Software Foundation; either version 2, or (at your option) any |
|
|
11 later version. |
|
|
12 |
|
|
13 Octave is distributed in the hope that it will be useful, but WITHOUT |
|
|
14 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
|
|
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
|
|
16 for more details. |
|
|
17 |
|
|
18 You should have received a copy of the GNU General Public License |
|
|
19 along with Octave; see the file COPYING. If not, write to the Free |
|
|
20 Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. |
|
|
21 |
|
|
22 */ |
|
|
23 |
|
240
|
24 #ifdef HAVE_CONFIG_H |
|
|
25 #include "config.h" |
|
1
|
26 #endif |
|
|
27 |
|
|
28 #include "tree-const.h" |
|
|
29 #include "f77-uscore.h" |
|
|
30 #include "error.h" |
|
|
31 #include "utils.h" |
|
7
|
32 #include "f-rand.h" |
|
1
|
33 |
|
|
34 // Possible distributions of random numbers. |
|
|
35 enum rand_dist { uniform, normal }; |
|
|
36 |
|
|
37 // Current distribution of random numbers. |
|
|
38 static rand_dist current_distribution = uniform; |
|
|
39 |
|
|
40 extern "C" |
|
|
41 { |
|
|
42 int *F77_FCN (dgennor) (double*, double*, double*); |
|
|
43 int *F77_FCN (dgenunf) (double*, double*, double*); |
|
|
44 int *F77_FCN (setall) (int*, int*); |
|
|
45 int *F77_FCN (getsd) (int*, int*); |
|
|
46 } |
|
|
47 |
|
|
48 #ifdef WITH_DLD |
|
497
|
49 Octave_object |
|
506
|
50 builtin_rand_2 (const Octave_object& args, int nargout) |
|
1
|
51 { |
|
506
|
52 return rand_internal (args, nargout); |
|
1
|
53 } |
|
|
54 #endif |
|
|
55 |
|
|
56 static double |
|
|
57 curr_rand_seed (void) |
|
|
58 { |
|
|
59 union d2i { double d; int i[2]; }; |
|
|
60 union d2i u; |
|
|
61 F77_FCN (getsd) (&(u.i[0]), &(u.i[1])); |
|
|
62 return u.d; |
|
|
63 } |
|
|
64 |
|
|
65 static int |
|
|
66 force_to_fit_range (int i, int lo, int hi) |
|
|
67 { |
|
|
68 assert (hi > lo && lo >= 0 && hi > lo); |
|
|
69 |
|
|
70 i = i > 0 ? i : -i; |
|
|
71 |
|
|
72 if (i < lo) |
|
|
73 i = lo; |
|
|
74 else if (i > hi) |
|
|
75 i = i % hi; |
|
|
76 |
|
|
77 return i; |
|
|
78 } |
|
|
79 |
|
|
80 static void |
|
|
81 set_rand_seed (double val) |
|
|
82 { |
|
|
83 union d2i { double d; int i[2]; }; |
|
|
84 union d2i u; |
|
|
85 u.d = val; |
|
|
86 int i0 = force_to_fit_range (u.i[0], 1, 2147483563); |
|
|
87 int i1 = force_to_fit_range (u.i[1], 1, 2147483399); |
|
|
88 F77_FCN (setall) (&i0, &i1); |
|
|
89 } |
|
|
90 |
|
|
91 static char * |
|
|
92 curr_rand_dist (void) |
|
|
93 { |
|
|
94 if (current_distribution == uniform) |
|
|
95 return "uniform"; |
|
|
96 else if (current_distribution == normal) |
|
|
97 return "normal"; |
|
|
98 else |
|
|
99 { |
|
|
100 panic_impossible (); |
|
|
101 return (char *) NULL; |
|
|
102 } |
|
|
103 } |
|
|
104 |
|
497
|
105 Octave_object |
|
506
|
106 rand_internal (const Octave_object& args, int nargout) |
|
1
|
107 { |
|
|
108 // Assumes that we have been given the correct number of arguments. |
|
|
109 |
|
497
|
110 Octave_object retval; |
|
1
|
111 |
|
506
|
112 int nargin = args.length (); |
|
|
113 |
|
1
|
114 static int initialized = 0; |
|
|
115 if (! initialized) |
|
|
116 { |
|
|
117 // Make the random number generator give us a different sequence every |
|
|
118 // time we start octave unless we specifically set the seed. |
|
|
119 #if 0 |
|
|
120 int s0 = 1234567890; |
|
|
121 int s1 = 123456789; |
|
|
122 #else |
|
|
123 time_t now; |
|
|
124 struct tm *tm; |
|
|
125 |
|
|
126 time (&now); |
|
|
127 tm = localtime (&now); |
|
|
128 |
|
|
129 int s0 = tm->tm_min * 60 + tm->tm_sec; |
|
|
130 int s1 = (tm->tm_mday - 1) * 24 * 3600 + tm->tm_hour * 3600 + s0; |
|
|
131 #endif |
|
|
132 s0 = force_to_fit_range (s0, 1, 2147483563); |
|
|
133 s1 = force_to_fit_range (s1, 1, 2147483399); |
|
|
134 |
|
|
135 F77_FCN (setall) (&s0, &s1); |
|
|
136 initialized = 1; |
|
|
137 } |
|
|
138 |
|
|
139 int n = 0; |
|
|
140 int m = 0; |
|
|
141 if (nargin == 1) |
|
|
142 { |
|
|
143 n = 1; |
|
|
144 m = 1; |
|
|
145 goto gen_matrix; |
|
|
146 } |
|
|
147 else if (nargin == 2) |
|
|
148 { |
|
497
|
149 switch (args(1).const_type ()) |
|
1
|
150 { |
|
|
151 case tree_constant_rep::string_constant: |
|
497
|
152 char *s_arg = args(1).string_value (); |
|
1
|
153 if (strcmp (s_arg, "dist") == 0) |
|
|
154 { |
|
497
|
155 retval.resize (1); |
|
1
|
156 char *s = curr_rand_dist (); |
|
497
|
157 retval(0) = tree_constant (s); |
|
1
|
158 } |
|
|
159 else if (strcmp (s_arg, "seed") == 0) |
|
|
160 { |
|
497
|
161 retval.resize (1); |
|
1
|
162 double d = curr_rand_seed (); |
|
497
|
163 retval(0) = tree_constant (d); |
|
1
|
164 } |
|
|
165 else if (strcmp (s_arg, "uniform") == 0) |
|
|
166 current_distribution = uniform; |
|
|
167 else if (strcmp (s_arg, "normal") == 0) |
|
|
168 current_distribution = normal; |
|
|
169 else |
|
497
|
170 error ("rand: unrecognized string argument"); |
|
1
|
171 break; |
|
|
172 case tree_constant_rep::scalar_constant: |
|
|
173 case tree_constant_rep::complex_scalar_constant: |
|
497
|
174 n = NINT (args(1).double_value ()); |
|
1
|
175 m = n; |
|
|
176 goto gen_matrix; |
|
|
177 case tree_constant_rep::range_constant: |
|
|
178 { |
|
497
|
179 Range r = args(1).range_value (); |
|
1
|
180 n = 1; |
|
|
181 m = NINT (r.nelem ()); |
|
|
182 } |
|
|
183 goto gen_matrix; |
|
|
184 case tree_constant_rep::matrix_constant: |
|
|
185 case tree_constant_rep::complex_matrix_constant: |
|
497
|
186 n = NINT (args(1).rows ()); |
|
|
187 m = NINT (args(1).columns ()); |
|
1
|
188 goto gen_matrix; |
|
|
189 default: |
|
|
190 panic_impossible (); |
|
|
191 break; |
|
|
192 } |
|
|
193 } |
|
|
194 else if (nargin == 3) |
|
|
195 { |
|
497
|
196 if (args(1).is_string_type () |
|
|
197 && strcmp (args(1).string_value (), "seed") == 0) |
|
1
|
198 { |
|
497
|
199 double d = args(2).to_scalar (); |
|
1
|
200 set_rand_seed (d); |
|
|
201 } |
|
|
202 else |
|
|
203 { |
|
497
|
204 n = NINT (args(1).to_scalar ()); |
|
|
205 m = NINT (args(2).to_scalar ()); |
|
1
|
206 goto gen_matrix; |
|
|
207 } |
|
|
208 } |
|
|
209 |
|
|
210 return retval; |
|
|
211 |
|
|
212 gen_matrix: |
|
|
213 |
|
|
214 if (n == 0 || m == 0) |
|
|
215 { |
|
497
|
216 retval.resize (1); |
|
1
|
217 Matrix m (0, 0); |
|
497
|
218 retval(0) = tree_constant (m); |
|
1
|
219 } |
|
|
220 else if (n > 0 && m > 0) |
|
|
221 { |
|
497
|
222 retval.resize (1); |
|
1
|
223 Matrix rand_mat (n, m); |
|
|
224 for (int j = 0; j < m; j++) |
|
|
225 for (int i = 0; i < n; i++) |
|
|
226 { |
|
|
227 double d_zero = 0.0; |
|
|
228 double d_one = 1.0; |
|
|
229 double val; |
|
|
230 switch (current_distribution) |
|
|
231 { |
|
|
232 case uniform: |
|
|
233 F77_FCN (dgenunf) (&d_zero, &d_one, &val); |
|
|
234 rand_mat.elem (i, j) = val; |
|
|
235 break; |
|
|
236 case normal: |
|
|
237 F77_FCN (dgennor) (&d_zero, &d_one, &val); |
|
|
238 rand_mat.elem (i, j) = val; |
|
|
239 break; |
|
|
240 default: |
|
|
241 panic_impossible (); |
|
|
242 break; |
|
|
243 } |
|
|
244 } |
|
|
245 |
|
497
|
246 retval(0) = tree_constant (rand_mat); |
|
1
|
247 } |
|
|
248 else |
|
216
|
249 error ("rand: invalid negative argument"); |
|
1
|
250 |
|
|
251 return retval; |
|
|
252 } |
|
|
253 |
|
|
254 /* |
|
|
255 ;;; Local Variables: *** |
|
|
256 ;;; mode: C++ *** |
|
|
257 ;;; page-delimiter: "^/\\*" *** |
|
|
258 ;;; End: *** |
|
|
259 */ |
