04 June 2010

Loop-Unrolling Macros in Clojure

I've been fooling around with clojure a bit lately. Something fun: a dotimes-like macro, but it completely unrolls your loop at compile time:
(defmacro dotimes-unroll
  "Just like dotimes, but completely unrolls the given loop.  The
number of iterations must read as a number."
  [[i n] & body]
  (assert (number? n))
  (if (= n 0)
    `nil
    `(let [~i ~(- n 1)]
       (dotimes-unroll [~i ~(- n 1)] ~@body)
       ~@body)))
Or, unroll by a given number of repetitions:
(defmacro dotimes-partial-unroll
  "Like dotimes, but partially unrolls the loop.  The number of
  repetitions attempted is the first argument to the macro, which must
  read as a literal number."
  [nstep [i nexpr] & body]
  (assert (number? nstep))
  `(let [n# (int ~nexpr)
         nmain# (quot n# ~nstep)
         nextra# (rem n# ~nstep)
         extrainc# (* nmain# ~nstep)]
     (loop [iouter# (int 0)]
       (if (>= iouter# nmain#)
         (dotimes [j# nextra#]
           (let [~i (int (+ j# extrainc#))]
             ~@body))
         (let [inc# (int (* iouter# ~nstep))]
           (dotimes-unroll [j# ~nstep]
             (let [~i (int (+ j# inc#))]
               ~@body))
           (recur (int (+ iouter# 1))))))))
I have no idea whether (or to what extent) this sort of loop unrolling is done by the JVM. So, these macros may be completely useless and redundant, but they were fun to write.

04 March 2010

Distribution for the Ratio of Two Uniform Deviates

Weird: the distribution for the ratio of two uniform deviates is constant for ratios smaller than 1, and falls like r^2 for larger ratios. p(r) = 1/2 if r < 1, and 1/(2r^2) if r >= 1 when r = a/b, where a,b ~ U(0, A). What's up with that?

09 January 2010

Another Library: 1-D Root Finding in Ocaml

...and here's another library for 1-D root finding (non-linear equation solving) in OCaml. Get it with
git clone git://github.com/farr/ocaml-solve.git
I'm particularly proud of this one for two reasons:
  1. All the algorithms, even the ones requiring derivatives, use bracketing. If a step in the algorithm would fall outside the brackets for the root, the algorithms perform a bisection step instead, and then try again. In this way, all algorithms are guaranteed to converge on a root. (It's always annoyed me that the GSL doesn't do this with its "algorithms that require derivatives".)
  2. One of the algorithms in the library is a higher-order Newton's method. Given an arbitrary number of derivatives, this algorithm takes advantage of the additional information to implement higher-order convergence. I learned this trick from Danby's book, where he is trying to solve Kepler's equation using a fourth-order method.

Polynomials in OCaml

I just posted some code on GitHub to manipulate polynomials in one variable in OCaml. It's pretty simplistic, and only works for polynomials whose coefficients are floating-point real numbers (so beware of roundoff error!), but I have found it to be useful in many quick-and-dirty situations. To grab it, use
git clone git://github.com/farr/ocaml-poly.git