hilbert2.clj

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

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

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

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

(defn red [o]
  (let [m (unchecked-int (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 (unchecked-int (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 (hilbert-to-point-unchecked o d)]
      (if (< d (unchecked-int (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 30)]
  (time (red-faster 10)))