hilbert2.clj

;; algorithm: see http://en.wikipedia.org/wiki/Hilbert_curve

(set! *warn-on-reflection* true)
(set! *unchecked-math* true)

;; TODO rewrite with *unchecked-math*

(defrecord Point [^long x ^long y])

(defn hilbert-to-point-unchecked [^long order ^long d]
  (let [n (bit-shift-left (unchecked-long 1) (unchecked-long order))]
    (loop [s 1, rx 0 ry 0, t (long d), x 0, y 0]
      (let [rx (bit-and 1 (bit-shift-right t 1))
            ry (bit-and 1 (bit-xor t rx))
            sm (unchecked-dec s)
            ;; TODO: reduce the number of ==1 checks
            x (if (== 1 ry) x (if (== 1 rx) (unchecked-subtract sm y) y))
            y (if (== 1 ry) y (if (== 1 rx) (unchecked-subtract sm x) x))
            x (unchecked-add x (unchecked-multiply s rx))
            y (unchecked-add y (unchecked-multiply s ry))
            t (bit-shift-right t 2)
            s (bit-shift-left s 1)]
        (if (< s n) (recur s rx ry t x y)
            (Point. x y))))))

(defn red [o]
  (let [m (bit-shift-left 1 (* 2 o))]
    (reduce (fn [a b] (Point. (unchecked-add (:x a) (:x b)) (unchecked-add (:y a) (:y b))))
            (Point. 0 0)
            (for [d (range 0 m)] (hilbert-to-point-unchecked o d)))))

(defn red-fast [o]
  (loop [^Point sum (Point. 0 0)
         r (for [d (range 0 (bit-shift-left 1 (* 2 o)))]
             (hilbert-to-point-unchecked o d))]
    (if r (recur (let [a ^Point (first r)] (Point. (unchecked-add (.x a) (.x sum))
                                                   (unchecked-add (.y a) (.y sum))))
                 (next r))
        sum)))

(defn red-faster [o]
  (loop [^Point sum (Point. 0 0)
         d 0]
    (let [a ^Point (hilbert-to-point-unchecked o d)]
      (if (< d (bit-shift-left 1 (* 2 o)))
        (recur (Point. (unchecked-add (.x a) (.x sum)) (unchecked-add (.y a) (.y sum))) (inc d) )
        sum))))

(prn (red-faster 10))

(doseq [x (range 0 100)]
  (time (red-faster 10)))