:

## Cartesian Products of Lists

The Cartesian product of a sequence of m lists L1, ..., Ln, with list i having ni elements, is given by the sequence of Πi ni lists [L11,...,Lm1], ..., [L1n1,...,Lmnm], where Ljk is the k-th element of the j-th list.

If the right-most element in the lists varies the fastest, the sequence is in row-major order, if the left-most elements varies the fastest, the sequence is in column-major order. With one exception (noted) all the routines given below use row-major order.

Because a Cartesian product can be large, Maple provides a procedure, combinat[cartprod], that generates an iterator for the Cartesian product. Each call to this iterator returns the next value in the sequence. Here is a procedure that uses cartprod to iterate through each list in the Cartesian product of its arguments, passing each generated list to the function f.

Note: This procedure, as well as most of the following, uses the seq parameter modifier, a mechanism introduced in Maple 11. The parameter L is assigned an expression sequence consisting of the following lists. Here it is equivalent to using the older args[2..-1], provided the remaining arguments are lists.

```use_combinat := proc(f, L::seq(list))
local T;
uses combinat;
T := cartprod([L]);
while not T['finished'] do
f(T['nextvalue']());
end do;
NULL;
end proc:
```

By assigning f to print its argument, we can demonstrate that this works as intended:

```f := curry(printf,"  %a"):
lsts := [\$1..3],[a,b]:
use_combinat(f,lsts): printf("\n"):
[1, a]  [1, b]  [2, a]  [2, b]  [3, a]  [3, b]
```

Using combinat[cartprod] is a straight-forward and effective way to step through the entire set without having to allocate memory for the set. The iterator itself, however, requires some time to execute and if the loop itself is otherwise fast—say it calls a compiled procedure—then the time spent evaluating the iterator can be significant.

The Maple built-in procedure rtable_scanblock applies a function to each index of an rtable. If each dimension of the rtable is a range that corresponds to the number of elements in each of the lists, then there is a simple mapping from the indices of the rtable to the sublists in the Cartesian product. We can use this, then, to effectively iterate over the Cartesian product. Because the rtable is created with the storage=sparse option, the memory allocation is minimal. Here is a test procedure that implements this technique.

```use_scanblock := proc(f)
local a,dims,n,L;
L := args[2..-1]:
n := nargs-1;
dims := seq(1..nops(a), a in L);
rtable_scanblock(rtable(dims, 'storage=sparse')
, [dims]
, proc(val,indx,data)
local i;
f([seq(L[i][indx[i]], i=1..n)]);
end proc
);
NULL;
end proc:
```

```use_scanblock(f,lsts): printf("\n"):
[1, a]  [1, b]  [2, a]  [2, b]  [3, a]  [3, b]
```

The following tabulates, for each procedure, the time and memory required to iterate through each element of a Cartesian product of lists. Four different inputs are tested. In each case, the rtable_scanblock technique is faster (about 2×) and uses less memory than combinat[cartprod].

```Iterate through the Cartesian product of N lists of E elements,
generating E^N iterations.  At each step, call the (unassigned)
global function f with the generated list.

gcfreq = 10^7

+-------- timing -------+--------- memory ---------+
256^2          | maple (s) | total (s) | used (MiB) | alloc (MiB) | gctimes
---------------+-----------+-----------+------------+-------------+--------
use_combinat   |     1.77  |     1.90  |    29.98   |     23.75   |    0
use_scanblock  |     0.93  |     0.95  |    19.23   |     13.87   |    0
---------------+-----------+-----------+------------+-------------+--------
10^5           |           |           |            |             |
---------------+-----------+-----------+------------+-------------+--------
use_combinat   |     3.38  |     3.47  |    47.84   |     31.74   |    1
use_scanblock  |     1.84  |     1.93  |    28.49   |     22.68   |    0
---------------+-----------+-----------+------------+-------------+--------
3^10           |           |           |            |             |
---------------+-----------+-----------+------------+-------------+--------
use_combinat   |     2.73  |     2.81  |    37.70   |     31.49   |    0
use_scanblock  |     1.40  |     1.45  |    21.51   |     16.06   |    0
---------------+-----------+-----------+------------+-------------+--------
2^16           |           |           |            |             |
---------------+-----------+-----------+------------+-------------+--------
use_combinat   |     3.99  |     4.13  |    49.50   |     32.62   |    1
use_scanblock  |     2.06  |     2.08  |    25.63   |     21.00   |    0
---------------------------------------------------------------------------
```

While iterators are useful for reducing the memory footprint, at times it is necessary to instantiate the Cartesian product as a single entity. Maple, alas, does not provide a function for doing so (this might make be a nice addition to the ListTools package). It is not difficult to create a procedure, the interesting part is making it fast. The following site has a few examples of implementations; I have used them below with minor modifications, along with some of my own contributions.

The first example uses combinat[cartprod], called Πi ni times, to generate the Cartesian product.

```CartProdDeiterated := proc(L::seq(list))
local C,i,S;
C := combinat['cartprod']([L]);
[seq(C['nextvalue'](), i=1..mul(nops(S), S=[L]))]
end proc:
```

The next example uses a modification of the rtable_scanblock technique. Instead of calling rtable_scanblock, it generates an rtable with an initializer that maps each index to the corresponding entry in the Cartesian product. The elements of the rtable are then converted to a list by using seq.

```CartRtable := proc(L::seq(list))
local i,a,n;
n := nargs;
[seq(a, a in rtable(seq(1..nops(a), a in args)
,() -> [seq(L[i][args[i]],i=1..n)]
,'order = C_order'
))];
end proc:
```

Here is an offering from Carl DeVore (aka Carl Love). I've modified it slightly to use the seq parameter modifier previously mentioned. This routine produces the Cartesian product in column-major order; that is, the left-most element in the lists varies the fastest. All other routines produce the elements in row-major order.

```CartProdAsList := proc(L::seq(list))
local q,i,N;
N := map(nops, [L]);
[seq([seq(L[i][1+irem(q,N[i],'q')], i= 1..nargs)], q=0..`*`(N[])-1)]
end proc:
```

Here is a suggestion from Helmut Kahovec, also slightly modified.

```CartProdRecursive := proc(L::seq(list))
local i,j;
if nargs=2 then
[ seq(seq([`if`(i::list,i[],i),
`if`(j::list,j[],j)
] , j in args )
, i in L ) ];

elif nargs>2 then
procname( procname(args[1..2]), args[3..-1] );
else
args;
end if;
end proc:
```

The previous code has a few minor flaws and permits a few optimizations. Here is a non-recursive version that corrects the flaws and improves its performance.

```CartProdNonrecursive := proc(L::seq(list))
local i,j,X,lst;
if nargs=0 then return [];
elif nargs=1 then return map(`[]`,L);
end if;
X := [seq(seq([i,j], j in L ), i in L )] ;
for lst in L[3..-1] do
X := [seq(seq([i[],j], j in lst ), i in X )];
end do;
X;
end proc:
```

Here is an interesting version that I originally devised as a way to write an add-like procedure that summed over multiple dimensions, converted here to use seq to generate the Cartesian product. Because foldl is not a built-in, it is possible to improve it further, however, the resulting speed gains were trivial (less than five percent).

```CartProdSeq := proc(L::seq(list))
local Seq,i,j;
eval([subs(Seq=seq, foldl(Seq
, [cat(i,1..nargs)]
, seq(cat(i,j)=L[j],j=nargs..1,-1)
))]);
end proc:
```
To better understand what it does, run it without the eval statement, or do
```test := subs(eval=''eval'',op(CartProdSeq)):
test([1,2,3],[a,b]); %;
eval([seq(seq([i1, i2], i2 = [a, b]), i1 = [1, 2, 3])])
[[1, a], [1, b], [2, a], [2, b], [3, a], [3, b]]
```

The following tabulates, for each procedure, the time and memory required to generate the Cartesian product of a sequence of lists. Four different inputs are tested. For all cases the CartProdSeq is the fastest and uses the least memory.

```Compute the cartesian product of N lists of E elements
to get a list of E^N sublists of N elements

gcfreq = 10^7

+-------- timing -------+--------- memory ---------+
256^2                | maple (s) | total (s) | used (MiB) | alloc (MiB) | gctimes
---------------------+-----------+-----------+------------+-------------+--------
CartProdDeiterated   |     1.30  |     1.33  |    23.06   |     19.50   |    0
CartProdRtable       |     0.54  |     0.55  |    10.00   |      7.50   |    0
CartProdAsList       |     0.77  |     0.86  |    14.81   |     10.81   |    0
CartProdRecursive    |     0.25  |     0.26  |     7.79   |      5.75   |    0
CartProdNonrecursive |     0.13  |     0.13  |     5.77   |      3.75   |    0
CartProdSeq          |     0.12  |     0.13  |     5.76   |      3.75   |    0
---------------------+-----------+-----------+------------+-------------+--------
10^5                 |           |           |            |             |
---------------------+-----------+-----------+------------+-------------+--------
CartProdDeiterated   |     2.82  |     2.94  |    43.31   |     32.74   |    1
CartProdRtable       |     1.30  |     1.36  |    19.85   |     14.93   |    0
CartProdAsList       |     1.94  |     1.98  |    27.67   |     21.93   |    0
CartProdRecursive    |     0.74  |     0.74  |    18.41   |     13.56   |    0
CartProdNonrecursive |     0.47  |     0.49  |    14.78   |      9.94   |    0
CartProdSeq          |     0.32  |     0.35  |    13.53   |      9.37   |    0
---------------------+-----------+-----------+------------+-------------+--------
3^10                 |           |           |            |             |
---------------------+-----------+-----------+------------+-------------+--------
CartProdDeiterated   |     2.28  |     2.34  |    31.04   |     27.99   |    0
CartProdRtable       |     1.08  |     1.10  |    12.54   |     10.56   |    0
CartProdAsList       |     1.85  |     1.87  |    24.33   |     20.68   |    0
CartProdRecursive    |     0.66  |     0.73  |    18.49   |     13.37   |    0
CartProdNonrecursive |     0.43  |     0.46  |    14.54   |     10.19   |    0
CartProdSeq          |     0.24  |     0.27  |    10.56   |      8.25   |    0
---------------------+-----------+-----------+------------+-------------+--------
2^16                 |           |           |            |             |
---------------------+-----------+-----------+------------+-------------+--------
CartProdDeiterated   |     3.67  |     3.80  |    44.03   |     35.24   |    1
CartProdRtable       |     1.64  |     1.67  |    17.02   |     14.56   |    0
CartProdAsList       |     3.06  |     3.18  |    35.65   |     30.74   |    0
CartProdRecursive    |     1.15  |     1.31  |    28.41   |     23.31   |    0
CartProdNonrecursive |     0.78  |     0.83  |    23.88   |     19.25   |    0
CartProdSeq          |     0.44  |     0.48  |    18.21   |     14.00   |    0
---------------------------------------------------------------------------------
``` ﻿