We start with the definition of the problem: given a set 'l' and a
number 'n', return the set of all subsets of 'l' of cardinality
'n'. Sets are represented by lists.

The key was to choose the right definition.

Let ps(L) is a powerset of L that does not include {empty set} (that
is, a singular set whose element is the empty set).
        length(ps(L)) = (- (expt 2 (length L)) 1)

Let L be a non-empty set and let L = A U B where A and B are two
disjoint subsets of L.

Obviously,
        ps(L) = ps(A U B) =
        ps(A) U ps(B) U { y U x | x <- ps(A), y <- ps(B) }
Let (subsets L n) = (filter (lambda (el) (= n (length el))) ps(L) )

Thus the desired function 'subsets' is a filtered powerset ps. This
seems to be a stupid definition. Let's not be hasty however. Note that
filter and union commute: the filter of a union of sets is the union
of filtered sets. Therefore, from the previous expression

        (subsets L n) = (subsets (union A B) n)
        = (subsets A n) U (subsets B n)
        Union{ y U x | x <- (subsets A k), y <- (subsets B (- n k)),
               k=1,n-1 }

Well, this is it. Note, we didn't say how to split L into two disjoint
subsets A and B. We can do as we wish. For example, we can choose to
split in such a way so that (length B) is n. In the most difficult
case where n = (/ (length L) 2), this corresponds to a "divide and
conquer" strategy, so to speak (at least in the first stages).

The Scheme code below implements this idea verbatim. It uses the
accumulator-passing style that was expounded earlier.

(define (subsets-v5 l n)

  ; The initialization function. Check the boundary conditions
  (define (loop l ln n accum)
    (cond
     ((<= n 0) (cons '() accum))
     ((< ln n) accum)
     ((= ln n) (cons l accum))
     ((= n 1)
      (let fold ((l l) (accum accum))
        (if (null? l) accum
            (fold (cdr l) (cons (cons (car l) '()) accum)))))
     (else
      (split l ln n accum))))

  ; split l in two parts a and b so that (length b) is n
  ; Invariant: (equal? (append a b) l)
  ; ln is the length of l
  (define (split l ln n accum)
    (let loop  ((a '()) (b l) (lna 0) (lnb ln))
      (if (= lnb n) (cont a lna b lnb n accum)
          (loop (cons (car b) a) (cdr b) (+ 1 lna) (- lnb 1)))))

  ; The main body of the algorithm
  (define (cont a lna b lnb n accum)
    (let* ((accum
            (loop a lna n accum))
           (accum       ; this is actually (loop b lnb n accum)
            (cons b accum))
           )
      (let inner ((k 1) (accum accum))
        (if (> k (min lna (- n 1))) ; don't loop past meaningful boundaries
            accum
            (let ((as (loop a lna k '()))
                  (bs (loop b lnb (- n k) '())))
              (inner (+ 1 k)
                     ; compute the cross-product of as and bs
                     ; and append it to accum
                     (let fold ((bs bs) (accum accum))
                       (if (null? bs) accum
                           (fold (cdr bs)
                                 (append
                                  (map (lambda (lst) (append lst (car bs))) as)
                                  accum))))))))))

  (loop l (length l) n '()))

The benchmark runs in 993 ms of user time and allocates only 36.5 MB
of memory, on Gambit-C interpreter. This is the absolute, incredible
record. Under SCM:

subsets-v3 (called combos by John David Stone)
;Evaluation took 1596 mSec (98 in gc) 657662 cells work, 4721364 env, 97 bytes other

subsets-v5:
;Evaluation took 700 mSec (322 in gc) 1708112 cells work, 610264 env, 105 bytes
other                                                                          
That is, more than twice as fast.

Continuing the table from the previous post:

Procedure       Gambit-C               Bigloo 2.4b      Bigloo 2.4b
                interpreter, s         interpreter, s   compiler, s
subsets-v0        285.0                   11.59          5.62    3.14
subsets-v1          6.3                    5.45          2.22    0.34
subsets-v3          8.1                    4.78          0.96    0.27
subsets-v20        14.1                    5.53          0.96    0.26    
subsets-v21         7.7                    4.88          0.66    0.26
subsets-v22         5.0                    3.18          0.62    0.25
subsets-v23         4.1                    2.86          0.82    0.25
subsets-v5          0.9                    1.56          1.10    0.76

Well, compiled code isn't that fast -- but that because subsets-v5
isn't too optimal. The append and map ought to be converted into the
accumulation-passing style. That will reduce the amount of garbage as
well. But it's past midnight.

    But what _is_ subsets function?  What does it
compute, and what is it used for?

    I would be grateful for additional comments.

        That is:  Apart from garbage collection, the two procedures are
very closely comparable -- but SUBSETS-V3 creates only a third as much
garbage.

        My guess is that, in my testing environment, SUBSETS-V3 is getting
a bigger advantage from Chez Scheme's optimizations.  The results I get under
MzScheme are similar to Oleg's:

And he's right.

  http://sardine.cs.columbia.edu:8080/subsets/

All graphs are timings of different lengths out of a set of 22
elements, msubset is the memoization of Oleg's first version, msubset1
is even simpler, removing the special case for len=1.  The code is at

  http://sardine.cs.columbia.edu:8080/subsets/subsets.ss

but this is just for reference, it contains all the code pieces with
some junk that makes gnuplot doo all the graphs.

Obviously a lazy functional language like Haskell shows this
implementation in its best light.  As Oleg posted, we can define the
subsets function mathematically as

subsets(S, 0)       = {{}}
subsets({}, n)      = {}
subsets({a} U S, n) = {{a}} x subsets(S, n-1)  U  subsets(S, n)

This translates directly into Scheme or Haskell.  Scheme code has
already been posted.  In Haskell it looks like this:

subsets _      0 = [[]]
subsets []     _ = []
subsets (x:xs) n = map (x :) (subsets xs (n - 1)) ++ (subsets xs n)

This is slow, partly because it conses too much but mostly because the
same subproblems are solved repeatedly.  Memoization can be used to
attack this performance problem.  The lazy streams approach will give
many the benefits of memoization while avoiding most of the
disadvantages.

Memoization pros

 - can be applied automatically without changing the naive
   implementation

 - memoized results can be saved across top level invocations of the
   function

Memoization cons

 - not always faster

 - doesn't work well if function arguments are complex (the required
   hash compare might be too slow)

 - usually implemented using imperative machinery under the hood
   (although the resulting memoized function has functional behavior)

 - can lead to horrendous space leaks in the memo tables; recovering
   this space might require manual control partly negating the first
   pro listed above

In a lazy stream formulation, the nth element of the lazy stream will
be a list of the subsets of cardinality n.  In other words, the stream
will be [subsets(S,0), subsets(S,1), subsets(S,2), ...].  The first
stream element will always be the list of the empty set.  For all n
greater than length of the input list the element will be null.  This
is very smooth in Haskell:

subsets s n = (subsets_stream s) !! n

subsets_stream []     = [[]] : repeat []
subsets_stream (x:xs) =
    let r = subsets_stream xs
        s = map (map (x:)) r
    in [[]] : zipWith (++) s (tail r)

As a small bonus We can easily use subsets_stream to implement the
powerset function:

powerset s = concat (takeWhile (not . null) (subsets s))

Unfortunately this isn't as pretty in Scheme.  Comparing the Scheme
implementation below to Haskell, Scheme is far more verbose since we
have to explicitly request laziness, we have to provide lazy
equivalents of standard library functions (CONS, CAR, CDR, LIST-REF
and MAP, although these could be put in an external library) and
finally currying is more convenient to use in Haskell than it is in
Scheme.  (Often Haskell's pattern matching would be an advantage as
well, but here it plays no important part.)  One tricky part about
explicitly requesting laziness is that it's easy to get wrong.  Look
at the mix of uses of STREAM-MAP vs MAP and STREAM-CONS vs CONS in
SUBSET-STREAMS below.

Since most Scheme implementations don't natively support R5RS macros,
I have simply replaced all uses of STREAM-CONS with
(CONS (DELAY a) (DELAY b)).  Notice that I use fully lazy streams.
The streams in SICP have lazy tails but eager heads.

; Most Scheme implementations don't have R5RS macros, alas.
;(define-syntax stream-cons
;  (syntax-rules ()
;    ((stream-cons hd tl) (cons (delay hd) (delay tl)))))

(define (subsets s n)
  (stream-ref (subsets-stream s) n))

;; Return a lazy stream of all subsets of list XS.  The nth element of
;; the stream is a list of the subsets of cardinality n (counting from 0).
(define (subsets-stream xs)
  (if (null? xs)
      (cons (delay (list '())) (delay (stream-repeat '())))
      (let* ((r (subsets-stream (cdr xs)))
             (s (stream-map (lambda (ss)
                              (map (lambda (y) (cons (car xs) y)) ss))
                            r)))
        (cons (delay (list '()))
              (delay (stream-map append s (stream-cdr r)))))))

;; The remainder should be part of a lazy streams library.
;; These are the lazy list equivalents of CAR, CDR, LIST-REF and MAP.
;; STREAM-REPEAT is inspired by Haskell's repeat function.

(define (stream-car s)
  (force (car s)))

(define (stream-cdr s)
  (force (cdr s)))

;; Create an infinite stream of X's.
(define (stream-repeat x)
  (cons (delay x) (delay (stream-repeat x))))

(define (stream-ref s n)
  (if (= n 0)
      (stream-car s)
      (stream-ref (stream-cdr s) (- n 1))))

;; Answer to an easy exercise in SICP.  Most SICP exercises are harder ;-)
(define (stream-map f . streams)
  (if (null? (car streams))
      '()
      (cons (delay (apply f (map stream-car streams)))
            (delay (apply stream-map
                          (cons f (map stream-cdr streams)))))))

  http://sardine.cs.columbia.edu:8080/subsets/

Saving results across top level invocations is at odds with the space
leak problem I mention later.  For a restricted set of arguments for
this specific problem, you can encapsulate the answers for all of the
integer subsets in one place.  This is a bit ugly in Scheme, and it
also changes the arguments to the subsets function.  (Instead of
passing a list and a number we need two numbers.)

(define (ints-from n) (cons (delay n) (delay (ints-from (+ n 1)))))

(define (int-subsets-stream-from n)
  (cons (delay (subset-stream (ints-from 1))
        (delay (int-subsets-stream-from (+ n 1))))))

;; Return a list of all subsets of the integers {1, ..., N} of
;; cardinality R.
(define subsets
  (let ((all-subsets (cons (delay '(())
                           (delay (int-subsets-stream-from 1))))))
    (lambda (n r)
      (stream-ref (stream-ref all-subsets n) r))))

Space leaks are a problem here and would require manual intervention.
(No better than memoization on this point.)

Ignoring the space issues, this is a one-liner in Haskell:

int_subsets = [subsets_stream [1 .. n] | n <- [0 ..]]
subsets n r = (int_subsets !! n) !! r

I didn't mention the performance benefits of memoization from avoiding
repeated computation of the same subproblems since it is an obvious
feature.  I was thinking of it when I wrote "many", but perhaps many
wasn't the right word when I meant just two out of three.

Of course laziness can certainly hurt performance as well.  In Haskell
this problem shows up in the other direction -- you have to decide
when you want to be strict instead of lazy.  My guess is that this is
normally harder than deciding when memoization will help.

I agree completely on the advantage of memoization over dynamic
programming.  The memoized function treated as a black box is a
functional unit, and standard dynamic programming techniques can't be
treated that way.  The dynamic programming techniques are normally a
lot harder to code correctly too, so bugs are much more likely.

Does anyone have a good dynamic programming implementation of the
classic count_change problem?  I looked at this, but it's hard to get
the dynamic programming implementation correct.  This problem memoizes
very nicely, and there is an elegant lazy streams implementation.
"Concrete Mathematics" gives the solution in closed form, but the
closed form looks to be about as hard (or harder!) to compute than
just using efficient ways to enumerate the combinations.

Weak pointers are an excellent solution with one caveat: they aren't a
standard part of Scheme so you can't write your code in R5RS.  Weak
pointers need gc support so they aren't available in all Scheme
implementations .  (But I think they should be.  Hash consing is a
specific application of memoization that is very useful but really
requires weak pointers to be practical.)  Of course lazy evaluation is
subject to space leaks also, and I think those leaks are usually
harder to track down, understand and fix than memo table leaks.

I like memoization and I didn't intend to try to bury it.  (I should
have simply provided the lazy streams code and avoided
editorializing!)  I only wanted to compare memoization to laziness
since both techniques are often useful in the same problem domains.
When both techniques apply to a problem, I think memoization is
generally easier and more natural to use in strict languages like
Scheme and SML.  Lazy streams are shown to better advantage in lazy
languages like Haskell.  (Of course I first saw lazy streams in Peter
Henderson's book which uses them in a lisp environment.  At that time
Haskell didn't yet exist and I was unfamiliar with KRC and Miranda
etc.  Actually I must have seen Burge's book "Recursive Programming
Techniques" around that time, but I didn't understand it!)

subsets [1 .. 1000] 10

in an interactive Haskell system like Hugs and it will start printing
solutions immediately.  I don't think the computation would complete
before the universe ends (there are 1000!/((1000 - 10)! x 10!)
different subsets in the result list), but you can immediately start
processing as many of them as you time for.  (The Scheme
implementation could be made to do this too.  It would require using a
special writer to display the resulting lazy list, forcing each
element just before printing it.  Fully lazy languages like Haskell do
this automatically.)

[I'll continue very briefly due to time constraints and avoiding "me
too"s.]

| * Some implementations may implement "implicit forcing," where the
|   value of a promise is forced by primitive procedures like cdr and +:
|     (+ (delay (* 3 7)) 13)  ==>  34

I wonder if it is possible to hack something up to make it do so
without a full evaluator implementation...

bytes.  Interpreted, the computation takes 38298824 bytes,
compiled, it takes 31401552 bytes using subsets-v5. That's
not so bad.

(define (count! l)
  (define c 0)
  (define !! (list '!!))
  (let loop ((l l))
    (if (and (pair? l) (not (eq? (car l) !!)))
      (let ((x (car l)) (y (cdr l)))
        (set-car! l !!)
        (set! c (+ 1 c))
        (loop x) (loop y))
      0))
  c)

and then:

| > (define x (subsets-v5 '(a b c d e f g h i j k l m n o p q r s t) 10))
| > (length x)
| 184756
| > (count! x)
| 1111280

That is, only ~1MB was actually used.

  537471 = [ C(20,10) + C(19,9) + C(18,8) + ... + C(11,1) ]  +  C(20,10)

where C(n,m) = n!/(m!n-m!).  (The extra C(20,10) are the pairs
connecting the subsets into a list.)  I'll leave the proof of this
claim as an exercise for Oleg.

537471 is the number of conses done by the accumulator-based version:

(define (subsets l n)
  (let ss ((l l) (len (length l)) (n n) (little-acc '()) (big-acc '()))
    (cond ((zero? n) (cons little-acc big-acc))
          ((> n len) big-acc)
          (else (ss (cdr l) (- len 1) (- n 1) (cons (car l) little-acc)
                    (ss (cdr l) (- len 1) n little-acc big-acc))))))

You can save an additional whopping ten pairs by reusing parts of the
argument.  I proclaim the following the STINGIEST subsets function
ever known to man:

(define (subsets l n)
  (let ss ((len (length l)) (n n) (little-acc '()) (big-acc '()))
    (cond ((zero? n) (cons little-acc big-acc))
          ((> n len) big-acc)
          (else (let ((last-pair (list-tail l (- len 1))))
                  (ss (- len 1) (- n 1)
                      (if (eq? (cdr last-pair) little-acc)
                          last-pair
                          (cons (car last-pair) little-acc))
                      (ss (- len 1) n little-acc big-acc)))))))

-al

  subsets      537471
  subs         537471
  subs2        537471
  combos       537471
  subset      2032316
  msubset      537526
  msubset1     537471
  subsets-v23  537471
  subsets-v5  1111280
  ssubsets     537471

An elegant proof still eludes me.

Bob