Revision 706, 1.4 kB
(checked in by smidl, 15 years ago)
|
eol-native
|
-
Property svn:eol-style set to
native
|
Rev | Line | |
---|
[579] | 1 | function r=randun(m,n) |
---|
| 2 | %RANDUN generate sample from Uniform distribution |
---|
| 3 | % |
---|
| 4 | % r = randun |
---|
| 5 | % r = randun(m) |
---|
| 6 | % r = randun(m,n) |
---|
| 7 | % |
---|
| 8 | % r : generated sample |
---|
| 9 | % SEED : global variable SEED |
---|
| 10 | % m : rows |
---|
| 11 | % n : columns |
---|
| 12 | % |
---|
| 13 | % Design : J. Andrysek |
---|
| 14 | % Project : BadDyr |
---|
| 15 | % Comments: used to generate exactly the same numbers in .m and .dll |
---|
| 16 | % versions |
---|
| 17 | % |
---|
| 18 | % Patched by Ludvik Tesar to output matrix instead of scalar |
---|
| 19 | |
---|
| 20 | |
---|
| 21 | global SEED; |
---|
| 22 | if(isempty(SEED)) SEED = 1111111; end; |
---|
| 23 | %SEED = randn('seed'); |
---|
| 24 | |
---|
| 25 | global RANDUN_COUNTER; |
---|
| 26 | if(isempty(RANDUN_COUNTER)) RANDUN_COUNTER = 0; end; |
---|
| 27 | |
---|
| 28 | |
---|
| 29 | if ~exist('m','var'); m = 1; end; |
---|
| 30 | if ~exist('n','var'); n = 1; end; |
---|
| 31 | |
---|
| 32 | r = zeros(m,n); |
---|
| 33 | |
---|
| 34 | for i=1:(m*n); |
---|
| 35 | |
---|
| 36 | % Ahrens linear congruential pseudo-random generator with A=663608941 B=0 M=2^32 |
---|
| 37 | % A = 663608941; |
---|
| 38 | % M = 4294967296; |
---|
| 39 | % SEED = soucinek(A,SEED); |
---|
| 40 | % r(i) = SEED/M; |
---|
| 41 | % replaced by: |
---|
| 42 | |
---|
| 43 | % Park-Miller linear congruential pseudo-random generator with A=16807 B=0 M=2^31-1 |
---|
| 44 | % (it spans all 2^31-1 numbers and has good statistical properties) |
---|
| 45 | A = 16807; |
---|
| 46 | M = 2147483647; |
---|
| 47 | SEED = mod(A*SEED,M); |
---|
| 48 | r(i) = SEED/M; |
---|
| 49 | end; |
---|
| 50 | |
---|
| 51 | %randn('seed',SEED); |
---|
| 52 | |
---|
| 53 | RANDUN_COUNTER = RANDUN_COUNTER+m*n; |
---|
| 54 | |
---|
| 55 | function r=soucinek(a,b); |
---|
| 56 | z=2^16; |
---|
| 57 | ha(1)=floor(a/z); |
---|
| 58 | ha(2)=mod(a,z); |
---|
| 59 | |
---|
| 60 | hb(1)=floor(b/z); |
---|
| 61 | hb(2)=mod(b,z); |
---|
| 62 | |
---|
| 63 | c=ha(2)*hb(2); |
---|
| 64 | res(2)=mod(c,z); |
---|
| 65 | res(1)=mod(floor(c/z)+hb(2)*ha(1)+hb(1)*ha(2),z); |
---|
| 66 | |
---|
| 67 | r=res(1)*z+res(2); |
---|
| 68 | |
---|