bugGNU Octave - Bugs: bug #58800, BIST for rng sometimes fails


bug #58800: BIST for rng sometimes fails

Submitter:  Hg200 <hg200>
Submitted:  Tue 21 Jul 2020 03:56:01 PM UTC
Category:  Test Suite Severity:  3 - Normal
Priority:  5 - Normal Item Group:  None
Status:  Fixed Assigned to:  None
Originator Name:  hg200 Open/Closed:  * Closed
Release:  * dev Operating System:  * Microsoft Windows
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment (Rich Markup)

( Jump to the original submission )

Mon 27 Jul 2020 12:02:15 PM UTC, comment #29: 

I pushed a slightly modified version of your documentation update here:

Closing as fixed.

Markus Mützel <mmuetzel>
Project Member
Sat 25 Jul 2020 07:30:20 PM UTC, comment #28: 

Markus, thank you. I have tested your cset on Windows and Linux and the generator initialization is fine. However, there is one location in the Octave-Manual that should be fix. Please find attached a proposal.

Sorry if i have caused confusion but i want to know.

From my side the report can be closed.

(file #49545)

Hg200 <hg200>
Sat 25 Jul 2020 10:04:01 AM UTC, comment #27: 

Thank you to everyone involved for giving valuable feedback.
I pushed the last changeset (minus two errors) here:

Marking as ready for test.

Markus Mützel <mmuetzel>
Project Member
Fri 24 Jul 2020 11:58:51 AM UTC, comment #26: 

Sorry, we crossed our comments. Yes, the init_key is distributed over all state vector elements, this is perfectly okay.

Michael Leitner <mleitner>
Fri 24 Jul 2020 11:57:44 AM UTC, comment #25: 

I just read now the code how this entropy vector is actually used, and I am very confident that it seems that it indeed mixes entropy[i] into all entries of the state vector higher than i. So I am quite sure that always all the entropy is used, and whatever the sequence of entropy sources in the entropy vector, provided at least one bit is different, the stream should be unrelated. Yes, I would think that already with current time, CPU time, fractional part of current time and PID this should be okay for all practical purposes.

Michael Leitner <mleitner>
Fri 24 Jul 2020 11:57:26 AM UTC, comment #24: 

Sorry, I can't follow. Are you proposing there needs to change something else in the algorithm?
IIUC, init_mersenne_twister (const uint32_t *init_key, const int key_length) already takes care of distributing the init_key over all state vector elements. That function is used for years now.

Markus Mützel <mmuetzel>
Project Member
Fri 24 Jul 2020 11:40:50 AM UTC, comment #23: 

If we always re-initialized with the same state, the sequence of numbers would always be the same. As we have established, this is not what we want. (Always re-initializing the random generator with the same state basically corresponds to zero entropy.)

Later on we raised the requirements and wanted to have a "non-predictable" sequence of random numbers. Even if that might not be easy, the PID and clock at the time of initialization of the initial state can be guesstimated.
That is why we add some additional "entropy" from random_device.
While it might be that we add a little bit less entropy on Linux after the change, 1024 bit should be more than enough if I correctly got Michael Leitner's point.
The new approach is definitely better on Windows where we didn't add any additional "entropy" before.

Markus Mützel <mmuetzel>
Project Member
Fri 24 Jul 2020 11:36:47 AM UTC, comment #22: 

"Reliable random patterns" in the sense of k-equidistribution and so on should never be a problem -- there is quite a number of transformations done on the entries of the state vector before the result is returned. So there should be no problems as to how it is seeded, provided the seed is different. Yes, even if you use a constant seed of 0 (which is what I do always if the Monte Carlo simulations are self-averaging, using random seeds only if I really need different, independent streams) the numbers should be as pseudo-random as they can be.

As to a close Hamming distance of two state vectors: as I read the wikipedia entry, x_k, x_k+1 and x_k+m (with m in the C++ reference implementation being 397) determine x_k+n. So if you put all your entropy into x_0 and initialize all the other entries of the state vector with (say) 0, you would still carry over some entropy (around 32/2 bits) into x_n+1 and thus x_n+2 and so on. In the proposed initialization from a single seed value, its entropy is propagated into all entries of the initial state vector. Putting your entropy at the end of the state vector  would indeed be a bad thing to do. Something that should be okay is to first initialize the state vector with a constant seed (but  of course doing the initialization bit shuffling and so on, that is, not having 625 times the same entries!), then XOR the clock time into the first entry of the state vector, the PID into the second, and all other sources of entropy you can think of in default of /dev/urandom in the following ones, and then draw a corresponding number of pseudo-random numbers from the generator and discard the results. This should mix all your entropy into the state vector, so that we should then really have a very large number of possible streams (and the initial bug that for equal clock time the results are equal would affect only the first pseudo-random number, which is discarded anyway).

Michael Leitner <mleitner>
Fri 24 Jul 2020 11:04:19 AM UTC, comment #21: 

Correction: Earlier we had urandom or (if not available) the clock. Now we have a mix between clock + some other stuff and random numbers.

Hg200 <hg200>
Fri 24 Jul 2020 10:54:47 AM UTC, comment #20: 

It is difficult to find literature about correctly MT-seeding and then we also have the Octave specific implementation.

My nightmare is that if we setup the generator with too few non-zero entropy numbers, that the generator will take several cycles until it produces reliable random patterns. I use Octave for Monte Carlo simulations and I don't want any degradation in the random number quality. On closer investigation of randmtzig.cc, it looks as if the MT-state is set up with the seed 19650218UL. Entropy is added later. So my understanding is that we still should get reliable random numbers even if we don't add entropy at all. Therefore, there should be no loss of statistical number quality if we remove "urandom" and use, for example, only the "system clock" instead.

The question remains whether the order or the position at which the entropy is given matters. In the patch, clock and random numbers are exchanged compared to earlier. I don't believe it has any impact, but i don't know fur sure.

Finally I agree to Markus' argument that the number of possible "entropy-realizations" determines the number of different "random-number realizatins". However: ".. if the seeds are different, the generated pseudo-random number streams will have no detectable relation ..". Here i would ask the question if this also holds for the hamming distance of the output if two different seeds come close together (but i don't have an answer for this as well).

Hg200 <hg200>
Fri 24 Jul 2020 09:37:56 AM UTC, comment #19: 

The attached patch tries to add 1024 bit of entropy from random_device. Together with the fact that the clock is also used as a seed, that is probably "secure enough" for the foreseeable future.

Rik is right: If not enough entropy is available that might throw an exception according to the standard. But that is implementation specific. It's not clear to me from the documentation I found if that can happen on Linux or Windows.
In that case, the available entropy (in blocks of 32 bit) would be added to the initial state vector.

Because I forgot earlier: I don't think it is possible to use a pseudo RNG to generate more entropy.

(file #49543)

Markus Mützel <mmuetzel>
Project Member
Fri 24 Jul 2020 08:25:28 AM UTC, comment #18: 

The Mersenne Twister has a practically infinite period (in the context of computing power available today or in the foreseeable future) of 2^19937−1, and it of course has a corresponding number of states. That means, if the seeds are different, the generated pseudo-random number streams will have no detectable relation. However, if we would only have 2^16 different possible seeds, we would have only 2^16 different streams. That's too low, you can easily enumerate that. Even 2^32 is too small. I think that we should have at least 2^64 different seeds.

And further, if it was only for linux with /dev/urandom, I see nothing that would speak against directly using that: when I request 4*625 bytes from /dev/urandom, /proc/sys/kernel/random/entropy_avail shows increasing rather than decreasing content of true entropy (which, if it fell to zero, would make /dev/random block). Of course, I do that from the command line, and the exact time when I do the call constitutes a source of entropy, while it could be that calling that in a program is considered as deterministic and does not increase entropy. But in any case, I would say /dev/urandom drains the available entropy only very slowly. However, I do not know about other systems.

Michael Leitner <mleitner>
Fri 24 Jul 2020 07:04:39 AM UTC, comment #17: 

Thank you for the review.
I agree that it is not necessary to pass the lower bound to the uniform_int_distribution constructor.

I also agree that it is enough to add only one element from the non-deterministic source to the initialization vector to ensure that the seed changes.
While Octave's random number generator is probably not meant for encryption purposes, imho we should nevertheless attempt to make it more "secure" if it is easily possible.
The result of random_devices ()-operator is an unsigned int. That might be just 16 bit worth of "entropy" depending on the implementation.
I'm not an expert. But that feels like it might make it easier to "guess" the state of the rng.

What about the attached change that - if possible - adds at least 32bit worth of "entropy" to the initialization vector (which isn't excessive imho)? I don't know if we should try to add a bit more.

(file #49539)

Markus Mützel <mmuetzel>
Project Member
Thu 23 Jul 2020 11:19:31 PM UTC, comment #16: 


Most of the patch looks good.  I have a question about the random number generation though.  The patch is

+    try
+      {
+        // The standard doesn't *guarantee* that random_device provides
+        // non-deterministic random numbers. So add entropy from this source
+        // last to make sure we gathered at least some entropy from the earlier
+        // sources.
+        std::random_device rd;
+        std::uniform_int_distribution<uint32_t> dist (0);
+        while (n < MT_N)
+          entropy[n++] = dist (rd);
+      }

First, the default lower bound is already 0 for the uniform_int_distribution constructor.  It might be clearer to just write

        std::uniform_int_distribution<uint32_t> dist;

Second, truly non-deterministic entropy is usually a limited resource.  The while loop, however, is harvesting (MT_N - n) x 32 bits.  MT_N is 624 while n is 4 so this is 19,840 bits.  If that amount is not available the program will likely throw an exception.

Instead, you could just add one 32-bit value worth of entropy with

if (n < MT_N)
  entropy[n++] = static_cast<uint32_t> rd ();

That is likely enough given the other seedings with 1) current time in seconds, 2) CPU time in microseconds, 3) fractional part of current time, 4) PID.

Or, if you really want fill out the remaining 620 entries then use a random number generator seeded from a single call to harvest entropy.

std::random_device rd;    // non-deterministic random data
std::mt19937 gen (rd ()); // Standard mersenne_twister_engine seeded with rd()
std::uniform_int_distribution<uint32_t> dist;
while (n < MT_N)
  entropy[n++] = dist (gen);

Rik <rik5>
Project Administrator
Thu 23 Jul 2020 11:14:19 PM UTC, comment #15: 

Thank you all for the feedback. The win64 build with Markus' patch from comment #11 has been finished and no problems with Rik's code or with "test("rng")" occur on Windows 7. I also ran the test suite with 0 fails.

I will built the win32 target overnight and do following tests tomorrow.

Hg200 <hg200>
Thu 23 Jul 2020 10:45:31 PM UTC, comment #14: 

Long term, it would be good to use industry standard, cross platform libraries for random numbers.  The Octave Projects page (https://wiki.octave.org/Projects) has this one under the Numerical topic:

Use C++11 <random> libraries for random number generation. Write link between Octave functions (rand, randi, randn, rande) and C++ API. Implement RandStream objects as Matlab does.

However, this isn't going to happen in the short term due to the size of the change and a lack of coders.

I ran the following test on Matlab

N = 100;

r = zeros (N,1);
for i = 1:N
  rng ("shuffle"), s1 = rng ();
  rng ("shuffle"), s2 = rng ();
  r(i) = ~ isequal (s1, s2);

All random number states were different so, either Matlab is using a higher quality time source or they have implemented a different method to ensure that every call to "rng ('shuffle')" results in a different state.  Octave should have equivalent behavior. 

Rik <rik5>
Project Administrator
Thu 23 Jul 2020 10:08:50 PM UTC, comment #13: 

@Hg200: Yes, all you are saying about init_mersenne_twister() is correct. The point I want to make is however that rand("state","reset") should not call init_mersenne_twister() if we cannot guarantee that it will get a new seed (different from the last invocation, however recent that was -- /dev/urandom would guarantee exactly that), but that in this case something like shuffle_mersenne_twister() should be called, which is guaranteed (in the sense of pseudo-random numbers) to change the state so that the resulting stream of random numbers has no recognizable relation. And I personally do not know, but I am quite sure that it is known in which way seeds have to differ for this to be fulfilled (or rather, it is definitely known which differences do not suffice). But of course, if we can guarantee that we have a source (on all supported platforms) for different seeds (even if we call it at the same point in time) then it would be safer to use that.

Michael Leitner <mleitner>
Thu 23 Jul 2020 06:13:59 PM UTC, comment #12: 

With Markus' patch (i mean that one from comment #10), Rik's code does not error (tested with N = 10000) and the "tst("rng")" passes wo errors.

The Matlab statement for "rng("shuffle")" is: "... initializes the generator based on the current time ...". What I want to know is the output of Matlab on Windows 7 for Rik's code in comparison to Octave. If anyone can do that, I would be interested.

@Michael:  " ... to do some fixed operation on the state (in the simplest case, just add 1)". The implementation in randmtzig.cc looks "memoryless" to me. That is, "init_mersenne_twister()" returns the same state if the same entropy is given at initialization. It is not intended to add information about the past. I would not recommend touching the implementation unless the original algorithm (i.e. written in papers) does provide such a feature.

Hg200 <hg200>
Thu 23 Jul 2020 05:41:03 PM UTC, comment #11: 

Reading on a little bit more about it, I believe it is save to replace reading from /dev/urandom with std::random_device on Linux. It should do essentially the same.
Since the standard doesn't guarantee that the random numbers generated by random_device are non-deterministic (or cryptographically secure), I moved the other sources to the front.
AFAICT, random_device is a "good" source of entropy on Linux and Windows.
Worst case on other platforms might be that it doesn't add any more entropy. But like Michael Leitner pointed out, the clock should be enough anyway.

The patch also adds the PID to the initialization vector like Michael Leitner suggested.

It looks like the changes didn't break the rng on Linux or Windows for me. It would be interesting to see if it fixes the original issue.

(file #49536)

Markus Mützel <mmuetzel>
Project Member
Thu 23 Jul 2020 01:24:12 PM UTC, comment #10: 

@Michael Leitner: Thanks for the feedback.
The C++11 <random> header seems to provide an interface to a "real" (possibly non-deterministic) random source:
IIUC, it uses /dev/urandom on Linux.
For mingw-w64, it should produce non-deterministic values (using rand_s) starting with gcc 9.2 (MXE Octave currently uses 9.3):

The attached patch adds the C++11 "entropy source" to the Mersenne-Twister initialization.

It might also be possible that that code block could be used as a replacement for the code block (just above in the same function) that reads from /dev/urandom.

If I correctly understand the documentation, the constructor of random_device might throw an exception. That isn't handled yet. The patch is only for testing.

@hg200: Could you try if you still see the issue with this patch? (I haven't tested yet. So I hope it cross-builds and doesn't fail.)

Doing something different on a reset (instead of re-initializing) also sounds like a good idea to me.

(file #49535)

Markus Mützel <mmuetzel>
Project Member
Thu 23 Jul 2020 12:53:21 PM UTC, comment #9: 

The "rng" function looks new. See bug #57564. Adding Guillaume to CC.

Hg200 <hg200>
Thu 23 Jul 2020 12:26:27 PM UTC, comment #8: 

Markus, both your suggestions won't do: with respect to this bug here, it is about initialization, not about generation random numbers. And a random number generation framework typically does not deal with generating its seed on its own from thin air, you have to provide some initial randomness.

I think that it is good to have a private implementation of random number generation in octave, for reasons of reproducibility between versions.

Yes, it is definitely important to have different seeds if one explicitly calls for a re-initialization of the random number generator. But if the generator has any worth, this should not be a problem: as I see it, it should suffice to initialize the state when octave is started with the clock time, and then, whenever a reset is called, to do some fixed operation on the state (in the simplest case, just add 1). As I said, if the generator is good, this should give you a new stream that has no recognizable relation to the previous one. If you initialize at startup with the clock time plus what would be the PID in linux, you would have different streams even if you started a number of octave processes in parallel.

Michael Leitner <mleitner>
Thu 23 Jul 2020 11:41:13 AM UTC, comment #7: 

Or probably better and cross-platform "uniform_int_distribution" from C++11: https://en.cppreference.com/w/cpp/numeric/random/uniform_int_distribution

Maybe the complete implementation of the Mersenne twister algorithm could be replace with the standard implementation from C++11?

Markus Mützel <mmuetzel>
Project Member
Thu 23 Jul 2020 10:51:00 AM UTC, comment #6: 

Maybe we could use the WinAPI to get entropy on Windows (which doesn't have /dev/urandom).
That way we wouldn't have to rely on the clock (in function init_mersenne_twister in liboctave/numeric/randmtzig.cc).

Possible sources could be:
https://docs.microsoft.com/en-us/windows/win32/api/wincrypt/nf-wincrypt-cryptgenrandom (Available since Windows XP. But now deprecated.)
https://docs.microsoft.com/en-us/windows/win32/api/bcrypt/nf-bcrypt-bcryptgenrandom (Current interface. But only available since Windows Vista.)

Markus Mützel <mmuetzel>
Project Member
Thu 23 Jul 2020 08:10:17 AM UTC, comment #5: 

error: assert (any (s1 != s2)) failed
error: called from
    assert at line 99 column 11
    tst_rand_reset at line 6 column 3

and here is the "distribution" of "any (s1 != s2)" over those 100 iterations:

Columns 1 through 40:

  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0

 Columns 41 through 80:

  0  0  1  0  0  0  0  1  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  1  0  0  0  0  1  0  0  0  0  1  0  0  0

 Columns 81 through 100:

  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0

Hg200 <hg200>
Thu 23 Jul 2020 03:38:44 AM UTC, comment #4: 

This is probably due to the low precision timer on Windows.  If I remember correctly it is only good to .01 seconds.

To verify, could the original reporter try this simplified example?

N = 100;

for i = 1:N
  rand ("state", "reset"), s1 = rand ("state");
  rand ("state", "reset"), s2 = rand ("state");
  assert (any (s1 != s2));

I've attached the same code as an m-file.

(file #49532)

Rik <rik5>
Project Administrator
Wed 22 Jul 2020 08:00:30 PM UTC, comment #3: 

I can't reproduce the bug here. If I call "test rng" in a loop all 20 tests pass each time.

I'm not an expert in that field. And I don't know if it is important that the seed is different each time the rng is shuffled (or the state or seed of "rand" is reset). Maybe a short pause in the BIST is good enough.

The documentation of rand contains this:

> The state or seed of the generator can be reset to a new random value using the @qcode{"reset"} keyword.

The word "new" sounds to me like the seed should be different after consecutive calls.

If I recall correctly, Rik was working on the rng in the past.
@Rik: Could you please comment on whether it is worth diving deeper into this?

Markus Mützel <mmuetzel>
Project Member
Wed 22 Jul 2020 03:51:49 PM UTC, comment #2: 

The documentation is a bit vague about 'rng ("shuffle")' and 'rand ("state", "reset")': 'By default, the generator is initialized from /dev/urandom if it is available, otherwise from CPU time, wall clock time, and the current fraction of a second.'

In deed on my linux system 'rng ("shuffle")' forces a read from urandom in 'randmtzig.cc':

void init_mersenne_twister (void)
    FILE *urandom = std::fopen ("/dev/urandom", "rb");
    if (urandom)

I don't know what is happening on Win. Probably there is no urandom and then entropy is gathered from the clock. But i can't influence this. This is hard - C code. Only idea i have is to add a 'pause()' into the BIST.

Hg200 <hg200>
Wed 22 Jul 2020 11:55:30 AM UTC, comment #1: 

I didn't have a look at the code yet. But from what you describe, it might be that the system time is used to generate a seed for the RNG.
If that is the case, maybe it would help to use a clock that returns the time with higher precision. (I haven't looked yet if such a clock is available on 32bit Windows.)

Markus Mützel <mmuetzel>
Project Member
Tue 21 Jul 2020 03:56:01 PM UTC, original submission:  

I am testing mxe-octave hg id a65cacd05892 as win-32 and win-64 target on a windows 7 laptop with a T5500. This bug report is closely related to bug #58795 and bug #58790.


test ('rng')

sometimes passes, sometimes it fails in both 32-bit and 64-bit builds. Some further investigations show that the initialization of the random number generator sometimes returns the same seed after three consecutive rng('shuffle') calls. Inserting of a "pause(0.1)" command into the BIST as follows solves the problem.

%! state = rng ();
%! unwind_protect
%!   rng ("shuffle");
%!   rng ("shuffle", "twister");
%!   pause (0.1);
%!   s = rng ("shuffle");
%!   assert (! isequal (s, rng ("shuffle")));
%!   s = rng ("shuffle", "twister");
%!   assert (! isequal (s, rng ("shuffle", "twister")));
%! unwind_protect_cleanup
%!   rng (state);
%! end_unwind_protect


Hg200 <hg200>


(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:

Attached Files
file #49545:  patch_rand.diff added by hg200 (2KiB - text/x-patch)
file #49544:  bug58800_random_device_v4.patch added by mmuetzel (3KiB - application/octet-stream)
file #49539:  bug58800_random_device_v3.patch added by mmuetzel (3KiB - application/octet-stream)
file #49536:  bug58800_random_device_v2.patch added by mmuetzel (3KiB - application/octet-stream)
file #49535:  bug58800_random_device.patch added by mmuetzel (1023B - application/octet-stream)
file #49532:  tst_rand_reset.m added by rik5 (152B - text/x-matlab)


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -email is unavailable- added by hg200
  • -email is unavailable- added by mleitner (Posted a comment)
  • -email is unavailable- added by rik5 (Updated the item)
  • -email is unavailable- added by mmuetzel
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -email is unavailable- added by hg200 (Submitted the item)
  • -email is unavailable- added by hg200

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only project members can vote.





    Follow 16 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2020-07-27 mmuetzel StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2020-07-25 hg200 Attached File- Added patch_rand.diff, #49545
    2020-07-25 mmuetzel StatusPatch Submitted Ready For Test
    2020-07-24 mmuetzel Attached File- Added bug58800_random_device_v4.patch, #49544
    2020-07-24 mmuetzel Attached File#49543 Removed
    2020-07-24 mmuetzel Attached File- Added bug58800_random_device_v4.patch, #49543
    2020-07-24 mmuetzel Attached File- Added bug58800_random_device_v3.patch, #49539
        StatusConfirmed Patch Submitted
    2020-07-23 rik5 StatusNone Confirmed
    2020-07-23 mmuetzel Attached File- Added bug58800_random_device_v2.patch, #49536
    2020-07-23 mmuetzel Attached File- Added bug58800_random_device.patch, #49535
    2020-07-23 hg200 Carbon-Copy- Added gyom
    2020-07-23 rik5 Attached File- Added tst_rand_reset.m, #49532
    2020-07-22 mmuetzel Carbon-Copy- Added rik5
    2020-07-21 hg200 Carbon-Copy- Added mmuetzel

    Back to the top

    Powered by Savane 3.11