How can I generate random single precision (float32) numbers ?

A SP float may be generated in a range from -Infinity to +Infinity or including NANs.

Full rage would contain all representable numbers 0x00000000 to 0xffffffff.

like typecast(uint32(hex2dec('00000000')), 'single') to typecast(uint32(hex2dec('ffffffff')), 'single')

What is the best approach ?

How about generate without NAN, or INF, or Denormalized numbers (subnormals) ?

Tried:

size_x = 1000;
size_y = 1;

% realmin('single') 1.1755e-38
% realmax('single') 3.4028e+38
a = realmin('single') + ( realmax('single') - realmin('single') ).*rand(size_x, size_y, 'single');
a_range = [min(a) mean(a) std(a) max(a)]


b = typecast(randi([0, intmax('uint32')], size_x, size_y, 'uint32'), 'single');
b_range = [min(b) mean(b) std(b) max(b)]

[Bruno Luong]

Just curious, what you are trying to do with such randomize process?

The distribution of numbers will be complex on the real axis measurement.

If you do some study downstream you might wonder the effect of that at first.

[Firan Lucian]

I expected that MATLAB would have some floating point number generators (single or double).

I wonder why there is no build-in rand function for floats, like randf.

With some parameters like:

• fullRange, default
• no NAN (positive & negative option)
• no quiet NAN (positive & negative option)
• no signaling NaN (positive & negative option)
• no Inf (positive & negative option)
• no subnormals (positive & negative option)
• no normals (positive & negative option)

Similar to is* functions:

already: isfinite, isnan, isinf, ismissing

but not: isnormal, issubnormal, isnanq, isnans,...

Similar the NAN's may have sign attached to them, and 2 flavors (signaling, quiet)

https://www.doc.ic.ac.uk/~eedwards/compsys/float/nan.html

Both types of NaNs are represented by the largest biased exponent allowed by the format (single- or double-precision) and a mantissa that is non-zero.

The bit pattern of the mantissa for a signalling NaN has the most significant digit set to zero and at least one of the remaining digits set to one.

The bit pattern of the mantissa for a quiet NaN has the most significant digit set to one

A signalling NaN (NANS) is represented by any bit pattern

between 7F800001 and 7FBFFFFF or between FF800001 and FFBFFFFF

A quiet NaN (NANQ) is represented by any bit pattern

between 7FC00000 and 7FFFFFFF or between FFC00000 and FFFFFFFF

[Stephen23]

@Firan Lucian: you can:

1. write one and upload it to FEX.
2. make a feature request: https://www.mathworks.com/matlabcentral/answers/460998-how-can-i-make-a-features-request

[Bruno Luong]

It is still not clear to me the pertinent of such thing.

I was teached to manipulate random variable with some "meaningful" probability distribution, uniform on finite interval, gaussian noise, poisson, gamma, ... you name it. Something that is at least integrable and have few low-order moment well defined.

Generate such distribution related to the binary coding mantissa/exponent,sign of 4-bytes floating number?

You 'll get likely generation of many huge numbers that does appear anywhere in practice life or computer reality. That is my prediction.

[Firan Lucian]

A random float generator may be used to test/compare some 2 implementation of the same algorithm:

• on same platform but different compile options (flush to zero, denormals are zeros, different truncation, IEEE compliance …)
• on same platform using different math libraries
• on same platform using different IDE’s ( MATLAB vs gcc C vs python ….)
• on different platforms (different OS, different hardware, …),
• on different chips with specific ALU functionality ….

Either smoke test, or corner case test some specific functionality like NAN compare functionality, subnormal …

[Guillaume]

I would echo bruno's question. What's the purpose of this? Since matlab is a high-level language the stated reasons don't seem very pertinent. Matlab makes no statement about types of NaN, so even you generate a qnan or an snan, there's no guarantee that it'll be preserved. Similarly, with +/- zero. While negative zeros sometimes appear, again the documentation has nothing to say about them, so it's really an implementation detail, subject to change. Same with rounding modes, matlab does not give you control over the fp rounding modes.

Note: "A signalling NaN (NANS) is represented by [...] A quiet NaN (NANQ) is represented by [...]" vs "on different chips with specific ALU functionality"

According to wikipedia, "Encodings of qNaN and sNaN are not completely specified in IEEE 754 and depend on the processor"

[Walter Roberson]

MATLAB preserves the sign of 0:

>> a=-0;1/(a*5)

ans =

  -Inf

There just are not many cases where it matters.

MATLAB only preserves exact nans for direct assignment. Many of the routines that have to be nan aware end up generating new nans instead of copying the existing ones. For example, sort() of a vector involving nans will generate new nans instead of carefully keeping track of the kind of nan that was input.

The support for rounding modes is undocumented: https://stackoverflow.com/questions/55682034/how-to-change-the-rounding-mode-for-floating-point-operations-in-matlab

[Guillaume]

Well, yes I did say that negative 0 appear sometimes. But it's undocumented, so maybe it won't always be the case (although I doubt it'll change) and there's no guarantee that some functions won't preserve it, as is the case with nan as you've pointed out.

And that's my point, the whole thing is undocumented on purpose. As a user of a high-level language you shouldn't care about these details. If you do, maybe matlab is not the language you should be using.

[Walter Roberson]

https://www.advanpix.com/2016/04/28/branch-cuts-and-signed-zeros-in-matlab/ talks about some difficulties with MATLAB negative 0s.

Answer 1

[Stephen23]

This generates all possible singles:

>> V = ['0':'9','A':'F'];
>> X = randi(16,1,8);
>> V(X)
ans = 7B3A5B5F
>> N = typecast(uint32(hex2dec(V(X))),'single')
N = 9.6762e+035

and back again:

>> num2hex(N)
ans = 7b3a5b5f

"How about generate without NAN, or INF, or Denormalized numbers (subnormals) ?"

Read the specification for single and make sure that you do not generate those random numbers (you will probably need several randi calls, each with a different range, or to use setdiff or something similar).

See also:

https://www.mathworks.com/matlabcentral/answers/364485-convert-32-bit-hex-to-signed-floating-point

[Firan Lucian]

Send range blocks with gaps to randi might be an issue:

full SP range [0x00000000 0xffffffff] dec [0 4294967295]

+zero hex 0x00000000 dec 0

+subnormals range hex [0x00000001 0x007fffff] dec [1 8388607]

+ normals

+inf hex 0x7f800000 dec 2139095040

+NAN range [0x7f800001 0x7fffffff] dec [2139095041 2147483647]

-zero hex 0x80000000 dec 2147483648

-subnormals range [0x80000001 0x807fffff] dec [2147483649 2155872255]

- normals

-inf hex 0xff800000 dec 4286578688

-NAN range [0xff800001 0xffffffff] dec [4286578689 4294967295]

Removing NAN would mean to remove 2 ranges from full SP range, like: [0 4294967295] - [2139095041 2147483647] - [4286578689 4294967295] any thoughts on that ?

[Stephen23]

"any thoughts on that ?"

One straighforward solution would be to generate hex strings from the entire number range, and then use set or logical comparisons to detect the unwanted hex digits/combinations, and then remove those "numbers". As there are only a few special cases this would not be so difficult.

[Walter Roberson]

hex2num() to do the conversion from char to numeric.

[Stephen23]

"hex2num() to do the conversion from char to numeric."

@Walter Roberson: can you please give an example of using hex2num, which (currently) only generates double class numerics, to generate a single class numeric (as the question requests)?

I based my answer on the code given here:

https://www.mathworks.com/matlabcentral/answers/92132-how-can-i-create-a-single-precision-number-given-an-ascii-hex-string-in-matlab-7-8-r2009a#answer_101483

where MathWorks Support Team wrote: "Since the HEX2NUM only supports DOUBLE precision data, the following is a workaround for using only built-in functions to allow HEX2NUM type functionality to extent to SINGLE data". As far as I can tell, there has been no change to hex2num since that answer was written.

[Walter Roberson]

Ah good point.

Answer 2

[James Tursa]

For the generic answer with all bit patterns possible and selected with equal probability (including inf & nan & denorm) your "b" method is direct and straightforward and is the one I would use. Note that while there is only one bit pattern for -inf and one bit pattern for +inf, there are many bit patterns for nan values (any bit pattern with all 1-bits exponent and any non-zero bits mantissa will be a nan value). So it is very much more likely to get a nan with this method that it is to get a -inf or +inf. Also note that this method will generate +0 and -0 as two distinct bit patterns.

For the answer not including the special bit patterns, this can get tricky. Can we assume that you want all bit patterns except the special ones selected uniformly? If so, one could oversample and then downsize using isfinite( ) function and isdenorm( ) function. E.g., a brute force isdenorm( ) function could be something like this for single precison:

% Function returns true (element-wise) if element is denormalized number.
% s must be single precision float
% 2^(-126) is the smallest normalized single float number, realmin('single')
function g = isdenorm(s)
g = ( abs(s) < single(2^(-126)) ) & ( s ~= 0 );
end

Since the special bit patterns have exponent bits all 1's (inf and nan) or all 0's (denorm), you would simply have to oversample by about 1% or so on average (two of the 2^8 number of possible exponent bit patterns are mostly unwanted).

CAUTION: The above just shows the algorithm. A well written function would also include input argument checks which I haven't done.

*** EDIT ***

Here is some generation code based on Walter's suggestion:

% Generates random numbers for entire range of normalized single floats
% (Based on an idea by Walter Roberson)
% Note: Numbers are distributed uniformly across all possible bit patterns,
%       they are NOT distributed uniformly across the range since floating
%       point values are not uniformly distributed accross their range.
% dims = A vector containing the dimensions of the result
% Programmer: James Tursa
function r = rand_single_range(dims)
neg_sign_bit = typecast(-single(0),'uint32'); % Only the sign bit is set
smin = realmin('single'); % Smallest normalized single float
smax = realmax('single'); % Largest normalized single float
imin = typecast(smin,'uint32'); % uint32 containing bit pattern of smin
imax = typecast(smax,'uint32'); % uint32 containing bit pattern of smax
ishift = imax - imin + 1; % The amount to shift the neg values for sampling
kmax = imax + ishift + 1; % Encompass the pos + neg ranges + an extra 1 for zero
k = randi([imin kmax],prod(dims),1,'uint32'); % Generate the bit patterns
k(k==kmax) = 0; % Map the 0's first
kneg = k > imax; % Logical indexes of the (eventual) negative values
k(kneg) = k(kneg) + (neg_sign_bit - ishift); % Shift the negative values back in range and apply sign bit
r = reshape(typecast(k,'single'),dims); % Reshape result to desired dimensions
end

It generates about twice the range of bit patterns needed and then maps about half of them back into "negative" bit patterns. It does allow for a 0 bit pattern although it would be quite rare to actually get it in practice.

[Walter Roberson]

It works! ;-)

[Firan Lucian]

Nice idea.

(…The brute force approach - try generate and filter outputs - may not be practical as may require random time to complete.)

Another approach might be to have a domains <-> weights dictionary;

associate to each domain some sort of realtive probability.

Knowing the output vector (matrix) element size,

the function can compute how many elements in each domain should be generated.

Than generate for each domain.

(additionally reshuffle all generated data if needed)

domain =[‘zero+’, ‘zero-’, ‘subnormals+’, ‘subnormals-’, ‘normals+’, ‘normals-’, ‘Inf+’, ‘Inf-’, ‘NANq+’ , ‘NANq-’ , ‘NANs+’ , ‘NANs-’];

and weights=[0.0001, 0.0001, 0, 0, 0.49, 0.49, 0, 0, 0, 0, 0, 0 ];

(sum of all weights should be 1)

[Bruno Luong]

That comes down more or less to specify the probability distribution you want to generate.

[Firan Lucian]

Assemble a list of all weights:

Range list

full SP range

count 4294967296 weight 1.000000e+00 range hex [ 0x00000000 .. 0xffffffff ] dec [ 0 .. 4294967295 ] float [ 0.000000e+00 .. NaN ]

1 zero +

count 1 weight 2.328306e-10 range hex [ 0x00000000 .. 0x00000000 ] dec [ 0 .. 0 ] float [ 0.000000e+00 .. 0.000000e+00 ]

2 subnormals +

count 8388607 weight 1.953125e-03 range hex [ 0x00000001 .. 0x007fffff ] dec [ 1 .. 8388607 ] float [ 1.401298e-45 .. 1.175494e-38 ]

3 normals +

count 2130706432 weight 4.960938e-01 range hex [ 0x00800000 .. 0x7f7fffff ] dec [ 8388608 .. 2139095039 ] float [ 1.175494e-38 .. 3.402823e+38 ]

4 infinity +

count 1 weight 2.328306e-10 range hex [ 0x7f800000 .. 0x7f800000 ] dec [ 2139095040 .. 2139095040 ] float [ Inf .. Inf ]

5 NAN +

count 8388607 weight 1.953125e-03 range hex [ 0x7f800001 .. 0x7fffffff ] dec [ 2139095041 .. 2147483647 ] float [ NaN .. NaN ]

6 zero -

count 1 weight 2.328306e-10 range hex [ 0x80000000 .. 0x80000000 ] dec [ 2147483648 .. 2147483648 ] float [ -0.000000e+00 .. -0.000000e+00 ]

7 subnormals -

count 8388607 weight 1.953125e-03 range hex [ 0x80000001 .. 0x807fffff ] dec [ 2147483649 .. 2155872255 ] float [ -1.401298e-45 .. -1.175494e-38 ]

8 normals -

count 2130706432 weight 4.960938e-01 range hex [ 0x80800000 .. 0xff7fffff ] dec [ 2155872256 .. 4286578687 ] float [ -1.175494e-38 .. -3.402823e+38 ]

9 infinity -

count 1 weight 2.328306e-10 range hex [ 0xff800000 .. 0xff800000 ] dec [ 4286578688 .. 4286578688 ] float [ -Inf .. -Inf ]

10 NAN -

count 8388607 weight 1.953125e-03 range hex [ 0xff800001 .. 0xffffffff ] dec [ 4286578689 .. 4294967295 ] float [ NaN .. NaN ]

Weights list 1 2.3283e-10 0.0019531 0.49609 2.3283e-10 0.0019531

2.3283e-10 0.0019531 0.49609 2.3283e-10 0.0019531

Size list 4294967296 1 8388607 2130706432 1 8388607

1 8388607 2130706432 1 8388607

Size list 4.294967e+09 1.000000e+00 8.388607e+06 2.130706e+09 1.000000e+00 8.388607e+06 1.000000e+00 8.388607e+06 2.130706e+09 1.000000e+00 8.388607e+06

All zeros cnt 2, weight 4.6566e-10

All subnormals denormalized cnt 16777214, weight 0.0039062

All normals cnt 4261412864, weight 0.99219

All infs cnt 2, weight 4.6566e-10

All NANs cnt 16777214, weight 0.0039062

NAN subrange

5.1 signalling NAN +

count 4194303 weight 9.765623e-04 range hex [ 0x7f800001 .. 0x7fbfffff ] dec [ 2139095041 .. 2143289343 ] float [ NaN .. NaN ]

10.1 quiet NAN +

count 4194304 weight 9.765625e-04 range hex [ 0x7fc00000 .. 0x7FFFFFFF ] dec [ 2143289344 .. 2147483647 ] float [ NaN .. NaN ]

5.2 signalling NAN -

count 4194303 weight 9.765623e-04 range hex [ 0xff800001 .. 0xffbfffff ] dec [ 4286578689 .. 4290772991 ] float [ NaN .. NaN ]

10.2 quiet NAN -

count 4194304 weight 9.765625e-04 range hex [ 0xffc00000 .. 0xffffffff ] dec [ 4290772992 .. 4294967295 ] float [ NaN .. NaN ]

All signalling NANS cnt 8388608 (8.388608e+06), weight 0.0019531

All quiet NANQ cnt 8388606 (8.388608e+06), weight 0.0019531

And the code:

% code generates single precision 32bit float 
% special ranges - IEEE 754 - zero, infinty, NAN
% weigths, sign ranges 

close all
clear all
tic

w_v = zeros(1,10, 'double');
sz_v = zeros(1,10, 'double');

disp(' Range list ');

% 'full SP range [0x00000000 0xffffffff] '
[txt, w, sz] = parse_txt( '   full SP range', '00000000', 'ffffffff');
disp(txt);
w_v(1) = w; 
sz_v(1) = sz; 
%'+zero hex 0x00000000 dec  0''
[txt, w, sz] = parse_txt( ' 1       zero + ', '00000000', '00000000');
disp(txt);
w_v(2) = w;
sz_v(2) = sz; 
%+subnormals range hex [0x00000001 0x007fffff]
[txt, w, sz] = parse_txt( ' 2 subnormals + ', '00000001', '007fffff'); 
disp(txt);
w_v(3) = w;
sz_v(3) = sz; 
%+normals range hex [0x00800000 0x7f7fffff] 
[txt, w, sz] = parse_txt( ' 3    normals + ', '00800000', '7f7fffff');
disp(txt);
w_v(4) = w;
sz_v(4) = sz; 
%+infinity 0x7f800000
[txt, w, sz] = parse_txt( ' 4   infinity + ', '7f800000', '7f800000');
disp(txt);
w_v(5) = w;
sz_v(5) = sz; 
%+nan range hex [0x7f800001 0x7fffffff] 
[txt, w, sz] = parse_txt( ' 5        NAN + ', '7f800001', '7fffffff');
disp(txt);
w_v(6) = w;
sz_v(6) = sz; 

%'-zero hex 0x80000000 '
[txt, w, sz] = parse_txt( ' 6       zero - ', '80000000', '80000000');
disp(txt);
w_v(7) = w;
sz_v(7) = sz; 
%-subnormals range hex [0x80000001 0x807fffff]
[txt, w, sz] = parse_txt( ' 7 subnormals - ', '80000001', '807fffff'); 
disp(txt);
w_v(8) = w;
sz_v(8) = sz; 
%-normals range hex [0x80800000 0xff7fffff] 
[txt, w, sz] = parse_txt( ' 8    normals - ', '80800000', 'ff7fffff');
disp(txt);
w_v(9) = w;
sz_v(9) = sz; 
%-infinity 0xff800000
[txt, w, sz] = parse_txt( ' 9   infinity - ', 'ff800000', 'ff800000');
disp(txt);
w_v(10) = w;
sz_v(10) = sz; 
%-nan range hex [0x7f800001 0xffffffff] 
[txt, w, sz] = parse_txt( '10        NAN - ', 'ff800001', 'ffffffff');
disp(txt);
w_v(11) = w;
sz_v(11) = sz; 

disp(' ');
disp([ 'Weights list ' num2str(w_v) ]);
disp([ 'Size list ' num2str(sz_v) ]);
disp([ 'Size list  ' num2str(sz_v, '%e ') ]);
disp(' ');
disp([ 'All zeros cnt ' num2str(sz_v(2) + sz_v(7)) ', weight ' num2str(w_v(2) + w_v(7))]);
disp([ 'All subnormals denormalized cnt ' num2str(sz_v(3) + sz_v(8)) ', weight ' num2str(w_v(3) + w_v(8))]);
disp([ 'All normals cnt ' num2str(sz_v(4) + sz_v(9)) ', weight ' num2str(w_v(4) + w_v(9))]);
disp([ 'All infs cnt ' num2str(sz_v(5) + sz_v(10)) ', weight ' num2str(w_v(5) + w_v(10))]);
disp([ 'All NANs cnt ' num2str(sz_v(6) + sz_v(11)) ', weight ' num2str(w_v(6) + w_v(11))]);
disp(' ');
disp(' ');

% special nan's 
disp(' NAN subrange ');

%+nan signalling range hex [0x7f800001 0x7fbfffff] 
[txt, w, sz] = parse_txt( ' 5.1 signalling  NAN + ', '7f800001', '7fbfffff');
disp(txt);
w_qnan = w;
sz_qnan = sz;


%+nan quiet range hex [0x7fc00000 0x7FFFFFFF] 
[txt, w, sz] = parse_txt( '10.1      quiet  NAN + ', '7fc00000', '7FFFFFFF');
disp(txt);
w_snan = w;
sz_snan = sz;

%-nan signalling range hex [0xff800001 0xffbfffff] 
[txt, w, sz] = parse_txt( ' 5.2 signalling  NAN - ', 'ff800001', 'ffbfffff');
disp(txt);
w_qnan = w_qnan + w;
sz_qnan = sz_qnan + sz;

%-nan quiet range hex [0xffc00000 0xffffffff] 
[txt, w, sz] = parse_txt( '10.2      quiet  NAN - ', 'ffc00000', 'ffffffff');
disp(txt);
w_snan = w_snan + w;
sz_snan = sz_snan + sz;

disp(' ');
disp([ 'All signalling NANS cnt ' num2str(sz_snan) ' (' num2str(sz_snan,'%e') '), weight ' num2str(w_snan)]);
disp([ 'All quiet NANQ cnt      ' num2str(sz_qnan) ' (' num2str(sz_snan,'%e') '), weight ' num2str(w_qnan)]);


toc

function [txt, w, sz] = parse_txt(iname, hll, hul)

    % number of all SP elements
    full_size = uint64(intmax('uint32')) + 1;
    
    
    dl = uint32(hex2dec(hll));
    du = uint32(hex2dec(hul));     
    sz = uint64(du) - uint64(dl) + 1;
    w = double(sz) / double(full_size);
    fl = typecast(dl,'single');
    fu = typecast(du,'single');
    
    txt = sprintf(' %s \n count %10d weight %e range hex [ 0x%s .. 0x%s ] dec [ %10d .. %10d ] float [ %e .. %e ] ', iname, sz, w, hll, hul, dl, du, fl, fu ); 

end

[Guillaume]

Not sure what the intent of that is

full_size = intmax('uint32') - intmin('uint32') + 1;

but it's the same as:

full_size = intmax('uint32')  %ie full_size = 2^32-1

Note that intmin('uint32') is 0, intmax returns a uint32, and adding 1 to that overflows uint32 so the result is capped at intmax.

Maybe

full_size = 2^32;

would be simpler (if that was the intent).

[Firan Lucian]

Thanks Guillaume, updated code.

Having clear domain makes easy to write range checkers:

function [out] = isSubnormal(val)
% check if value is subnormal or denormalized SP number
%+subnormals range hex [0x00000001 0x007fffff]
%-subnormals range hex [0x80000001 0x807fffff]
 
    int_v  = typecast(single(val), 'uint32' );
    
    pos_range = and (( int_v >= hex2dec('00000001') ) , ( int_v <= hex2dec('007fffff') ));
    neg_range = and (( int_v >= hex2dec('80000001') ) , ( int_v <= hex2dec('807fffff') ));
    
    out = or (pos_range , neg_range);
end

Answer 3

[the cyclist]

According to the documentation here,

x = rand(10,1,"single")

will do it.

Your question has sufficient nuance that there could be some discrepancies "at the edges" that you might want to investigate further. But the above is almost certainly adequate for the vast majority of purposes.

[Stephen23]

Due to the precision of rand's output, scaling those values will miss some numbers, e.g.:

>> one = single(1);
>> rmx = realmax('single');
>> A = rmx*(one) % largest
A =
   3.4028235e+38
>> B = rmx*(one-eps(one)) % scaled second largest
B =
   3.4028231e+38
>> C = rmx-eps(rmx) % missed second largest
C =
   3.4028233e+38
>> num2hex(A) % largest
ans =
7f7fffff
>> num2hex(B) % scaled second largest
ans =
7f7ffffd
>> num2hex(C) % missed second largest
ans =
7f7ffffe

[Firan Lucian]

Added c variant rand(size_x, size_y, "single") on large size elements, but range seems smaller. Best variant yet seems b.

size_x = 100000000;
size_y = 1;

a = realmin('single') + ( realmax('single') - realmin('single') ).*rand(size_x, size_y, 'single');
disp( [ sprintf("Var a Min: %e, Max: %e,  Mean: %e, Std: %e", min(a), max(a), mean(a), std(a) ) ] );
t1a = kstest(a);
t2a = adtest(a);
t3a = jbtest(a);
disp( [ sprintf("      ktest: %d, adtest: %d, jbtest: %d ", t1a, t2a, t3a ) ] );
% remove nan and inf
a = a( ~isnan( a ) );
a = a( ~isinf( a ) );
a = double(a);
disp( [ sprintf("      Mean: %e, Std: %e ", mean( a ), std( a ) ) ] ); 

b = typecast(randi([0, intmax('uint32')], size_x, size_y, 'uint32'), 'single');
disp( [ sprintf("Var b Min: %e, Max: %e,  Mean: %e, Std: %e ", min(b), max(b), mean(b), std(b) ) ] );
t1b = kstest(b);
t2b = adtest(b);
t3b = jbtest(b);
disp( [ sprintf("      ktest: %d, adtest: %d, jbtest: %d ", t1b, t2b, t3b ) ] );
% remove nan and inf
b = b( ~isnan( b ) );
b = b( ~isinf( b ) );
b = double(b);
disp( [ sprintf("      Mean: %e, Std: %e ", mean( b ), std( b ) ) ] ); 

c = rand(size_x, size_y, "single");
disp( [ sprintf("Var c Min: %e, Max: %e,  Mean: %e, Std: %e ", min(c), max(c), mean(c), std(c) ) ] );
t1c = kstest(c);
t2c = adtest(c);
t3c = jbtest(c);
disp( [ sprintf("      ktest: %d, adtest: %d, jbtest: %d ", t1c, t2c, t3c ) ] );
c = c( ~isnan( c ) );
c = c( ~isinf( c ) );
c = double(c);
disp( [ sprintf("      Mean: %e, Std: %e ", mean( c ), std( c ) ) ] ); 

Var a Min: 5.572502e+29, Max: 3.402823e+38, Mean: Inf, Std: Inf

ktest: 1, adtest: 1, jbtest: 1

Mean: 1.701168e+38, Std: 9.822943e+37

Var b Min: -3.402817e+38, Max: 3.402815e+38, Mean: NaN, Std: NaN

ktest: 1, adtest: 1, jbtest: 1

Mean: 2.128702e+32, Std: 1.878338e+37

Var c Min: 1.400272e-09, Max: 9.999999e-01, Mean: 4.999776e-01, Std: 2.886557e-01

Warning: P is less than the smallest tabulated value, returning 0.0005.

> In adtest (line 279)

In SP_rand (line 32)

ktest: 1, adtest: 1, jbtest: 1

Mean: 4.999774e-01, Std: 2.886558e-01

[Firan Lucian]

Another issue of variant a and c (simple rand) is that it does not output negative numbers.

Var a Min: 5.023105e+30, Max: 3.402823e+38, Mean: Inf, Std: Inf

rel count: is+ 100.00, is- 0.00, is0 0.00, Inf 0.000000e+00, NAN 0.000000e+00

ktest: 1, adtest: 1, jbtest: 1

Mean: 1.701390e+38, Std: 9.822651e+37

Var b Min: -3.402819e+38, Max: Inf, Mean: NaN, Std: NaN

rel count: is+ 49.80, is- 49.81, is0 0.00, Inf 0.000000e+00, NAN 0.000000e+00

ktest: 1, adtest: 1, jbtest: 1

Mean: 6.011336e+32, Std: 1.879643e+37

Var c Min: 2.171108e-08, Max: 9.999999e-01, Mean: 4.999863e-01, Std: 2.886831e-01

rel count: is+ 100.00, is- 0.00, is0 0.00, Inf 0.000000e+00, NAN 0.000000e+00

Warning: P is less than the smallest tabulated value, returning 0.0005.

> In adtest (line 279)

In SP_rand (line 52)

ktest: 1, adtest: 1, jbtest: 1

Mean: 4.999864e-01, Std: 2.886831e-01

ipositif = sum(sign(a(:))==1) /size_x * 100;
inegatif = sum(sign(a(:))==-1) /size_x * 100;
izero = sum(a(:)==0) /size_x * 100;
iinf = sum(isinf(a)) /size_x * 100;
inan = sum(isnan(a)) /size_x * 100;
disp( [ sprintf("      rel count: is+ %.02f, is- %.02f, is0 %.02f, Inf %e, NAN %e", ipositif, inegatif, izero, iinf, inan ) ] );

Answer 4

[Walter Roberson]

For restricting to finite normalized values, randi() from the decimal representation of the smallest positive normalized number, to the largest positive normalized number times 2 plus 1. Values up to the largest positive get converted directly with typecast to single. Map the next group to negative by reflection and typecast and multiply by -1. The final upper value you map to 0 exactly.

This method should not require any rejection.

Answer 5

[Firan Lucian]

Using bitwise operations:

>> rand_SP_real

init 00000000000000000000000000000000

gen exponent 00010101

bitshift exponent 00001010100000000000000000000000

Add exponent 00001010100000000000000000000000

gen mantissa 00111010100111100011110

Add mantissa 00001010100111010100111100011110

ans =

single

1.5148e-32

>> rand_SP_real

init 00000000000000000000000000000000

Add negative sign 10000000000000000000000000000000

gen exponent 00100111

bitshift exponent 00010011100000000000000000000000

Add exponent 10010011100000000000000000000000

gen mantissa 11010011011010001011110

Add mantissa 10010011111010011011010001011110

ans =

single

-5.8995e-27

>> rand_SP_real

init 00000000000000000000000000000000

Add negative sign 10000000000000000000000000000000

gen exponent 11111111

bitshift exponent 01111111100000000000000000000000

Add exponent 11111111100000000000000000000000

gen mantissa 00010100000000110100111

Add mantissa 11111111100010100000000110100111

ans =

single

NaN

function [out] = rand_SP_real()
% generate a random SP number 
% does generate zero, subnormals, nan, inf

    %full range 
    % sp_sign  = randi(2^1)-1;    
    % sp_exp   = randi(2^8)-1;    
    % sp_mant  = randi(2^23)-1;
    
    sp = uint32(0);
    disp(['init              ' dec2bin(sp, 32) ]);
    
    if randi(2^1) == 2
        %negative sign
        sp = uint32(bitset(sp, 32));  
        disp(['Add negative sign ' dec2bin(sp, 32)]);
    end
    
    sp_exp   = randi(2^8)-1; 
           
    disp(['gen exponent       ' dec2bin(sp_exp, 8)]);
    % bitshift exponent     
    sp_exp  = uint32(bitshift(sp_exp, 23, 'uint32'));
    disp(['bitshift exponent ' dec2bin(sp_exp, 32)]);
    % add exponent
    sp = uint32(sp + sp_exp) ;
    disp(['Add exponent      ' dec2bin(sp, 32)]);
    % rand mantissa
    sp_mant  = randi(2^23)-1;
    disp(['gen mantissa               ' dec2bin(sp_mant, 23)]);
    sp = uint32(sp + sp_mant) ;
    disp(['Add mantissa      ' dec2bin(sp, 32) ]);
    
    
    % convert to single
    out = typecast( uint32(sp), 'single');    
end

>> rand_SP_real_0normals

init 00000000000000000000000000000000

Add negative sign 10000000000000000000000000000000

gen exponent 10110010

bitshift exponent 01011001000000000000000000000000

Add exponent 11011001000000000000000000000000

gen mantissa 11100100000100100011110

Add mantissa 11011001011100100000100100011110

ans =

single

-4.2579e+15

function [out] = rand_SP_real_0normals()
% generate a random SP normal number and zero
% does not generate subnormals, nan, inf

    %full range 
    % sp_sign  = randi(2^1)-1;    
    % sp_exp   = randi(2^8)-1;    
    % sp_mant  = randi(2^23)-1;

    
    sp = uint32(0);
    disp(['init              ' dec2bin(sp, 32) ]);
    
    if randi(2^1) == 2
        %negative sign
        sp = uint32(bitset(sp, 32));  
        disp(['Add negative sign ' dec2bin(sp, 32)]);
    end
    
    % remove inf nan from exponent (-1)
    sp_exp   = randi(2^8-1)-1; 
    
    % remove subnormals (denormalized), if exponent == 0
    if sp_exp ~= 0             
        disp(['gen exponent       ' dec2bin(sp_exp, 8)]);
        % bitshift exponent     
        sp_exp  = uint32(bitshift(sp_exp, 23, 'uint32'));
        disp(['bitshift exponent ' dec2bin(sp_exp, 32)]);
        % add exponent
        sp = uint32(sp + sp_exp) ;
        disp(['Add exponent      ' dec2bin(sp, 32)]);
        % rand mantissa
        sp_mant  = randi(2^23)-1;
        disp(['gen mantissa               ' dec2bin(sp_mant, 23)]);
        sp = uint32(sp + sp_mant) ;
        disp(['Add mantissa      ' dec2bin(sp, 32) ]);
    end
    
    % convert to single
    out = typecast( uint32(sp), 'single');    
 end 

>> rand_SP_real_infnan

init 00000000000000000000000000000000

gen exponent 11111111

bitshift exponent 01111111100000000000000000000000

Add exponent 01111111100000000000000000000000

gen mantissa 00011011010011011001100

Add mantissa 01111111100011011010011011001100

ans =

single

NaN

>> rand_SP_real_infnan

init 00000000000000000000000000000000

Add negative sign 10000000000000000000000000000000

gen exponent 11111111

bitshift exponent 01111111100000000000000000000000

Add exponent 11111111100000000000000000000000

gen mantissa 00000001001011111011010

Add mantissa 11111111100000001001011111011010

ans =

single

NaN

function [out] = rand_SP_real_infnan()
% generate a random SP number in nan range and inf
% does generate subnormals, normals, zero 

    %full range 
    % sp_sign  = randi(2^1)-1;    
    % sp_exp   = randi(2^8)-1;    
    % sp_mant  = randi(2^23)-1;
    
    sp = uint32(0);
    disp(['init              ' dec2bin(sp, 32) ]);
    
    if randi(2^1) == 2
        %negative sign
        sp = uint32(bitset(sp, 32));  
        disp(['Add negative sign ' dec2bin(sp, 32)]);
    end
    
    % only one exponent
    sp_exp   = 255; 
           
    disp(['gen exponent       ' dec2bin(sp_exp, 8)]);
    % bitshift exponent     
    sp_exp  = uint32(bitshift(sp_exp, 23, 'uint32'));
    disp(['bitshift exponent ' dec2bin(sp_exp, 32)]);
    % add exponent
    sp = uint32(sp + sp_exp) ;
    disp(['Add exponent      ' dec2bin(sp, 32)]);
    % rand mantissa
    sp_mant  = randi(2^23)-1;
    disp(['gen mantissa               ' dec2bin(sp_mant, 23)]);
    sp = uint32(sp + sp_mant) ;
    disp(['Add mantissa      ' dec2bin(sp, 32) ]);
    
    
    % convert to single
    out = typecast( uint32(sp), 'single');    
 end

>> rand_SP_real_subnormals0

init 00000000000000000000000000000000

Add negative sign 10000000000000000000000000000000

gen exponent 00000000

bitshift exponent 00000000000000000000000000000000

Add exponent 10000000000000000000000000000000

gen mantissa 11010001001110101100100

Add mantissa 10000000011010001001110101100100

ans =

single

-9.6074e-39

function [out] = rand_SP_real_subnormals0()
% generate a random SP number in subnormals range and zero 
% does generate normals, nan, inf range

    %full range 
    % sp_sign  = randi(2^1)-1;    
    % sp_exp   = randi(2^8)-1;    
    % sp_mant  = randi(2^23)-1;
    
    sp = uint32(0);
    disp(['init              ' dec2bin(sp, 32) ]);
    
    if randi(2^1) == 2
        %negative sign
        sp = uint32(bitset(sp, 32));  
        disp(['Add negative sign ' dec2bin(sp, 32)]);
    end
    
    % only one exponent
    sp_exp   = 0; 
           
    disp(['gen exponent       ' dec2bin(sp_exp, 8)]);
    % bitshift exponent     
    sp_exp  = uint32(bitshift(sp_exp, 23, 'uint32'));
    disp(['bitshift exponent ' dec2bin(sp_exp, 32)]);
    % add exponent
    sp = uint32(sp + sp_exp) ;
    disp(['Add exponent      ' dec2bin(sp, 32)]);
    % rand mantissa
    sp_mant  = randi(2^23)-1;
    disp(['gen mantissa               ' dec2bin(sp_mant, 23)]);
    sp = uint32(sp + sp_mant) ;
    disp(['Add mantissa      ' dec2bin(sp, 32) ]);
    
    
    % convert to single
    out = typecast( uint32(sp), 'single');    
end

[Firan Lucian]

function [out] = isSubnormal(val)
% check if value is subnormal or denormalized SP number, zero excluded
%+subnormals range hex [0x00000001 0x007fffff]
%-subnormals range hex [0x80000001 0x807fffff]

    int_v  = typecast(single(val), 'uint32' );
    
    pos_range = and (( int_v >= hex2dec('00000001') ) , ( int_v <= hex2dec('007fffff') ));
    neg_range = and (( int_v >= hex2dec('80000001') ) , ( int_v <= hex2dec('807fffff') ));
    
    out = or (pos_range , neg_range);
end

[Steven Lord]

There's no need to convert to uint32 and call hex2dec. Just use realmin directly. As a bonus, this will work for either double or single inputs.

function tf = isSubnormal(value)
smallestPositiveNormalized = realmin(class(value));
tf = abs(value) < smallestPositiveNormalized;
end

[James Tursa]

You need to check for 0 also, since 0 is not a denorm.

[Steven Lord]

As soon as I saw there was a comment on this question I realized the point James made. Yes, you need to exclude 0.