diff --git a/README.md b/README.md
index c344ae61..32889858 100644
--- a/README.md
+++ b/README.md
@@ -13,12 +13,38 @@ Professional Random Number Generators
 
 Documentation: http://docs.random.dlang.io/latest/index.html
 
+
+### Example (3 seconds)
+```d
+#!/usr/bin/env dub
+/+ dub.json:
+{
+    "name": "test_random",
+    "dependencies": {"mir-random": "~>0.3.3"}
+}
++/
+
+import std.range, std.stdio;
+
+import mir.random.variable: NormalVariable;
+import mir.random.algorithm: range;
+
+void main()
+{
+    auto sample = NormalVariable!double(0, 1).range.take(10).array;
+
+    sample[$.randIndex].writeln; // prints random element from the sample
+}
+```
+
+
+### Example (10 seconds)
 ```d
 #!/usr/bin/env dub
 /+ dub.json:
 {
     "name": "test_random",
-    "dependencies": {"mir-random": "~>0.3.1"}
+    "dependencies": {"mir-random": "~>0.3.3"}
 }
 +/
 
@@ -31,12 +57,12 @@ import mir.random.algorithm: range;
 
 void main(){
     auto rng = Random(unpredictableSeed);        // Engines are allocated on stack or global
-    auto sample = range!rng                      // Engines can passed by alias to algorithms
+    auto sample = range!rng                      // Engines can passed to algorithms by alias or by pointer
         (NormalVariable!double(0, 1))            // Random variables are passed by value
-        .take(1000)                              // Fix sample length to 1000 elements (Input Range API)
+        .take(10)                                // Fix sample length to 10 elements (Input Range API)
         .array;                                  // Allocates memory and performs computation
 
-    writeln(sample);
+    sample[$.randIndex].writeln;                 // prints random element from the sample
 }
 ```
 
diff --git a/source/mir/random/algorithm.d b/source/mir/random/algorithm.d
index 4540622a..4fb6223c 100644
--- a/source/mir/random/algorithm.d
+++ b/source/mir/random/algorithm.d
@@ -11,6 +11,7 @@ import mir.math.common;
 
 import mir.random;
 public import mir.random.engine;
+import mir.random.variable: isRandomVariable;
 
 private enum bool isPassByRef(P) = !is(P == struct) && !is(P == union) && !isPointer!P && !isScalarType!P
     && (isFunction!P || isDelegate!P || is(P == class) || is(P == interface));
@@ -140,28 +141,28 @@ struct RandomField(alias gen)
 
 /// ditto
 RandomField!(G, D, T) field(T, G, D)(ref G gen, D var)
-    if (isReferenceToSaturatedRandomEngine!G)
+    if (isReferenceToSaturatedRandomEngine!G && isRandomVariable!D)
 {
     return typeof(return)(gen, var);
 }
 
 /// ditto
 RandomField!(gen, D, T) field(alias gen, D, T)(D var)
-    if (isSaturatedRandomEngine!(typeof(gen)))
+    if (isSaturatedRandomEngine!(typeof(gen)) && isRandomVariable!D)
 {
     return RandomField!(gen,D,T)(var);
 }
 
 /// ditto
 auto field(G, D)(ref G gen, D var)
-    if (isReferenceToSaturatedRandomEngine!G)
+    if (isReferenceToSaturatedRandomEngine!G && isRandomVariable!D)
 {
     return RandomField!(G, D, Unqual!(typeof(var(deref(gen)))))(gen, var);
 }
 
 /// ditto
 auto field(alias gen, D)(D var)
-    if (isSaturatedRandomEngine!(typeof(gen)))
+    if (isSaturatedRandomEngine!(typeof(gen)) && isRandomVariable!D)
 {
     return RandomField!(gen,D,Unqual!(typeof(var(gen))))(var);
 }
@@ -268,7 +269,7 @@ Range interface for random distributions and uniform random bit generators.
 Note: $(UL $(LI The structure holds a pointer to a generator.) $(LI The structure must not be copied (explicitly or implicitly) outside from a function.))
 +/
 struct RandomRange(G, D)
-    if (isReferenceToSaturatedRandomEngine!G)
+    if (isReferenceToSaturatedRandomEngine!G && isRandomVariable!D)
 {
     private D _var;
     private G _gen;
@@ -288,7 +289,7 @@ struct RandomRange(G, D)
 
 /// ditto
 struct RandomRange(alias gen, D)
-    if (isSaturatedRandomEngine!(typeof(gen)))
+    if (isSaturatedRandomEngine!(typeof(gen)) && isRandomVariable!D)
 {
     private D _var;
     private Unqual!(typeof(_var(gen))) _val;
@@ -364,21 +365,21 @@ struct RandomRange(alias gen)
 }
 
 /// ditto
-RandomRange!(G, D) range(G, D)(ref G gen, D var)
-    if (isReferenceToSaturatedRandomEngine!G)
+RandomRange!(G, D) range(G, D)(G gen, D var)
+    if (isReferenceToSaturatedRandomEngine!G && isRandomVariable!D)
 {
     return typeof(return)(gen, var);
 }
 
 /// ditto
-RandomRange!(gen, D) range(alias gen, D)(D var)
-    if (isSaturatedRandomEngine!(typeof(gen)))
+RandomRange!(gen, D) range(alias gen = rne, D)(D var)
+    if (isSaturatedRandomEngine!(typeof(gen)) && isRandomVariable!D)
 {
     return typeof(return)(var);
 }
 
 /// ditto
-RandomRange!G range(G)(ref G gen)
+RandomRange!G range(G)(G gen)
     if (isReferenceToSaturatedRandomEngine!G)
 {
     return typeof(return)(gen);
@@ -743,6 +744,20 @@ void shuffle(Range, G)(scope ref G gen, scope Range range)
     }
 }
 
+/// ditto
+void shuffle(Range, G)(scope G* gen, scope Range range)
+    if (isSaturatedRandomEngine!G && isRandomAccessRange!Range && hasLength!Range)
+{
+    return .shuffle(*gen, range);
+}
+
+/// ditto
+void shuffle(Range)(scope Range range)
+    if (isRandomAccessRange!Range && hasLength!Range)
+{
+    return .shuffle(rne, range);
+}
+
 ///
 nothrow @safe version(mir_random_test) unittest
 {
@@ -750,10 +765,9 @@ nothrow @safe version(mir_random_test) unittest
     import mir.ndslice.topology: iota;
     import mir.ndslice.sorting;
 
-    auto gen = Random(unpredictableSeed);
     auto a = iota(10).slice;
 
-    gen.shuffle(a);
+    shuffle(a);
 
     sort(a);
     assert(a == iota(10));
@@ -768,15 +782,14 @@ These will be in an undefined order, but will not be random in the sense that th
 `shuffle` returns will not be independent of their order before
 `shuffle` was called.
 Params:
-    gen = random number engine to use
+    gen = (optional) random number engine to use
     range = random-access range with length whose elements are to be shuffled
     n = number of elements of `r` to shuffle (counting from the beginning);
         must be less than `r.length`
 Complexity: O(n)
 +/
 void shuffle(Range, G)(scope ref G gen, scope Range range, size_t n)
-    if ((isSaturatedRandomEngine!G || isReferenceToSaturatedRandomEngine!G)
-        && isRandomAccessRange!Range && hasLength!Range)
+    if (isSaturatedRandomEngine!G && isRandomAccessRange!Range && hasLength!Range)
 {
     import std.algorithm.mutation : swapAt;
     assert(n <= range.length, "n must be <= range.length for shuffle.");
@@ -789,6 +802,20 @@ void shuffle(Range, G)(scope ref G gen, scope Range range, size_t n)
     }
 }
 
+/// ditto
+void shuffle(Range, G)(scope G* gen, scope Range range, size_t n)
+    if (isSaturatedRandomEngine!G && isRandomAccessRange!Range && hasLength!Range)
+{
+    return .shuffle(*gen, range, n);
+}
+
+/// ditto
+void shuffle(Range)(scope Range range, size_t n)
+    if (isRandomAccessRange!Range && hasLength!Range)
+{
+    return .shuffle(rne, range, n);
+}
+
 ///
 nothrow @safe version(mir_random_test) unittest
 {
@@ -796,10 +823,9 @@ nothrow @safe version(mir_random_test) unittest
     import mir.ndslice.topology: iota;
     import mir.ndslice.sorting;
 
-    auto gen = Random(unpredictableSeed);
     auto a = iota(10).slice;
 
-    gen.shuffle(a, 4);
+    shuffle(a, 4);
 
     sort(a);
     assert(a == iota(10));
diff --git a/source/mir/random/package.d b/source/mir/random/package.d
index 375319f4..b9a96793 100644
--- a/source/mir/random/package.d
+++ b/source/mir/random/package.d
@@ -85,6 +85,7 @@ T rand(T, G)(scope ref G gen)
         return cast(T) gen();
     }
 }
+
 /// ditto
 T rand(T, G)(scope G* gen)
     if (isSaturatedRandomEngine!G && isIntegral!T && !is(T == enum))
@@ -92,6 +93,20 @@ T rand(T, G)(scope G* gen)
     return rand!(T, G)(*gen);
 }
 
+/// ditto
+T rand(T)()
+    if (isIntegral!T && !is(T == enum))
+{
+    return rand!T(rne);
+}
+
+///
+@nogc nothrow @safe version(mir_random_test) unittest
+{
+    auto s = rand!short;
+    auto n = rand!ulong;
+}
+
 ///
 @nogc nothrow pure @safe version(mir_random_test) unittest
 {
@@ -113,6 +128,7 @@ bool rand(T : bool, G)(scope ref G gen)
     import std.traits : Signed;
     return 0 > cast(Signed!(EngineReturnType!G)) gen();
 }
+
 /// ditto
 bool rand(T : bool, G)(scope G* gen)
     if (isSaturatedRandomEngine!G)
@@ -120,6 +136,18 @@ bool rand(T : bool, G)(scope G* gen)
     return rand!(T, G)(*gen);
 }
 
+/// ditto
+bool rand(T : bool)()
+{
+    return rand!T(rne);
+}
+
+///
+@nogc nothrow @safe version(mir_random_test) unittest
+{
+    auto s = rand!bool;
+}
+
 ///
 @nogc nothrow pure @safe version(mir_random_test) unittest
 {
@@ -170,6 +198,7 @@ T rand(T, G)(scope ref G gen)
         return members[gen.randIndex($)];
     }
 }
+
 /// ditto
 T rand(T, G)(scope G* gen)
     if (isSaturatedRandomEngine!G && is(T == enum))
@@ -177,6 +206,20 @@ T rand(T, G)(scope G* gen)
     return rand!(T, G)(*gen);
 }
 
+/// ditto
+T rand(T)()
+    if (is(T == enum))
+{
+    return .rand!T(rne);
+}
+
+///
+@nogc nothrow @safe version(mir_random_test) unittest
+{
+    enum A { a, b, c }
+    auto e = rand!A;
+}
+
 ///
 @nogc nothrow pure @safe version(mir_random_test) unittest
 {
@@ -344,6 +387,7 @@ T rand(T, G)(scope ref G gen, sizediff_t boundExp = 0)
     /// TODO: quadruple
     else static assert(0);
 }
+
 /// ditto
 T rand(T, G)(scope G* gen, sizediff_t boundExp = 0)
     if (isSaturatedRandomEngine!G && isFloatingPoint!T)
@@ -351,6 +395,32 @@ T rand(T, G)(scope G* gen, sizediff_t boundExp = 0)
     return rand!(T, G)(*gen, boundExp);
 }
 
+/// ditto
+T rand(T)(sizediff_t boundExp = 0)
+    if (isFloatingPoint!T)
+{
+    return rand!T(rne, boundExp);
+}
+
+
+///
+@nogc nothrow @safe version(mir_random_test) unittest
+{
+    import mir.math.common: fabs;
+    
+    auto a = rand!float;
+    assert(-1 < a && a < +1);
+
+    auto b = rand!double(4);
+    assert(-16 < b && b < +16);
+    
+    auto c = rand!double(-2);
+    assert(-0.25 < c && c < +0.25);
+    
+    auto d = rand!real.fabs;
+    assert(0.0L <= d && d < 1.0L);
+}
+
 ///
 @nogc nothrow pure @safe version(mir_random_test) unittest
 {
@@ -371,7 +441,6 @@ T rand(T, G)(scope G* gen, sizediff_t boundExp = 0)
     assert(0.0L <= d && d < 1.0L);
 }
 
-
 /// Subnormal numbers
 @nogc nothrow pure @safe version(mir_random_test) unittest
 {
@@ -410,9 +479,10 @@ Params:
 Returns:
     Uniformly distributed integer for interval `[0 .. m$(RPAREN)`.
 +/
-T randIndex(T, G)(scope ref G gen, T m)
+T randIndex(T, G)(scope ref G gen, T _m)
     if(isSaturatedRandomEngine!G && isUnsigned!T)
 {
+    immutable m = _m + 0u;
     static if (EngineReturnType!G.sizeof >= T.sizeof * 2)
         alias MaybeR = EngineReturnType!G;
     else static if (uint.sizeof >= T.sizeof * 2)
@@ -507,6 +577,7 @@ T randIndex(T, G)(scope ref G gen, T m)
         return u.high;
     }
 }
+
 /// ditto
 T randIndex(T, G)(scope G* gen, T m)
     if(isSaturatedRandomEngine!G && isUnsigned!T)
@@ -514,6 +585,20 @@ T randIndex(T, G)(scope G* gen, T m)
     return randIndex!(T, G)(*gen, m);
 }
 
+/// ditto
+T randIndex(T)(T m)
+    if(isUnsigned!T)
+{
+    return randIndex!T(rne, m);
+}
+
+///
+@nogc nothrow @safe version(mir_random_test) unittest
+{
+    auto s = randIndex(100u);
+    auto n = randIndex!ulong(-100);
+}
+
 ///
 @nogc nothrow pure @safe version(mir_random_test) unittest
 {
@@ -577,6 +662,7 @@ size_t randGeometric(G)(scope ref G gen)
         if(auto val = gen.rand!T())
             return count + bsf(val);
 }
+
 /// ditto
 size_t randGeometric(G)(scope G* gen)
     if(isSaturatedRandomEngine!G)
@@ -584,18 +670,25 @@ size_t randGeometric(G)(scope G* gen)
     return randGeometric!(G)(*gen);
 }
 
-@nogc nothrow pure @safe version(mir_random_test) unittest
+/// ditto
+size_t randGeometric()()
 {
-    import mir.random.engine.xorshift;
-    auto gen = Xoroshiro128Plus(1);
-    size_t s = gen.randGeometric;//Merely verify the call is @safe etc.
+    return randGeometric(rne);
 }
 
+///
 @nogc nothrow @safe version(mir_random_test) unittest
 {
-    //Coverage. Impure because uses thread-local.
-    Random* gen = threadLocalPtr!Random;
-    size_t s = gen.randGeometric;//Merely verify the call is @safe etc.
+    size_t s = randGeometric;
+}
+
+///
+@nogc nothrow pure @safe version(mir_random_test) unittest
+{
+    import mir.random.engine.xorshift;
+    auto gen = Xoroshiro128Plus(1);
+
+    size_t s = gen.randGeometric;
 }
 
 /++
@@ -649,6 +742,7 @@ T randExponential2(T, G)(scope ref G gen)
     else
         return -log2(x) + y;
 }
+
 /// ditto
 T randExponential2(T, G)(scope G* gen)
     if (isSaturatedRandomEngine!G && isFloatingPoint!T)
@@ -656,18 +750,24 @@ T randExponential2(T, G)(scope G* gen)
     return randExponential2!(T, G)(*gen);
 }
 
+/// ditto
+T randExponential2(T)()
+    if (isFloatingPoint!T)
+{
+    return randExponential2!T(rne);
+}
+
 ///
 @nogc nothrow @safe version(mir_random_test) unittest
 {
-    import mir.random.engine.xorshift;
-    auto gen = Xorshift(cast(uint)unpredictableSeed);
-    auto v = gen.randExponential2!double();
+    auto v = randExponential2!double;
 }
 
-@nogc nothrow @safe version(mir_random_test) unittest
+///
+@nogc nothrow @safe pure version(mir_random_test) unittest
 {
-    //Coverage. Impure because uses thread-local.
-    Random* gen = threadLocalPtr!Random;
+    import mir.random.engine.xorshift;
+    auto gen = Xorshift(1);
     auto v = gen.randExponential2!double();
 }
 
@@ -742,6 +842,7 @@ struct PhobosRandom(Engine) if (isRandomEngine!Engine && !isPhobosUniformRNG!Eng
         return _engine;
     }
 }
+
 /// ditto
 template PhobosRandom(Engine) if (isRandomEngine!Engine && isPhobosUniformRNG!Engine)
 {