An Optimal Parallel Algorithm for
Generating Combinations
Selim G. Aki*
David Gries**
Ivan Stojmenovic***
IR 89-1029
August 1989
Department of Computer Science
Cornell University
Ithaca, NY 14853-7501
*Computing and Information Science, Oueen's tiniversity, Kingston, Ontario, Canada K7L 3N6.
**Computer Science, Upson Hall, Cornell University, Ithaca, NY 14853-7501, USA.
***Computer sdenoe, Unrversity of Ofawa, 0rrtar?, Canada Ki N 9B4.
An Optimal Parallel Algorithm for Generating Combinations
Selim G. Aki1, David Cries2, Ivan Stojmenovic3
Key words: concurrency, design of algorithms, parallel algorithms, generating combinations.
Introduction
We present a parallel algorithm for generating in lexicographically ascending order the C (m, n) combi-
nations of m objects chosen from l..n,where4 1 <?m <n and C(m, n) is the total number of combi-
nations. (Henceforth, we write C (rn, n) as Cmn.) We call these combinations m-combinations. The
algorithm is designed to run on a linear array of m processors, where each processor communicates only
with its (left and right) neighbors, each has constant size memory, and each is responsible for producing
one element of each combination. A combination takes constant time to produce from the preceding one,
so the algorithm requires time 0 (Cmn). This improves over previous algorithms in three respects: it
runs on a weaker model, it is cost optimal (assuming the time to output the combinations is counted), and
it requires the calculation of no integer greater than n.
Comparisons of the various sequential ?Agorithms of [8, 10, 5, 9, 13, 12] for this problem are given in
[11, 3]. In the sequential setting, since each of the Cmn combinations requires 0(m) time to produce
as output, an optimal algorithm requires time 0 (m*Cmn). However, there are sequential algorithms that
require only time 0 (Cmn) if the combinations need not actually be output.
In [6, 7], parallel algorithms are given to generate the combinations of at most m out of n objects.
However, when applied to m-combinations, neither algorithm is cost-optimal (i.e. the number of proces-
sors times the time complexity does not match the optimal sequential running time). The algorithm in [2,
1] uses an arbitrary number of independent processors, each producing an interval of consecutive combi-
nations. The algorithm is cost-optimal, but each processor requires memory of size 0(n) and has to deal
with large integers.
Finally, [4, 1] describe an 0 (Cmn*Iog.n)-time algorithm for generating m-combinations on an
Exclusive-Read Exclusive-Write Parallel Random Access Machine (EREW PRAM). In this model, a
number of processors share a common memory and execute the same instruction (on different data) syn-
chronously, but no two processors may access the same memory location simultaneously. The algorithm
of [4, 1] uses m processors and a memory of size 0(m). In contrast, our algorithm requires time
o (Cmn) and runs on a weaker model, the linear array of processors, and any aIgorithm running on a
linear array of processors can be simulated without penalty on the EREW PRAM [cf. 1].
Computing and Information Science, Qucen's University. Kingston. 0ntario, Canada K7L 3N6
2 Computer Science, Upson hall. Cornell University. Ithaca, NY l4853?7SOl, USA.
34 Computer Science. University7of 0ttowa, 0ntario, Canada KIN 9B4.
In the case m = n there is exactly one combination. For technical reasons, this case is different from the rest, and. in the interests of simplicity
and space conservation, it is omitted completely from the discussion.
2
2. The algorithm
We denote the ascending sequence of m-combinations of 1..n by d. 1 , d. 2d.Cmn, so d.j is an
ascending sequence of integers d.j. 1d.j.m. [Note, the period denotes application of a function,
array, or sequence to a value.] Fig. 1 contains the C (3, 5) combinations, in lexicographical order.
d.1;			1 2			3
d.2:			l 2			4
d.3:			1 2			5
d.4:			1 3			4
d.5:			1 3			5
d.6:			1 4			5
d.7:			2 3			4
d.8:			2 3			5
d.9:			2 4			5
d.lO:3 4			5
123
124
//
125
//
134
135
145
27/4
235
245
/ /
345
34
3/
Fig. 1			Fig. 2
Consider an algorithm to generate the m-combinations using m processors and m shared variables
D1Dm , whose meaning is to be decided upon laten Each processor i is responsible for maintaining
Dj and, we hope, reads only Dj?1 and Dj+1 . Each processor i has a while loop that at each iteration
prints d.j.i so that, since the processors act in lock-step fashion, each iteration will print d.j. We might
naively expect that processor i would maintain Dj = d.j.i. Thus, the body of the loop of each processor
would be
Print D?;
Dj :=d.(j+1).i;
j :=]+1
(using the fact that Dj = d.i.i)
However, the computation of d. kj+1).t would require too much time, too much space, or too many
references to shared memory. We therefore look for a slight variation on this theme.
It will help for processor i be one step ahead of processor i+l in the calculation, as shown in Fig. 2.
A sequence of lines connecting values indicates the values that are in the variables Dj at each step. For
example, initially, we have D1 = d.3.1 , D2 = d.2.2, and D3 = d. 1.3. This means that when Dm is
d.Cmn.m , processors 1..m--Hl will have calculated further values, as shown in Fig. 2. For later pulposes,
these further values, as well as d.j. 0 and d.j. (m +1), are defined as follows.
d.j.O=n--Hm,			for l<?j,
di. (m+l)=n			for I <--Hi.
d.j.i = d.Cmn.i = n--Hm+i
(1)
for I > Cmn and l < i < m
3
Thus, we give the first invariant PO of the loops of processors 1..m; in doing so, we introduce two
additional shared variabics D0 and Dm+i
PO: l<j<?Cmn+l and
the m-combinations d. 1d q--Hl) are printed and
Dj = d.?+m--Hi)J for 0 <i < m+l
The key to our algorithm is the efficient calculation of each Dj at each step, when j is increased by
1. Just how this calculation can be accomplished is given in
Lemma 1. The following formula describes d. (j+l).i for 1 < i < m:
%t$i=n--Hm+i			? d.(j+l).q--Hl)+l
d.(j+l).i =			d.fi <n--Hm+i A T,i.i			? dii +1
#f? <n--Hm+i A ?Tj.i ? d.j.i
where T.j.i = "value n--Hm+i is in position i+l of the m-combinations d. (i--Hl)d. (i--H(m?))":
T.j.i = V(k: 1 <k <m--Hi: d.(j--Hk).(i+l) = n--Hm+i),
Proof Consider the first case, d.j.i = n--Hm+i If j > Cmn, the formula follows from (1): d. (i+l).i =
n--Hm+i = d. (j+1).(i--Hl)+l . Suppose that j < Cmn. Since d.j is an ascending sequence of values from
l..n we have
dl. (im) = n--Hm+i, n--Hm+i+ln--Hl, n.
Hence, d. (i+l).i is at most d.j.i, and since it is as small as possible (so that the m-combinations are in
lexicographic order), it equals d. (i+l).(i--Hl)+l
Now consider the case d.j.i = v (say) where v < n--Hm+i By (1), I < Cmn. Let d.j. (l..i) = X v
for some sequence X of length i--Hi , where catenation of sequences and elements is denoted by juxtapo-
sition. Suppose j is the last combination (in lexicographic order) with prefix X v. We list below the
last of the combinations (in lexicographic order) with this prefix , as well as the following combination
d.(i+l):
d.(j--Ht)			=			x			v			n-m+i			Y1
d.kj--Hl)			=			x			v			n--Hm+i			Y,
d.j			=			x			v			n--Hm+i+l n-m+i+2
=			X			v+l			v+2
for some sequences Y1 ..., Y1. Note that n--Hm+i+l is the highest value that can be placed in position
i+l of a combination. Hence, in combination d. (j+l), we have v+l in position i, as shown. Also as
shown, there are t (say) previous combinations with n--Hm+i in position 1+1 . What is t? Since the
sequences Y1Y1 consist of m--Hi--Hl integers chosen from the m--Hi integers n--Hm+i+l..n, the
numberofsequencesis t= C(m--Hi,m--Hi--Hl) =m--Hi
Any combination with prefix X v, except the last, is not preceded by m --Hi combinations with
n--Hm+i in position i+l Hence, Lemma 1 is correct for the case d.j.i < n--Hm+i. El
4
Invariant PO defines the shared variables Do..m+i , with each processor i, 1 < i < m, maintaining
Dj and D0 and Dm+i remaining constant. We now determine other values that must be maintained by
processor i and introduce an efficient method of maintaining them. First, at iteration ?, pmeessor
prints d.j.i. However, in variable Dj is the value d. j+m--Hi).i. In order to print dj.i , pmcessor i will
maintain the whole range of values d. (j..j+m?).i
Second, we investigate what processor i+l may require from processor i in order to calculate its
next value Dj+1 at each iteration. By PO, since j is to be increased by 1, this next value is
d. (j+l+m--H(i+1)).(i+1) , i.e. d. (j+m--Hi).(i+l). By Lemma 1, this calculation may require reference to
d. (j+m--Hi).i , so processor i will maintain it.
Third, we investigate what processor i--Hi may require from processor i in order to calculate its next
value Dj?j . By PO , this next value is d. j+l+m--H(i--Hl)).(i--Hl). Evaluation of this may require evalua-
tion of T. j+m--Hi+l).(i--Hl) of Lemma 1, which references the values d. ?+l.j+m--Hi).i. Therefore,
processor t will maintain these values.
Combining the results of the last three paragraphs, we see that processor i must maintain
d. (I. j+m --Hi).i
Luckily, it turns out that d. ?..j+m--Hi).i contains at most three distinct values, for the following rea-
son. ?t position i, any element v, v < n --Hm +i, occurs C (n --Hv, mi) times in successive combina-
tions. For v = n--Hm+i, the number of occurrences in succession is only 1, but if v < n--Hm+t,
C(n--Hv, m--Hi) > m--Hi+l Hence, within any m--Hi+l successive combinations, position I has at most
three different values.
Processor I can maintain d. kj.j+mA).I in five local variables Dj , r? ,dij , r1? and d2j , as shown
in the following invariant Pi , the second invariant of the loop for processor I. It is now straightforward
to verify that the body of the while loop in the algorithm below maintains loop invariants PO and Pi
Pi: Th?dli?d2i and
d. (j..j+m--HI).I is a suffix of the sequence
d2j , (dlj repeated ri? times), (Dj repeated r? times)
Finally, we determine the loop-termination condition.
Lemma 2.? All the elements of d. (Cmn--Hi+l..Cmn).i equal n--Hm+i. This is the only place within the
Cmn combinations that the element n--Hm+i occurs I times in succession in position I
Proof The proof is by induction on I. For I = 1 , the last combination is n--Hm+l..n, and since each
combination is an ascending sequence of integers, n--Hm+l cannot appear in the first position of any
other combination.
We now assume the lemma is true for I < m and prove it true for 1+1 We first show that all ele-
ments of d. (Cmn--Hi..Cmn).(i+l) equal n--Hm+i+l . By the inductive assumption, d.j.i = n--Hm+i for I
in Cmn--Hi+l..Cmn. Hence df ??.n) = n--Hm+I..n, and df (1+1) = n--Hm+I+l for these values of I.
We now prove that J (Cmn--Hi).(i+l) = n--Hm+I+l . By the inductive assumption, the fact that n--Hm+i
)? ?L
temma 2 does not hold for n = m this is why we require s <?.
5
is the largest element that can appear in position i , and the lexicographic order, we have d. (Cmn--Hi).i =
n--Hm+i--Hl. Thus d.(Cmn--Hi).(i+l) > n--Hm+i, for if d.Cmn.(i+l) = n--Hm+i, then we have an unac-
counted for combination between d. (Cmn--Hi) and d. (Cmn--Hi+l), with the same first i elements as
d. (Cmn--Hi) and with n?n+i+l in position i+1 , which is absunl. Since d. (Cmn--Hi).(i+l) cannot be
larger than n ?n+i+1, our claim is established.
We now show that d. (Cmn--Hi..Cmn) is the only place within the sequence of combinations in which
n--Hm+i+l occurs i+l times in position i+l. Assume that the n?+i+l occurs in position i+l for
a sequence of i+l combinations not equal to (but that may partly overlap) the sequence
d. (Cmn--Hi..Cmn). Let d.j be any one of the first i combinations of this new sequence. Since
dj.(i+l) = d.(i+l).(i+l) = n?+i+l , it follows by Lemma 1 that d.(j+l).i = n--Hm+i. This means
that n--Hm+i occurs in position i for a sequence of i combinations other than d. (Cmn--Hi+l..Cmn).
This contradicts the inductive assumption. 0
Termination should occur when j = Cmn+l, for then, by PO , all the combinations have been printed.
For this value of j, by (1), all the elements of d. (I.j+m--Hi).i equal n--Hm+i. By Lemma 2, all the ele-
ments of d. (i--Hi..j--Hl).i equal n--Hm+i. Hence, the m+l elements of d. kj--Hi.j+m--Hi).i all equal
n--Hm+i. Further, by Lemma 2, for no other value of j in l..Cmn+l does this hold. Finally, by Pi
this condition can he written as Dj = n--Hm+i and r? > m.
We thus arrive at the following algorithm.
D0 :=n?; Dm+i,rm+i :=n,O;
forparallel i: 1 S i ? m do
Dt, r? := i, `n--Hi+l;
j := 1;
do ?(Dj = n--Hm+i and r? > m) ?
"Printd.j.i;" if4 > m--Hi then print (Dj)
else if 4+rij > m--Hi then print (dij)
else print ((t2?);
"Establish (PO and PJ)j?+1 :"
varx,y, ,y:int;
x := D??1; y := D?+1; ry := r?+1; (Note: 3 references to shared variables)
ifDj = n--Hm+i then ifDj = x+l then rj := 4+1
elsed2? :=d1?;
dij, rJ? := Dj, 4;
Dj, 4 :=x+l, 1
else if y = n--Hm+i A ry = m--Hi then d2j := d11;
d11, r11 :=D1, 4;
D1,4 :=D1+l, 1
else 4 := 4+1;
j :=j+l
od
6
Two simple modifications can be made to the algorithm. First, variable j is not used in a substantial
way; it is a 4'thought" variable, used only in discussing the algorithm, and may be eliminated. Second,
note that r? and r1? may become quite large. They both can be reduced to the range O..m+1 , without
effecting the algorithm in any other way, by replacing the statement r? r?+1 by if r? < m
then r? r?+1. The conjunct r1 > m was used in the loop condition, instead of r? = m+l , in order to
make this change possible.
Conclusion
We derived an 0 (Cmn)) algorithm for generating all combinations of m out of n objects in lexico-
graphic order on an EREW PRAM and on a linear array of m procesors, each using constant-space local
memory. The algorithin produces the combinations in constant time per combination and requires no
integer larger than n during the computation. This improves on the previously known algorithms and is
an optimal-cost algorithin, assuming that the time to output combinations is counted.
The algorithm can be made adaptive (i.e. to run on a linear array consisting of an arbitrary number k
of processors) if the array is divided into kim groups of m processors each such that each group pro-
duces an interval of consecutive m-combinations, as described in tl,2].
The problem remains open to fmd an analagous optimal algorithm for enumerating permutations of m
out of n objects.
References
[ ii S.G. Aki. The Design and Analysis ofParallel Algorithms. Prentice Hall, Englewood Cliffs. New
Jersey, 1989.
[ 2j 5.0. AkI. Adaptive and optimal parallel algorithms for enumerating permutations and combina-
tions. Computeri. 30, 5 (1987), 433436.
( 3] 5.0. Aki. A comparison of combination generation methods. ACM Trans. Mathematical Software
7, 1(1981), 4245.
[ 4] B. Chan and 5.6. Aki. Generating combinations in parallel. BIT 26, 1(1986), 2-6.
5] P.J. Chase. Combinations of m out of n objects (Mgorithm 382). Comm. ACM 13, 6 (1970), 368.
( 6] G.H. Chen and M.S. Chem. Parallel generation of permutations and combinations. BiT 26 (1986),
277-283.
[ 7] B. Djokic, M. Miyakawa, 5. Sekiguchi, I. Semba, and I. Stojmenovic. Parallel algorithms for gen-
erating subsets and set partitions. Computeri., to appear.
[ 8] J. Kutzberg. Combination (Algorithin 94). Comm. ACM 5, 6 (1962), 344.
[ 9] C.M. Liu and D.T. Tang. Enumerating combinations of m out of n objects (Algoritllm 452).
Comm. ACM 16, 8 (1973), 485.
[10] C.S. Misfud. Combination in lexicographic order (Algorithm, 154). Comm. ACM 6, 3 (1963), 103.
[11] W.H. Payne and F.M. Ives. Combination generators. ACM Trans. Mathematical Software 5, 2
(1972), 163-172.
[12] F. Ruskey. Adjacent interchange combinations. J. Algorithms 9 (1988), 162-180.
[13] I. 5emb? A note on enumerating combinations in lexicographic order. J. Information Processing
4, 1(1981), 35-37.