Just in case anyone is interested in using the final version of the PRNG_mod code that is discussed above in this thread, here it is. This includes both the division and geometric exponent distribution codes for uniform prngs in [0.0,1.0) and the normal distribution codes (mean=0.0 and standard devision=1.0). The normal distribution codes do not have the tail truncation problem. I also tested the code for default int64 integers, and it appears to work as expected. The documentation is included in the file, along with a main program to test the features. See the documentation for further details.
prng_mod.F90 (79.2 KB)
Here is the gfortran -O3 output for the code.
gfortran.txt (13.8 KB)
edit: I wanted to add a comment about the output file. The uniform distribution numbers in that file do not seem to depend on the optimization level or the compiler options. I did not experiment with every possible combination, but when I did make changes, those histograms don’t seem to change. The reason is that those prng codes actually do little or no actual floating point arithmetic. The bit strings that are generated with the xoshiro256** algorithm are all done with integer arithmetic operations, and those results do not depend on optimization level. Then with the GD algorithms, those bit strings are taken and placed into the fraction field of an integer that will eventually be the floating point number. The exponent is determined by computing leading zero counts from the bit strings, so that is not a floating point operation either. Then the biased exponent is computed and placed into the integer. Then that integer is converted to the floating point value with a transfer() operation, which is also not a true floating point operation since there is no chance of rounding occurring and there are no overflows, underflows, subnonrmals, or NaNs that are generated. So those results are all really just integer operations that do not depend on optimization levels or floating point compiler options. So far, they even appear to be the same when changing compilers and underlying hardware (provided the same IEEE floating point formats are used, the same endian addressing applies, etc.).
The division algorithm does have one floating point operation in the expression real(ix,real32). I think in principle that result could be affected by floating point options, but in practice on the compilers I’ve tested it does not seem to be. At that point the ix integer value can be represented exactly as a floating point value, so there is no rounding that is required by that conversion, and if rounding does occur at that point, one might even argue, from a quality of implementation perspective, that the compiler is making a mistake. The next step in the division algorithm is to divide that result by 2.0**24. However, instead of doing that division, or instead of multiplying by its reciprocal 2.0**(-24), where floating point rounding might occur, the code instead uses the intrinsic scale() function. I would argue that that is not an actual floating point operation either. The scale operation works by extracting the exponent field from the floating point value, adding the integer argument, testing for integer exponent underflow or underflow (which cannnot occur in this code), and then moving those new exponent bits back into the floating point exponent field. So even though the result is a floating point value, there was never any actual floating point operations performed by the floating point functional units, it is instead just some integer operations.
However, the normal distribution histograms do involve some floating point operations. The code uses the Box-Muller transform to change two uniform distribution values into two normal distribution values. This code actually does a vector version of that transform; i.e. it changes two vectors of uniform values into two vectors of normal values and returns the full set of results in a cache buffer. That transformation requires a log(), sqrt(), cos(), sin() and three multiplications, for each pair of output values. These operations all use the floating point units where rounding can affect the result. Thus compiler options that affect the rounding mode, the use of extended precision registers, the use of AVX vector hardware, or the use of GPU hardware can all affect the computed normal distribution values. When you compare the histogram values to those in the reference output, you might expect to see some small shifts where values are assigned to one histogram instead of its neighbor. I see such shifts in the real32 values when I change optimization levels with gfortran.