help-octave
[Top][All Lists]

## Re: Mersenne Twister range in rand.oct and randn.oct

 From: Mike Miller Subject: Re: Mersenne Twister range in rand.oct and randn.oct Date: Mon, 5 Dec 2005 00:17:13 -0600 (CST)

```On Sun, 4 Dec 2005, Paul Kienzle wrote:

```
```It does.  Expressed in powers of two the expression is:

( [0,2^27-1] * 2^26 + [0,2^26-1] + 0.4 ) / 2^53

So the min is 0.4/2^53 and the max is (2^53-1 + 0.4)/2^53 = 1 - 0.6/2^53
```
```
[snip]

```
Given x^2 < 2y and the largest possible value of 2y is -2 *log(min(RANDU)), that means |randn| is at most
```
sqrt( -2*log(min(RANDU)) ) + 3.6541528853610088

Plugging in 0.4/2^53 for min(RANDU), randn is in (-12.332,12.332).
```
```

```
Thanks, Paul. Sorry that I needed it to be so explicit. I like the choices you've made. It is good that we don't have to check that rand is zero or one in cases like 1/rand() or 1/(1-rand()). The range on the normals is also very good for most uses.
```
```
If we want more of a range on randn or we want numbers closer to zero for rand, I suppose we could use some kind of trick. I designed this to produce N random uniforms that hit every value between 0 and 1 at increments of 10^-20 starting at 5 x 10^-21 and ending at 1 minus 5 x 10^-21 :
```
floor(100*rand(N,10)) * 10.^(-2*[1:10]') + 5*10^-21

```
But I wonder if that fails in reality because machine precision doesn't allow that many digits.
```
```
Using a similar idea, it is also possible to get A and B separately for uniform (0,1) random numbers of the form A x 10^-B and have any desired upper limit for B. That could allow extension of the range of the normal deviates to include much larger possible values. It would be much slower than the rand.oct function, but there might be some uses for a function that does that. Maybe not.
```
Mike

-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.

Octave's home on the web:  http://www.octave.org
How to fund new projects:  http://www.octave.org/funding.html
Subscription information:  http://www.octave.org/archive.html
-------------------------------------------------------------

```