Permutations

permutron.core.clj

(ns permutron.core
  (:require [clojure.core.protocols :as ccp]
            [clojure.core.reducers :as r])
  (:import [clojure.lang Seqable Indexed Counted]))

;; Why are these private in clojure.core.reducers, but fjtask is public?
(def ^:private fjinvoke @#'r/fjinvoke)
(def ^:private fjfork @#'r/fjfork)
(def ^:private fjjoin @#'r/fjjoin)

(defn xrange
  ([] (xrange 0 Long/MAX_VALUE 1))
  ([end] (xrange 0 end 1))
  ([start end] (xrange start end 1))
  ([start end step]
     (reify
       Counted
       (count [_]
         (let [x (- end start), d (quot x step), r (mod x step)]
           (if (zero? r)
             d
             (inc d))))

       Indexed
       (nth [_ n]
         (let [r (+ start (* step n))]
           (if (< r end)
             r
             (throw (IndexOutOfBoundsException.)))))

       Seqable
       (seq [_]
         ((fn stepf [x]
            (lazy-seq
             (if (< x end)
               (cons x (stepf (+ x step))))))
          start))

       ccp/CollReduce
       (coll-reduce [this f] (ccp/coll-reduce this f (f)))
       (coll-reduce [_ f init]
         (loop [acc init, i start]
           (cond
            (reduced? acc) @acc
            (>= i end) acc
            :else (recur (f acc i) (+ i step)))))

       r/CollFold
       (coll-fold [r n combinef reducef]
         (let [nr (count r)]
           (cond
            (zero? nr) (combinef)
            (<= nr n) (reduce reducef (combinef) r)
            :else
            (let [split (-> nr (quot 2) (* step) (+ start))
                  r1 (xrange start split step)
                  r2 (xrange split end step)
                  fc (fn [r] #(r/coll-fold r n combinef reducef))]
              (fjinvoke
               #(let [f1 (fc r1)
                      t2 (r/fjtask (fc r2))]
                  (fjfork t2)
                  (combinef (f1) (fjjoin t2)))))))))))

(defn xrange-onto
  [mapf n]
  (reify
    Counted (count [_] n)
    Seqable (seq [_] (seq (map mapf (xrange n))))

    Indexed
    (nth [_ i]
      (if (< i n)
        (mapf i)
        (throw (IndexOutOfBoundsException.))))

    ccp/CollReduce
    (coll-reduce [this f] (ccp/coll-reduce this f (f)))
    (coll-reduce [_ f init]
      (r/reduce f (r/map mapf (xrange n))))

    r/CollFold
    (coll-fold [_ pn combinef reducef]
      (r/fold pn combinef reducef (r/map mapf (xrange n))))))

(def ^:const ^long fac-max-precalc
  "Maximum value of `fac` to pre-calculate."
  20)

(let [fac* (fn [^long n] (reduce *' 1 (xrange 1 (inc n))))
      memo (long-array (inc fac-max-precalc))
      nmemo (long (alength memo))]
  (aset-long memo 0 1)
  (doseq [i (xrange 1 nmemo)]
    (aset-long memo i (* i (aget memo (dec i)))))

  (defn fac
    "Factorial of `n`."
    ^long [^long n]
    (if (> nmemo n)
      (aget memo n)
      (throw (ex-info "Result too large for JVM long." {:n n}))))

  (defn fac'
    "Factorial of `n`."
    [^long n]
    (if (> nmemo n)
      (aget memo n)
      (fac* n))))

(defn permutation*
  [^long n ^long x]
  (let [p (long-array n)
        _ (loop [x (long x), f (long 2), i (- n 2)]
            (aset-long p i (mod x f))
            (if (pos? i)
              (recur (quot x f) (inc f) (dec i))))
        _ (loop [i 1, m (bit-shift-left 1 (aget p 0))]
            (let [x (loop [x (aget p i), x' 0]
                      (if-not (zero? (bit-and m (bit-shift-left 1 x')))
                        (recur x (inc x'))
                        (if (zero? x)
                          x'
                          (recur (dec x) (inc x')))))
                  _ (aset-long p i x)
                  i (inc i)]
              (if (< i n)
                (recur i (bit-or m (bit-shift-left 1 x))))))]
    p))

(defn permutation
  [xs ^long x]
  (let [n (count xs), ^longs p (permutation* n x)]
    (mapv (partial nth xs) p)))

(defn permutations
  [xs]
  (let [f (partial permutation xs)
        n (fac (count xs))]
    (xrange-onto f n)))

(defn cartprod*
  [xs i]
  (first
   (reduce (fn [[cp i j] f]
             (let [x (nth xs j)]
               [(conj cp (nth x (mod i f))) (quot i f) (inc j)]))
           [[] i 0] (r/map count xs))))

(defn cartprod
  [xs]
  (let [f (partial cartprod* xs)
        n (r/reduce * (r/map count xs))]
    (xrange-onto f n)))