# HG changeset patch
# User Philipp Kutin <philipp.kutin@gmail.com>
# Date 1376590470 -7200
# Node ID 213ee68b59dac9c8a6a79b25798e66512a102eaa
# Parent  ce55a85758742e9fef1d0758273f66e04bb56f2e
Handle out-of-range values consistently when initializing random number generator.

* oct-rand.cc: New function double2uint32().  It normalizes a finite double by
  taking it modulo 2^32, such that the result is in [0, 2^32).  Nonfinites give
  zero as result.  This prevents compiler-dependent casting of a double whose
  truncated value cannot be represented in a uint32_t.

* rand.cc: Add %!test block to verify behavior.

diff -r ce55a8575874 -r 213ee68b59da libinterp/corefcn/rand.cc
--- a/libinterp/corefcn/rand.cc	Thu Aug 15 14:17:20 2013 -0700
+++ b/libinterp/corefcn/rand.cc	Thu Aug 15 20:14:30 2013 +0200
@@ -533,6 +533,19 @@
 %! endif
 */
 
+/*
+%!# Test out-of-range values as rand() seeds.  See oct-rand.cc: double2uint32().
+%!function v = __rand_sample__ (initval)
+%!  rand ("state", initval);
+%!  v = rand (1, 6);
+%!endfunction
+%!
+%!assert (__rand_sample__ (0), __rand_sample__ (2^32))
+%!assert (__rand_sample__ (-2), __rand_sample__ (2^32-2))
+%!assert (__rand_sample__ (Inf), __rand_sample__ (NaN))
+%!assert (! isequal (__rand_sample__ (-1), __rand_sample__ (-2)))
+*/
+
 static std::string current_distribution = octave_rand::distribution ();
 
 DEFUN (randn, args, ,
diff -r ce55a8575874 -r 213ee68b59da liboctave/numeric/oct-rand.cc
--- a/liboctave/numeric/oct-rand.cc	Thu Aug 15 14:17:20 2013 -0700
+++ b/liboctave/numeric/oct-rand.cc	Thu Aug 15 20:14:30 2013 +0200
@@ -663,6 +663,30 @@
   return retval;
 }
 
+// Guarantee reproducible conversion of negative initialization values to
+// random number algorithm.  Note that Matlab employs slightly different rules.
+// 1) Seed saturates at 2^32-1 for any value larger than that.
+// 2) NaN, Inf are translated to 2^32-1.
+// 3) -Inf is translated to 0.
+static uint32_t
+double2uint32 (double d)
+{
+  uint32_t u;
+  static const double TWOUP32 = std::numeric_limits<uint32_t>::max() + 1.0;
+
+  if (! xfinite (d))
+    u = 0;
+  else
+    {
+      d = fmod (d, TWOUP32);
+      if (d < 0)
+        d += TWOUP32;
+      u = static_cast<uint32_t> (d);
+    }
+
+  return u;
+}
+
 void
 octave_rand::set_internal_state (const ColumnVector& s)
 {
@@ -672,7 +696,7 @@
   OCTAVE_LOCAL_BUFFER (uint32_t, tmp, MT_N + 1);
 
   for (octave_idx_type i = 0; i < n; i++)
-    tmp[i] = static_cast<uint32_t> (s.elem (i));
+    tmp[i] = double2uint32 (s.elem (i));
 
   if (len == MT_N + 1 && tmp[MT_N] <= MT_N && tmp[MT_N] > 0)
     oct_set_state (tmp);