Re: challenge: set partitions and power set

[adorf@eso.org (Hans-Martin Adorf)]

Hello MCLers,

as an addendum to my previous summary post, here is a small function for
computing the number of set partitions.

Enjoy.

Hans-Martin Adorf

;;;;----------------------------------------------------------------------------
;;;; Bell's number, the number of partitions of an n-element set
;;; 
;;;; Hans-Martin Adorf
;;;; ST-ECF/ESO
;;;; Karl-Schwarzschild-Str. 2
;;;; D-85748 Garching b. Muenchen
;;;; adorf@eso.org
;;;;
;;;; upon a suggestion by Mark McConnell <mmcconn@math.okstate.edu>
;;;;----------------------------------------------------------------------------

(defmacro sum ((i i-min i-max &optional (i-step 1)) expr)
  "Sum over indexed expression from i = i-min to i = i-max (inclusive)"
  `(do ((,i ,i-min (incf ,i ,i-step))
        (result 0))
       ((> ,i ,i-max) result)
     (incf result ,expr)))
;; (sum (j 1 3) (* j j))

(defmacro product ((i i-min i-max &optional (i-step 1)) expr)
"Compute the product over indexed expression from i = i-min to i = i-max
(inclusive)"
  `(do ((,i ,i-min (incf ,i ,i-step))
        (result 1))
       ((> ,i ,i-max) result)
     (setf result (* result ,expr))))
;; (product (j 1 3) (* j j))

(defun bell (n)
  "Bell's number B(n), the number of partitions of an n-element set"
  (if (= n 0)
    1
    (sum (k 0 (1- n)) 
         (* (binomial k (1- n)) (bell k)))))
#|
(bell 0)
(bell 1)
(bell 2)
(bell 3)
(bell 4)
(bell 5)
(bell 6)
|#

;; elegant implementation using rationals
(defun binomial (k n)
  "Binomial coefficient"
  (product (i 0 (1- k))
           (/ (- n i) (1+ i))))
;; (binomial 2 3)

#|
;; less elegant implementation using rationals
(defun binomial (k n &aux (result 1))
  "Binomial coefficient"
  (dotimes (i k result)
    (setf result (* result (/ (- n i) (1+ i))))))

;; leads to very large integer numbers
(defun binomial (k n)
  "Binomial coefficient"
  (when (<= k n)
    (/ (factorial n) 
       (factorial k) (factorial (- n k)))))
;; (binomial 0 2)

(defun factorial (n)
  (if (< n 1)
    1
    (* n (factorial (1- n)))))
;; (factorial 10)
|#