diff --git a/applications/test/Random/Test-Random.C b/applications/test/Random/Test-Random.C
index a724fe8357..a0635488b0 100644
--- a/applications/test/Random/Test-Random.C
+++ b/applications/test/Random/Test-Random.C
@@ -68,6 +68,22 @@ double randomFraction(const uint64_t bits)
 
 using namespace Foam;
 
+// Test uniformity of random
+void testPosition(const label n)
+{
+    List<label> samples(n, Zero);
+
+    Random rnd(123456);
+    for (label i=0; i < 100000*n; ++i)
+    {
+        ++samples[rnd.position<label>(0,n-1)];
+    }
+
+    Info<< nl << "uniform [0," << n << ")\n  "
+        << flatOutput(samples) << nl;
+}
+
+
 // Output with cout instead of Info to keep values unsigned on output
 using std::cout;
 using std::setw;
@@ -175,20 +191,24 @@ int main(int argc, char *argv[])
     }
 
     // Test uniformity of random
+    testPosition(20);
+    testPosition(3);
+
+    // This should fail (in FULLDEBUG)
+    const bool throwingError = FatalError.throwExceptions();
+    try
     {
-        List<label> samples(20, Zero);
-
-        Random rnd(123456);
-        for (label i=0; i < 1000*samples.size(); ++i)
-        {
-            ++samples[rnd.position<label>(0,19)];
-        }
-
-        Info<< nl << "uniform [0,20)" << nl << "  "
-            << flatOutput(samples) << nl;
+        Info<<"Random position(10,5): "
+            << Random().position<label>(10, 5) << endl;
+    }
+    catch (Foam::error& err)
+    {
+        Info<< "Caught FatalError " << err << nl << endl;
     }
 
-    Info<< nl << "Done." << endl;
+    FatalError.throwExceptions(throwingError);
+
+    Info<< "\nDone" << nl << endl;
 
     return 0;
 }
diff --git a/src/OpenFOAM/primitives/random/Random/Random.C b/src/OpenFOAM/primitives/random/Random/Random.C
index ed09267b12..338bd83bd6 100644
--- a/src/OpenFOAM/primitives/random/Random/Random.C
+++ b/src/OpenFOAM/primitives/random/Random/Random.C
@@ -121,9 +121,24 @@ Foam::scalar Foam::Random::position
 template<>
 Foam::label Foam::Random::position(const label& start, const label& end)
 {
-    // Extend range from [0, N-1] to (-0.5, N-0.5) to ensure that round()
-    // results in the same number density at the ends.
-    return start + round(scalar01()*((end - start) + 0.998) - 0.499);
+    #ifdef FULLDEBUG
+    if (start > end)
+    {
+        FatalErrorInFunction
+            << "start index " << start << " > end index " << end << nl
+            << abort(FatalError);
+    }
+    #endif
+
+    // Extend the upper sampling range by 1 and floor the result.
+    // Since the range is non-negative, can use integer truncation
+    // instead using floor().
+
+    const label val = start + label(scalar01()*(end - start + 1));
+
+    // Rare case when scalar01() returns exactly 1.000 and the truncated
+    // value would be out of range.
+    return min(val, end);
 }
 
 
diff --git a/src/OpenFOAM/primitives/random/Random/Random.H b/src/OpenFOAM/primitives/random/Random/Random.H
index 58d489b3e7..d09424c668 100644
--- a/src/OpenFOAM/primitives/random/Random/Random.H
+++ b/src/OpenFOAM/primitives/random/Random/Random.H
@@ -61,10 +61,10 @@ class Random
         //- Initial random number seed
         label seed_;
 
-        //- Random number generator on the int32 interval [0, 2^31)
+        //- Random number generator on the int32 interval [0,2^31)
         Rand48 generator_;
 
-        //- Uniform distribution on the scalar interval [0, 1]
+        //- Uniform distribution on the scalar interval [0,1]
         std::uniform_real_distribution<scalar> uniform01_;
 
         //- Is there a gaussian sample cached?
@@ -91,7 +91,7 @@ public:
         Random(const Random& r, const bool reset = false);
 
 
-    // Destructor
+    //- Destructor
     ~Random() = default;
 
 
@@ -111,20 +111,20 @@ public:
         //- Return a random bit
         inline int bit();
 
-        //- Return a sample whose components lie in the range 0-1
+        //- Return a sample whose components lie in the range [0,1]
         template<class Type>
         Type sample01();
 
         //- Return a sample whose components are normally distributed
-        //  with zero mean and unity variance N(0, 1)
+        //- with zero mean and unity variance N(0,1)
         template<class Type>
         Type GaussNormal();
 
-        //- Return a sample between start and end
+        //- Return a sample on the interval [start,end]
         template<class Type>
         Type position(const Type& start, const Type& end);
 
-        //- Randomise value in the range 0-1
+        //- Randomise value in the range [0,1]
         template<class Type>
         void randomise01(Type& value);
 
@@ -135,16 +135,16 @@ public:
 
     // Global random numbers - consistent across all processors
 
-        //- Return a sample whose components lie in the range 0-1
+        //- Return a sample whose components lie in the range [0,1]
         template<class Type>
         Type globalSample01();
 
         //- Return a sample whose components are normally distributed
-        //  with zero mean and unity variance N(0, 1)
+        //- with zero mean and unity variance N(0,1)
         template<class Type>
         Type globalGaussNormal();
 
-        //- Return a sample between start and end
+        //- Return a sample on the interval [start,end]
         template<class Type>
         Type globalPosition(const Type& start, const Type& end);