Add randomSubset method

Author: d-torrance

M2/Macaulay2/m2/exports.m2

@@ -1045,6 +1045,7 @@ export {
 	"random",
 	"randomKRationalPoint",
 	"randomMutableMatrix",
+	"randomSubset",
 	"rank",
 	"rays",
 	"read",

M2/Macaulay2/m2/lists.m2

@@ -237,6 +237,34 @@ random List := opts -> s -> (
 	  t := s#i ; s#i = s#j ; s#j = t;
 	  );
      new List from s)
+
+randomSubset = method()
+-- Knuth Algorithm S, Art of Computer Programming, Section 3.4.2
+randomSubset(ZZ, ZZ) := (N, n) -> (
+    if n < 0 or n > N then error("expected an integer between 0 and ", N);
+    t := 0;
+    apply(n, m -> (
+	    while (N - t) * rawRandomRRUniform defaultPrecision >= n - m
+	    do t += 1;
+	    first (t, t += 1))))
+-- when size not given, use binomial distribution first to determine size
+randomSubset ZZ := N -> (
+    if N < 0 then error "expected a nonnegative integer"
+    else if N < 30 then (
+	-- for small N, sum up N Bernoulli variates
+	n := 0;
+	scan(N, i -> n += random 2))
+    else (
+	-- for large N, use normal approximation
+	n = round(n/2 + sqrt(n/4)*rawRandomRRNormal defaultPrecision);
+	-- in the vanishingly unlikely case we end up outside [0, N]:
+	n = min(N, max(0, n)));
+    randomSubset(N, n))
+randomSubset(VisibleList, ZZ) := (x, n) -> x_(randomSubset(#x, n))
+randomSubset VisibleList := x -> x_(randomSubset(#x))
+randomSubset(Set, ZZ) := (x, n) -> set randomSubset(toList x, n)
+randomSubset Set := x -> set randomSubset toList x
+
 -----------------------------------------------------------------------------
 -- sublists
 -----------------------------------------------------------------------------

M2/Macaulay2/packages/Macaulay2Doc/functions.m2

@@ -219,6 +219,7 @@ load "./functions/pushForward-doc.m2"
 load "./functions/QRDecomposition-doc.m2"
 load "./functions/quotient-remainder-doc.m2"
 load "./functions/random-doc.m2"
+load "./functions/randomSubset-doc.m2"
 load "./functions/rank-doc.m2"
 load "./functions/readDirectory-doc.m2"
 load "./functions/readlink-doc.m2"

M2/Macaulay2/packages/Macaulay2Doc/functions/random-doc.m2

@@ -23,6 +23,7 @@ Node
     (random, Module)
     (random, Module, Module)
     setRandomSeed
+    randomSubset
 
 Node
   Key

M2/Macaulay2/packages/Macaulay2Doc/functions/randomSubset-doc.m2

@@ -0,0 +1,40 @@
+doc ///
+  Key
+     randomSubset
+    (randomSubset, VisibleList, ZZ)
+    (randomSubset, Set, ZZ)
+    (randomSubset, ZZ, ZZ)
+    (randomSubset, VisibleList)
+    (randomSubset, Set)
+    (randomSubset, ZZ)
+  Usage
+    randomSubset(x, n)
+    randomSubset x
+  Inputs
+    x:{VisibleList, Set, ZZ}
+    n:ZZ
+  Description
+    Text
+      When @VAR "n"@ is given, then a random subset @VAR "x"@ of cardinality
+      @VAR "n"@ is returned.
+    Example
+      randomSubset({2, 3, 5, 7, 11}, 2)
+    Text
+      Otherwise, a random subset of arbitrary cardinality is returned.
+    Example
+      randomSubset {2, 3, 5, 7, 11}
+      randomSubset {2, 3, 5, 7, 11}
+      randomSubset {2, 3, 5, 7, 11}
+    Text
+      If @VAR "x"@ is an integer, then a subset of $\{0,\ldots,x - 1\}$ is
+      returned.
+    Example
+      randomSubset(4, 2)
+      randomSubset 5
+  References
+    Knuth, Donald E.
+    @EM "The Art of Computer Programming: Seminumerical Algorithms, Volume 2"@.
+    Addison-Wesley Professional, 2014. (Algorithm S, Section 3.4.2)
+  SeeAlso
+    random
+///