0
Not a bug

randn() precision limited to 55 digits

dattorro 2 weeks ago updated by Michal Kvasnicka 39 minutes ago 13

%randn ignores high precision

mp.Digits(1034)

t=randn(1,'mp')
t =0.1048747160164939795645722142580780200660228729248046875

t - 0.1048747160164939795645722142580780200660228729248046875
ans = 0

MATLAB Version 9.14 (R2023a)
Advanpix Multiprecision Computing Toolbox Version 5.2.8.15537

Answer

Answer
Not a bug

The rand(...,'mp') generates numbers in double precision and then converts them to mp. Hence digits are not filled up to the 1034th digit. We do this in order to produce the same random numbers as built-in MATLAB's rand function (for reproducibility).

If you want random numbers to the full precision (all 1034 digits are filled), please use mprand/mprandn.

>> mp.Digits(100);
>> randn(1,'mp')
ans =
0.8621733203681205548463140075909905135631561279296875

>> mprandn(1,'mp')
ans =
0.1868724737876430395187949415904261118721051078748075419355524734700877850288747581541290481709356331
Answer
Not a bug

The rand(...,'mp') generates numbers in double precision and then converts them to mp. Hence digits are not filled up to the 1034th digit. We do this in order to produce the same random numbers as built-in MATLAB's rand function (for reproducibility).

If you want random numbers to the full precision (all 1034 digits are filled), please use mprand/mprandn.

>> mp.Digits(100);
>> randn(1,'mp')
ans =
0.8621733203681205548463140075909905135631561279296875

>> mprandn(1,'mp')
ans =
0.1868724737876430395187949415904261118721051078748075419355524734700877850288747581541290481709356331

mprand() and mprandn() appear nowhere in the documentation nor on the Advanpix site;

Google search is empty: mprandn site:advanpix.com

Syntactically, function calls to mprandn(n) are expected equivalent to randn(n,'mp').

Sure looks like a bug to me.

Agree, this is lack of documentation. The mprand(n) was an experimental name, I just left it for the time being.

Please suggest other name, unambiguous enough. Thank you.

I just found that the mprand is not connected with rng, so there is no way how to control random number generator mprand(n)? How to setup user defined seed for mprand(n)?

From my point of view, the following two commands are now fully equivalent: mp(rand(1)) and rand(1,'mp') and final random number is effectively the same, and in double precision.

So I think, that rand(1,'mp') should represent random number in full precision defined by mp.Digits(100) by default, and be equivalent to your experimental function mprand(n).

But in the 'mp' (full precision) case should be possible to control random generator by rng command, too.


So, the suggested way is elimination of mprand(n) at all, and use the rand(1,'mp') at full precision
regime.

+1

Thank you very much for sharing your opinion. Yes, this is the way. This only difficulty is how to monitor when user calls "rng" to set the new seed. The easiest way would be to overload the 'rng' and make it part of toolbox. Or keep track of the 'rng' state internally in toolbox.... In any case, will include this into next toolbox version.

I have just released new version with this feature - 5.3.4.15817. Please update your toolbox.

The "rng" still not working with rand(1,'mp'):

>> mp.Digits(100)

ans =

300

>> rng(1)
rmp = rand(1,'mp')
rng(1)
rmp = rand(1,'mp')

rmp =

0.1019434353836572895054626158415646843010934355718609316510899011279499964199600542084573004435709827


rmp =

0.7686462842760086723268602992488197377426292102713405281965513370345202336516897460693092237330618525

Good catch!. It is not working only if the seed is the same. Now toolbox is monitoring the case when seed is changing. This is kind of light approach - without deep intrusions to MATLAB core.

To reset the sequence for the same seed would require full-blown hook over 'rng' and actually monitoring its calls. This will make it slow(er). 

Especially for me, it is important to have full control over random number generation.

Could you add some clear and simple example how "rng" works with rand(x,'mp') now? Or better, could you describe this specific behavior at function reference WWW page?

And finally, what about rng('default') and rng('shuffle'), these two options are extremely useful? 

Current situation. If user sets new seed with 'rng' (using rng(seed), rng('default') or rng('shuffle')) - then toolbox resets the state of its random number generator to this new seed. The 'new seed' means different seed from what was previously set.

In your example, first call rng(1) resets the state of generator. The second call rng(1) doesn't trigger reset because the seed is the same. This is the bug and needs to be fixed (working on this). 

  

Here is generator reset when seeds are different:

>> rng(0);
>> randn(2,1,'mp')
ans = 
        0.1868724736939647133172833335569502
        0.8829593750178176873331742398432483
 
>> rng(1);
>> randn(2,1,'mp')
ans = 
        -0.06238805020702318459649989108952415
        -0.51254018310017212737352467081772230
 
>> rng(0);
>> randn(2,1,'mp')
ans = 
        0.1868724736939647133172833335569502    
        0.8829593750178176873331742398432483    
 
>> rng(1);
>> randn(2,1,'mp')
ans = 
        -0.06238805020702318459649989108952415
        -0.51254018310017212737352467081772230

In the next version, the toolbox will work with the rng the same as built-in rand/randn (call to rng with the same seed will trigger the generator reset).

Well, here is update after few hours of trying to monitor/overload 'rng' function calls in MATLAB.

As it turns out, there is no such a way :(


 If, say, 'rng' is overloaded, then it is impossible to call the native MATLAB's function 'rng' anymore (from the overloaded). Which basically prevents us from monitoring each call to 'rng' = impossible to track each call to 'rng' to handle the case when seed is the same.

So there are two ways possible (please decide what you like better):

1. Same as now. Toolbox relies on the seed from 'rng', but resets state only iff the seed is different (see previous comment for details).

2. New, special function is introduced to reset the state of random number generator in toolbox. Say mp.ResetRandom.
Then it can be used as: 

mp.ResetRandom(seed)
mp.ResetRandom('default')
mp.ResetRandom('shuffle')

So that user needs to use mp.ResetRandom if he/she wants to reset the full-precision rand/randn(...,'mp').

Let me know what you prefer.

I think mp.ResetRandom is the way to go, if we want a full control.

This is a hard question, because for a regular MCT user like me, it is not clear how difficult it is to create fully functional mp.ResetRandom function as an full-precision alternative to the "rng" function. Moreover, full-precision random generators (rand(n) + rng ... 'mp' functions) are so specific, that it is not clear how important these functions really are for other MCT users in general.


From my point of view, is the new "mp.ResetRandom" function optimal solution, I just prefer better name: like "mp.rng" (if is it possible, of course, due to the mentioned overloading problem), to be close to standard MATLAB as possible. I am working now on the full-precision random sampling, where is this topic crucial, including the 'mp' alternative of "rng" to get full control over random generation process.