00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #ifdef HAVE_CONFIG_H
00024 #include <config.h>
00025 #endif
00026
00027 #include "f77-fcn.h"
00028 #include "lo-error.h"
00029 #include "oct-rand.h"
00030 #include "oct-time.h"
00031
00032
00033
00034 #define uniform_dist 1
00035 #define normal_dist 2
00036
00037
00038 static int current_distribution = uniform_dist;
00039
00040
00041 static bool initialized = false;
00042
00043 extern "C"
00044 {
00045 F77_RET_T
00046 F77_FUNC (dgennor, DGENNOR) (const double&, const double&, double&);
00047
00048 F77_RET_T
00049 F77_FUNC (dgenunf, DGENUNF) (const double&, const double&, double&);
00050
00051 F77_RET_T
00052 F77_FUNC (setall, SETALL) (const int&, const int&);
00053
00054 F77_RET_T
00055 F77_FUNC (getsd, GETSD) (int&, int&);
00056
00057 F77_RET_T
00058 F77_FUNC (setsd, SETSD) (const int&, const int&);
00059
00060 F77_RET_T
00061 F77_FUNC (setcgn, SETCGN) (const int&);
00062 }
00063
00064 static int
00065 force_to_fit_range (int i, int lo, int hi)
00066 {
00067 assert (hi > lo && lo >= 0 && hi > lo);
00068
00069 i = i > 0 ? i : -i;
00070
00071 if (i < lo)
00072 i = lo;
00073 else if (i > hi)
00074 i = i % hi;
00075
00076 return i;
00077 }
00078
00079
00080
00081
00082
00083
00084 static void
00085 do_initialization (void)
00086 {
00087 octave_localtime tm;
00088
00089 int hour = tm.hour() + 1;
00090 int minute = tm.min() + 1;
00091 int second = tm.sec() + 1;
00092
00093 int s0 = tm.mday() * hour * minute * second;
00094 int s1 = hour * minute * second;
00095
00096 s0 = force_to_fit_range (s0, 1, 2147483563);
00097 s1 = force_to_fit_range (s1, 1, 2147483399);
00098
00099 F77_FUNC (setall, SETALL) (s0, s1);
00100
00101 initialized = true;
00102 }
00103
00104 static inline void
00105 maybe_initialize (void)
00106 {
00107 if (! initialized)
00108 do_initialization ();
00109 }
00110
00111 double
00112 octave_rand::seed (void)
00113 {
00114 maybe_initialize ();
00115
00116 union d2i { double d; int i[2]; };
00117 union d2i u;
00118 F77_FUNC (getsd, GETSD) (u.i[0], u.i[1]);
00119 return u.d;
00120 }
00121
00122 void
00123 octave_rand::seed (double s)
00124 {
00125 maybe_initialize ();
00126
00127 union d2i { double d; int i[2]; };
00128 union d2i u;
00129 u.d = s;
00130 int i0 = force_to_fit_range (u.i[0], 1, 2147483563);
00131 int i1 = force_to_fit_range (u.i[1], 1, 2147483399);
00132 F77_FUNC (setsd, SETSD) (i0, i1);
00133 }
00134
00135 std::string
00136 octave_rand::distribution (void)
00137 {
00138 maybe_initialize ();
00139
00140 if (current_distribution == uniform_dist)
00141 return "uniform";
00142 else if (current_distribution == normal_dist)
00143 return "normal";
00144 else
00145 {
00146 abort ();
00147 return "";
00148 }
00149 }
00150
00151 void
00152 octave_rand::distribution (const std::string& d)
00153 {
00154 if (d == "uniform")
00155 octave_rand::uniform_distribution ();
00156 else if (d == "normal")
00157 octave_rand::normal_distribution ();
00158 else
00159 (*current_liboctave_error_handler) ("rand: invalid distribution");
00160 }
00161
00162 void
00163 octave_rand::uniform_distribution (void)
00164 {
00165 maybe_initialize ();
00166
00167 current_distribution = uniform_dist;
00168
00169 F77_FUNC (setcgn, SETCGN) (uniform_dist);
00170 }
00171
00172 void
00173 octave_rand::normal_distribution (void)
00174 {
00175 maybe_initialize ();
00176
00177 current_distribution = normal_dist;
00178
00179 F77_FUNC (setcgn, SETCGN) (normal_dist);
00180 }
00181
00182 double
00183 octave_rand::scalar (void)
00184 {
00185 maybe_initialize ();
00186
00187 double retval = 0.0;
00188
00189 switch (current_distribution)
00190 {
00191 case uniform_dist:
00192 F77_FUNC (dgenunf, DGENUNF) (0.0, 1.0, retval);
00193 break;
00194
00195 case normal_dist:
00196 F77_FUNC (dgennor, DGENNOR) (0.0, 1.0, retval);
00197 break;
00198
00199 default:
00200 abort ();
00201 break;
00202 }
00203
00204 return retval;
00205 }
00206
00207 #define MAKE_RAND_MAT(mat, nr, nc, f, F) \
00208 do \
00209 { \
00210 double val; \
00211 for (volatile int j = 0; j < nc; j++) \
00212 for (volatile int i = 0; i < nr; i++) \
00213 { \
00214 OCTAVE_QUIT; \
00215 F77_FUNC (f, F) (0.0, 1.0, val); \
00216 mat(i,j) = val; \
00217 } \
00218 } \
00219 while (0)
00220
00221 Matrix
00222 octave_rand::matrix (int n, int m)
00223 {
00224 maybe_initialize ();
00225
00226 Matrix retval;
00227
00228 if (n >= 0 && m >= 0)
00229 {
00230 retval.resize (n, m);
00231
00232 if (n > 0 && m > 0)
00233 {
00234 switch (current_distribution)
00235 {
00236 case uniform_dist:
00237 MAKE_RAND_MAT (retval, n, m, dgenunf, DGENUNF);
00238 break;
00239
00240 case normal_dist:
00241 MAKE_RAND_MAT (retval, n, m, dgennor, DGENNOR);
00242 break;
00243
00244 default:
00245 abort ();
00246 break;
00247 }
00248 }
00249 }
00250 else
00251 (*current_liboctave_error_handler) ("rand: invalid negative argument");
00252
00253 return retval;
00254 }
00255
00256 #define MAKE_RAND_ND_ARRAY(mat, len, f, F) \
00257 do \
00258 { \
00259 double val; \
00260 for (volatile int i = 0; i < len; i++) \
00261 { \
00262 OCTAVE_QUIT; \
00263 F77_FUNC (f, F) (0.0, 1.0, val); \
00264 mat(i) = val; \
00265 } \
00266 } \
00267 while (0)
00268
00269 NDArray
00270 octave_rand::nd_array (const dim_vector& dims)
00271 {
00272 maybe_initialize ();
00273
00274 NDArray retval;
00275
00276 if (! dims.all_zero ())
00277 {
00278 retval.resize (dims);
00279
00280 int len = retval.length ();
00281
00282 switch (current_distribution)
00283 {
00284 case uniform_dist:
00285 MAKE_RAND_ND_ARRAY (retval, len, dgenunf, DGENUNF);
00286 break;
00287
00288 case normal_dist:
00289 MAKE_RAND_ND_ARRAY (retval, len, dgennor, DGENNOR);
00290 break;
00291
00292 default:
00293 abort ();
00294 break;
00295 }
00296 }
00297
00298 return retval;
00299 }
00300
00301 #define MAKE_RAND_ARRAY(array, n, f, F) \
00302 do \
00303 { \
00304 double val; \
00305 for (volatile int i = 0; i < n; i++) \
00306 { \
00307 OCTAVE_QUIT; \
00308 F77_FUNC (f, F) (0.0, 1.0, val); \
00309 array(i) = val; \
00310 } \
00311 } \
00312 while (0)
00313
00314 Array<double>
00315 octave_rand::vector (int n)
00316 {
00317 maybe_initialize ();
00318
00319 Array<double> retval;
00320
00321 if (n > 0)
00322 {
00323 retval.resize (n);
00324
00325 switch (current_distribution)
00326 {
00327 case uniform_dist:
00328 MAKE_RAND_ARRAY (retval, n, dgenunf, DGENUNF);
00329 break;
00330
00331 case normal_dist:
00332 MAKE_RAND_ARRAY (retval, n, dgennor, DGENNOR);
00333 break;
00334
00335 default:
00336 abort ();
00337 break;
00338 }
00339 }
00340 else if (n < 0)
00341 (*current_liboctave_error_handler) ("rand: invalid negative argument");
00342
00343 return retval;
00344 }
00345
00346
00347
00348
00349
00350