Cartesian Product

September 6, 2013

We start with the recursive version of the function. I think it looks prettiest in Haskell:

cprod [ ] = [[ ]]
cprod (xs : zss) = [ x : ys | x <- xs, ys <- cprod zss ]

Since our Standard Prelude includes list comprehensions, that’s easy to translate to Scheme. It gets a little bit uglier because Scheme assembles the arguments to a function into a list when called with define‘s dot-notation, which is the only way to pass a variable number of arguments, and the recursive call has to pick apart the pieces of the list with apply so the define can package them up again:

(define (cprod . xss)
  (if (null? xss) (list (list))
    (list-of (cons x ys)
      (x in (car xss))
      (ys in (apply cprod (cdr xss))))))

> (cprod '(1 2 3) '(4) '(5 6))
  ((1 4 5) (1 4 6) (2 4 5) (2 4 6) (3 4 5) (3 4 6))

If you object to the use of the list comprehension, here is a version that has the two inner loops written out explicitly:

(define (cprod . xss)
  (if (null? xss) (list (list))
    (let loop1 ((xs (car xss)) (zs (list)))
      (if (null? xs) (reverse zs)
        (let loop2 ((yss (apply cprod (cdr xss))) (zs zs))
          (if (null? yss) (loop1 (cdr xs) zs)
            (loop2 (cdr yss) (cons (cons (car xs) (car yss)) zs))))))))

> (cprod '(1 2 3) '(4) '(5 6))
  ((1 4 5) (1 4 6) (2 4 5) (2 4 6) (3 4 5) (3 4 6))

We turn now to the iterative version. The basic idea is similar to the odometer in your car, except that instead of counting from 0 to 9 the digits count from 0 to the length of the list; at each step the odometer digits give the index of the item in each list that appears in the current cross-product. I find it best to think about this as counting in a mixed-radix positional number system; one function converts a number to a list of mixed-radix digits, a second function counts up from zero and outputs the list elements corresponding to the mixed-radix digits:

(define (mixed-radix n bases)
  (let loop ((n n) (bases bases) (xs (list)))
    (cond ((zero? n)
            (let loop ((k (length bases)) (xs xs))
              (if (zero? k) xs (loop (- k 1) (cons 0 xs)))))
          ((null? bases) (list))
            (else (loop (quotient n (car bases)) (cdr bases)
                        (cons (remainder n (car bases)) xs))))))

(define (cprod . xss)
  (let ((bases (reverse (map length xss))))
    (let loop ((n 0) (zss (list)))
      (let ((odometer (mixed-radix n bases)))
        (if (null? odometer) (reverse zss)
          (let ((zs (map list-ref xss odometer)))
            (loop (+ n 1) (cons zs zss))))))))

> (cprod '(1 2 3) '(4) '(5 6))
  ((1 4 5) (1 4 6) (2 4 5) (2 4 6) (3 4 5) (3 4 6))

The cartesian product function in the Standard Prelude, called cross, uses a functional pearl collected by Mike Spivey that uses the higher-order function fold-right so that neither recursion nor iteration is apparent in the source code:

(define (cross . xss)
  (define (f xs yss)
    (define (g x zss)
      (define (h ys uss)
        (cons (cons x ys) uss))
      (fold-right h zss yss))
    (fold-right g '() xs))
  (fold-right f (list '()) xss))

> (cross '(1 2 3) '(4) '(5 6))
  ((1 4 5) (1 4 6) (2 4 5) (2 4 6) (3 4 5) (3 4 6))

You can run the program at http://programmingpraxis.codepad.org/WanFTb6R.

About these ads

Pages: 1 2

8 Responses to “Cartesian Product”

  1. Josef Svenningsson said

    A recursive Haskell version:


    cart [] = [[]]
    cart (xs:xss) = [ x:ys | x <- xs, ys <- cart xss ]

  2. Jussi Piitulainen said

    Simple-minded recursion:

    (define (X . sets)
      (if (null? sets) '(())
          (let ((tails (apply X (cdr sets))))
            (apply append
                   (map (lambda (h)
                          (map (lambda (t) (cons h t)) tails))
                        (car sets))))))

    ;(X '(1 2 3) '(4) '(5 6))
    ;==> ((1 4 5) (1 4 6) (2 4 5) (2 4 6) (3 4 5) (3 4 6))

    Differently iterative:

    (define (for-X proc . sets)
      (if (null? sets) (proc '())
          (for-each
           (lambda (head)
             (apply for-X (lambda (tail)
                            (proc (cons head tail)))
                    (cdr sets)))
           (car sets))))

    ;(for-X write '(1 2 3) '(4) '(5 6)) ==> unspecified
    ;displays: (1 4 5)(1 4 6)(2 4 5)(2 4 6)(3 4 5)(3 4 6)

  3. Paul said

    A Python version. In Python this functionality is already in library itertools (product).

    def product_iteration(*lists):
        """a simplified copy of the itertools product"""
        result = [[]]
        for alist in map(tuple, lists):
            result = [x + [y] for x in result for y in alist]
        for prod in result:
            yield tuple(prod)
            
    def product_recursion(*lists):
        if not lists:
            yield ()
        else:
            for alist in lists[0]:
                for prod in product_recursion(*lists[1:]):
                    yield (alist,) + prod
                    
    def cartesian(*lists, **kwds):
        """default: iteration method
           use method=product_recursion for recursion method
        """
        method = kwds["method"] if "method" in kwds else product_iteration
        return tuple(method(*lists))
    
    print cartesian([1,2,3], [4],[5,6])
    print cartesian([1,2,3], [4],[5,6], method=product_recursion)
    #->((1, 4, 5), (1, 4, 6), (2, 4, 5), (2, 4, 6), (3, 4, 5), (3, 4, 6))
    
  4. #lang racket
    
    (define/match (cartesian-product/r . xss)
      [('()) '(())]
      [((cons xs yss))
       (for*/list ([ys (apply cartesian-product/r yss)]
                   [x xs])
          (cons x ys))])
    
    (define (cartesian-product/i . xss)
      (map reverse
           (for/fold ([ps '(())]) ([xs xss])
             (for*/list ([p ps] [x xs])
               (cons x p)))))
    
  5. Peter Salvi said

    Another recursive implementation in Common Lisp:

    (defun cartesian-rec (&rest lists)
      (if (null lists)
          '(())
          (mappend (lambda (x)
                     (mapcar (lambda (lst) (cons x lst))
                             (apply #'cartesian-rec (cdr lists))))
                   (car lists))))
    

    And, for fun, here is one with macros, generating a separate loop for every level:

    (defun cartesian-loops (lst vars)
      (if (null lst)
          `(push (list ,@(reverse vars)) acc)
          (let ((v (gensym)))
            `(dolist (,v ,(car lst))
               ,(cartesian-loops (cdr lst) (cons v vars))))))
    
    (defmacro cartesian (&rest lists)
      `(let ((acc '()))
         ,(cartesian-loops lists '())
         (nreverse acc)))
    
  6. brooknovak said

    Here is an iterative and recursive C# solution:

    public interface ICartesianProduct
    {
    	int[][] GetProduct(int[][] sets);
    }
    
    public class CartesianProductIterative : ICartesianProduct
    {
    	public int[][] GetProduct(int[][] inputSets) {
    		var resultSetCount = 1;
    		foreach (var s in inputSets) {
    			resultSetCount *= s.Length;
    		}
    		var resultSets = new int[resultSetCount][];
    		var resultCounter = 0;
    
    		// Create initial set of cartesian products built from the last set from the input
    		for (int j = 0; j < inputSets[inputSets.Length - 1].Length; j++) {
    			resultSets [j] = new int[inputSets.Length];
    			resultSets [j] [inputSets.Length - 1] = inputSets [inputSets.Length - 1] [j];
    			resultCounter++;
    		}
    		// Build up all remaining sets in the cartesian product
    		for (int i = inputSets.Length - 2; i >= 0; i--) {
    			// Inject first item in set [i], into previous built cartesian product sets
    			for (int j = 0; j < resultCounter; j++) {
    				resultSets [j] [i] = inputSets [i] [0];
    			}
    
    			// create new copies of all cartesian sets previously built (up to sets [i-1]) - but swap out the
    			// items in the current set [i] so all items in set [i] have been enumerated
    			for (int j = 1; j < inputSets[i].Length; j++) {
    				for (int k = 0; k < resultCounter; k++) {
    					int newIdx = (resultCounter * j) + k;
    					resultSets [newIdx] = new int[inputSets.Length];
    					Array.Copy (resultSets [k], i + 1, resultSets [newIdx], i + 1, inputSets.Length - (i + 1));
    					resultSets [newIdx] [i] = inputSets [i] [j];
    				}
    			}
    
    			resultCounter *= inputSets[i].Length;
    		}
    
    		return resultSets;
    	}
    }
    
    public class CartesianProductRecursive : ICartesianProduct
    {
    	public int[][] GetProduct(int[][] sets) {
    		var resultSets = new List<int[]> ();
    		var indices = new int[sets.Length];
    		CreateProductSet (sets, indices, 0, resultSets);
    		return resultSets.ToArray();
    	}
    
    	private void CreateProductSet(int[][] inputSets, int[] indices, int currentSet, List<int[]> resultSets) {
    		for (int i = 0; i < inputSets[currentSet].Length; i++) {
    			indices [currentSet] = i;
    			if (currentSet == (inputSets.Length - 1)) {
    				var cartesianSet = new int[inputSets.Length];
    				for (int j = 0; j < inputSets.Length; j++) {
    					cartesianSet[j] = (inputSets[j][indices[j]]);
    				}
    				resultSets.Add (cartesianSet);
    			} else {
    				CreateProductSet (inputSets, indices, currentSet + 1, resultSets);
    			}
    		}
    	}
    }
    
  7. treeowl said

    Haskell doesn’t really have functions with variable numbers of arguments, and faking such is more for programmers of Okasaki’s caliber than mine, so I’ll just have it take a list of lists. I started out with a recursive solution, then realized I could make it a little prettier:

    cart2 xs ys = [x:y | x <- xs, y  [[a]]
    cart ll = foldr cart2 [[]] ll
    

    Then I looked on the Internet for advice. Oddly enough, it turns out that Haskell has Cartesian products in the Prelude, but not in an obvious way.

    cart :: [[a]] -> [[a]]
    cart = sequence
    

    does the trick. It took me a while to figure out how it works, so I figured I’d type up my reasoning below:

    sequence is actually a rather more general. Here’s how it’s defined at http://hackage.haskell.org/package/base-4.6.0.1/docs/src/Control-Monad.html#sequence

    -- | Evaluate each action in the sequence from left to right,
    -- and collect the results.
    sequence       :: Monad m => [m a] -> m [a] 
    {-# INLINE sequence #-}
    sequence ms = foldr k (return []) ms
                where
                  k m m' = do { x <- m; xs <- m'; return (x:xs) }
    

    The funky return [] can immediately be rewritten based on the instance declaration for lists in http://hackage.haskell.org/package/base-4.6.0.1/docs/src/GHC-Base.html#Monad partially copied below:

        m >>= k            = foldr ((++) . k) [] m
        return x            = [x]
    

    In particular, note that return [] is actually the same as [[]]. Also, for lists, we can translate the do notation directly to list comprehension notation, so the above definition for sequence, specialized to lists, is written much more transparently thus:

    sequence :: [[a]] -> [[a]]
    sequence ls = foldr k [] ms
       where
          k m m' = [x:xs | x <- m, xs <- m']
    

    Now it’s much clearer what is going on: this is just translating the Cartesian product of [l1,l2,...,ln] into l1*(l2*(l3*…)).
    Now how is that do notation or list comprehension treated under the hood? Well, it’s translated to

    m >>= \x -> m' >>= \xs -> return (x:xs)
    

    Now what does that funky >>= do? Well, m >>= k is foldr ((++).k) [] m, which is not immediately helpful. But letting m=[m1,m2,...,mn] and expanding the foldr simplifies things a drop:

    m1 `(++).k` (m2 `(++).k` (m3 `(++).k` (... [])))
    

    It’s still a tad weird, though. What the heck is (++).k? It’s a function, obviously, that takes a value, applies k to it, and then applies (++) to the result. Writing this out, it’s

    \x -> (++) k(x)
    

    or, equivalently,

    \x y -> k(x) ++ y
    

    That is, p `(++).k` q = k(p) ++ q.
    Substituting this in, m >>= k becomes much simpler:

    k(m1) ++ k(m2) ++ ... ++ []

    so we’re just applying k to each of the elements of m and concatenating the results!
    So

    m >>= \x -> m' >>= \xs -> return (x:xs)
    

    just means that for each element x of m, for each element xs of m’ we form the list [x:xs], concatenate these together to form [x:xs1++x:xs2++...], and concatenate those together to form [x1:xs1++x1:xs2++...++x2:xs1++x2:xs2++...], which is the Cartesian product.

  8. treeowl said

    Which is the Cartesian product of x with xs, I should say. Now that I understand what the helper function “k” in sequence does, the operation of sequence becomes clear.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 598 other followers

%d bloggers like this: