1 files changed, 58 insertions, 0 deletions
diff --git a/src/misc.c b/src/misc.c
index 183f0ce..88e4b77 100644
--- a/src/misc.c
+++ b/src/misc.c
@@ -684,3 +684,61 @@ void unique_vertices(int *id, size_t n, size_t npeaks, size_t *result, size_t *n
     free(ordered_results);
     free(unique_results);
 }
+
+static int is_sorted(size_t *a, size_t n)
+{
+    size_t i;
+
+    for (i = 1; i < n; i++)
+        if (a[i] < a[i-1]) return 0;
+
+    return 1;
+}
+
+void combinations_with_replacement(size_t n, size_t r, size_t *result, size_t *len)
+{
+    /* Returns the set of all unique combinations of length `r` from `n` unique
+     * elements.
+     *
+     * The result array can be indexed as:
+     *
+     *     result[i*r + j]
+     *
+     * where i is the ith combination and j is the jth element in the
+     * combination. For example:
+     *
+     *     size_t result[12];
+     *     size_t len;
+     *     product(3,2,result,&len);
+     *     for (i = 0; i < len; i++)
+     *         print("%zu %zu\n", result[i*2], result[i*2+1]);
+     *
+     * will print
+     *
+     *     0 0
+     *     0 1
+     *     0 2
+     *     1 1
+     *     1 2
+     *     2 2
+     *
+     * `result` should be an array with at least (n+r-1)!/r!/(n-1)! elements. */
+    size_t i, j;
+    size_t *tmp;
+
+    tmp = malloc(sizeof(size_t)*n*ipow(r,n));
+
+    product(n,r,tmp);
+
+    *len = 0;
+    for (i = 0; i < ipow(n,r); i++) {
+        if (is_sorted(tmp+i*r,r)) {
+            for (j = 0; j < n; j++) {
+                result[(*len)*r+j] = tmp[i*r+j];
+            }
+            *len += 1;
+        }
+    }
+
+    free(tmp);
+}