Implement uniform_real/0 and uniform_real_s/1

1 files changed, 113 insertions, 5 deletions

diff --git a/lib/stdlib/doc/src/rand.xml b/lib/stdlib/doc/src/rand.xml
index 89fb858823..21f680a0ee 100644
--- a/lib/stdlib/doc/src/rand.xml
+++ b/lib/stdlib/doc/src/rand.xml
@@ -133,8 +133,9 @@
       variable <c>rand_seed</c> to remember the current state.</p>
 
     <p>If a process calls
-      <seealso marker="#uniform-0"><c>uniform/0</c></seealso> or
-      <seealso marker="#uniform-1"><c>uniform/1</c></seealso> without
+      <seealso marker="#uniform-0"><c>uniform/0</c></seealso>,
+      <seealso marker="#uniform-1"><c>uniform/1</c></seealso> or
+      <seealso marker="#uniform_real-0"><c>uniform_real/0</c></seealso> without
       setting a seed first, <seealso marker="#seed-1"><c>seed/1</c></seealso>
       is called automatically with the default algorithm and creates a
       non-constant seed.</p>
@@ -168,10 +169,17 @@ R3 = rand:uniform(),</pre>
 S0 = rand:seed_s(exrop),
 {R4, S1} = rand:uniform_s(S0),</pre>
 
+   <p>Textbook basic form Box-Muller standard normal deviate</p>
+
+   <pre>
+R5 = rand:uniform_real(),
+R6 = rand:uniform(),
+SND0 = math:sqrt(-2 * math:log(R5)) * math:cos(math:pi() * R6)</pre>
+
    <p>Create a standard normal deviate:</p>
 
    <pre>
-{SND0, S2} = rand:normal_s(S1),</pre>
+{SND1, S2} = rand:normal_s(S1),</pre>
 
    <p>Create a normal deviate with mean -3 and variance 0.5:</p>
 
@@ -414,7 +422,8 @@ tests. We suggest to use a sign test to extract a random Boolean value.</pre>
 	    This function may return exactly <c>0.0</c> which can be
 	    fatal for certain applications.  If that is undesired
 	    you can use <c>(1.0 - rand:uniform())</c> to get the
-	    interval <c>0.0 &lt; <anno>X</anno> =&lt; 1.0</c>.
+	    interval <c>0.0 &lt; <anno>X</anno> =&lt; 1.0</c>, or instead use
+	    <seealso marker="#uniform_real-0"><c>uniform_real/0</c></seealso>.
 	  </p>
 	  <p>
 	    If neither endpoint is desired you can test and re-try
@@ -432,6 +441,42 @@ end.</pre>
     </func>
 
     <func>
+      <name name="uniform_real" arity="0"/>
+      <fsummary>Return a random float.</fsummary>
+      <desc><marker id="uniform_real-0"/>
+        <p>
+	  Returns a random float
+	  uniformly distributed in the value range
+	  <c>DBL_MIN =&lt; <anno>X</anno> &lt; 1.0</c>
+	  and updates the state in the process dictionary.
+	</p>
+	<p>
+	  Conceptually, a random real number <c>R</c> is generated
+	  from the interval <c>0 =&lt; R &lt; 1</c> and then the
+	  closest rounded down normalized number
+	  in the IEEE 754 Double precision format
+	  is returned.
+	</p>
+	<note>
+	  <p>
+	    The generated numbers from this function has got better
+	    granularity for small numbers than the regular
+	    <seealso marker="#uniform-0"><c>uniform/0</c></seealso>
+	    because all bits in the mantissa are random.
+	    This property, in combination with the fact that exactly zero
+	    is never returned is useful for algoritms doing for example
+	    <c>1.0 / <anno>X</anno></c> or <c>math:log(<anno>X</anno>)</c>.
+	  </p>
+	</note>
+	<p>
+	  See
+	  <seealso marker="#uniform_real_s-1"><c>uniform_real_s/1</c></seealso>
+	  for more explanation.
+	</p>
+      </desc>
+    </func>
+
+    <func>
       <name name="uniform" arity="1"/>
       <fsummary>Return a random integer.</fsummary>
       <desc><marker id="uniform-1"/>
@@ -460,7 +505,8 @@ end.</pre>
 	    This function may return exactly <c>0.0</c> which can be
 	    fatal for certain applications.  If that is undesired
 	    you can use <c>(1.0 - rand:uniform(State))</c> to get the
-	    interval <c>0.0 &lt; <anno>X</anno> =&lt; 1.0</c>.
+	    interval <c>0.0 &lt; <anno>X</anno> =&lt; 1.0</c>, or instead use
+	    <seealso marker="#uniform_real_s-1"><c>uniform_real_s/1</c></seealso>.
 	  </p>
 	  <p>
 	    If neither endpoint is desired you can test and re-try
@@ -478,6 +524,68 @@ end.</pre>
     </func>
 
     <func>
+      <name name="uniform_real_s" arity="1"/>
+      <fsummary>Return a random float.</fsummary>
+      <desc>
+        <p>
+	  Returns, for a specified state, a random float
+	  uniformly distributed in the value range
+	  <c>DBL_MIN =&lt; <anno>X</anno> &lt; 1.0</c>
+	  and updates the state in the process dictionary.
+	</p>
+	<p>
+	  Conceptually, a random real number <c>R</c> is generated
+	  from the interval <c>0 =&lt; R &lt; 1</c> and then the
+	  closest rounded down normalized number
+	  in the IEEE 754 Double precision format
+	  is returned.
+	</p>
+	<note>
+	  <p>
+	    The generated numbers from this function has got better
+	    granularity for small numbers than the regular
+	    <seealso marker="#uniform_s-1"><c>uniform_s/1</c></seealso>
+	    because all bits in the mantissa are random.
+	    This property, in combination with the fact that exactly zero
+	    is never returned is useful for algoritms doing for example
+	    <c>1.0 / <anno>X</anno></c> or <c>math:log(<anno>X</anno>)</c>.
+	  </p>
+	</note>
+	<p>
+	  The concept implicates that the probability to get
+	  exactly zero is extremely low; so low that this function
+	  is in fact guaranteed to never return zero.  The smallest
+	  number that it might return is <c>DBL_MIN</c>, which is
+	  2.0^(-1022).
+	</p>
+	<p>
+	  The value range stated at the top of this function
+	  description is technically correct, but
+	  <c>0.0 =&lt; <anno>X</anno> &lt; 1.0</c>
+	  is a better description of the generated numbers'
+	  statistical distribution.  Except that exactly 0.0
+	  is never returned, which is not possible to observe
+	  statistically.
+	</p>
+	<p>
+	  For example; for all sub ranges
+	  <c>N*2.0^(-53) =&lt; X &lt; (N+1)*2.0^(-53)</c>
+	  where
+	  <c>0 =&lt; integer(N) &lt; 2.0^53</c>
+	  the probability is the same.
+	  Compare that with the form of the numbers generated by
+	  <seealso marker="#uniform_s-1"><c>uniform_s/1</c></seealso>.
+	</p>
+	<p>
+	  Having to generate extra random bits for
+	  small numbers costs a little performance.
+	  This function is about 20% slower than the regular
+	  <seealso marker="#uniform_s-1"><c>uniform_s/1</c></seealso>
+	</p>
+      </desc>
+    </func>
+
+    <func>
       <name name="uniform_s" arity="2"/>
       <fsummary>Return a random integer.</fsummary>
       <desc>