Blog

2010-01-08: Clang self-hosts

The Clang frontends for LLVM have taken a big step forward; Clang now successfully self-hosts. This is a major improvement in the C++ support and I'm now very tempted to try it out with some of my projects.

2010-01-31: Cerebreduko

Me havos la okaziono ludar Big Brain Academy sur la Nintendo DS. La ludalo probas dilatar onua efektivocerebro par la ludaleto. Onu mesuras onua progreso en cerebrala pezo qua augmentas (assertos par rekalamanto). Me komencos a basa nivelo e sucesos elevar kelke facile.

Tamen, me demandas me quele evidenti existas suportanta ca aserturo. Naturale, me lektos la granda wikipedia qua dicos tala ludali havas ula ciencala evidento. Multe interesanta es la stimulanto par lernanta nova linguo es multe bona e exercas multa cerebroparti. Bona qua me lernas Ido :)

2010-01-30: La linguo Ido

Me deskovris recenta la linguo konstruktita Ido. Ido dizignis da Louis de Beaufront en 1907 kom Esperanto reformita do ol esas tre simpla e rapida lernar. La gramatiko esas kompleta asidua do la konjugi esas plu faculia kam angla.

Bona esas la povo konvertar inter la formi kom la adjektivi, substantivi e verbi. Ca permisas multa flexo e faciligas lernita.

Onu pouvas deskovrar a http://idolinguo.org.uk/. Me skribos hike en Ido tempope.

2010-01-04: Inverse searches between xdvi and vim

Two useful pages on the subject:

VIM wikia
a blog entry with example settings for VIM-LaTeX.

2009-11-18: Google's new programming language

A few days ago Google released a new programming language called go. Go was designed to combine the ease of use of interpreted languages with the speed and static typing of compiled languages.

As a language geek, I of course too a natural interest in it and had a quick look at the tutorial. Esthetically, the syntax appears to be a cross between C/C++ and Python, which is not surprising considering the design goal. Interestingly, there is only one type of loop construct: the for loop. A loop over a range 0..n-1 is written

for i:=0; i<n; i++ {
    // do something
}

This has a very familar feel for a C programmer, though the for parenthesis have been dropped (python-esq?). The brackets around the loop body are required.

Conditionals are similar, in that the parenthesis have vanished:

if i>5 {
    // do something
}

This looks a bit weird to me, but that's probably due to my familiarity with C.

One of the biggest differences with C is that there are pointers but no pointer arithmetic. The designers claim that this is to increase safety as referencing illegal addresses can never succeed incorrectly, and also to simplify the garbage collector (yes it has garbage collection). These are strong arguments against pointer arithmetic, but also somewhat unfortunate as some solutions can be succinctly expressed using pointer arithmetic and I like short programs.

The package system seems to be one of the biggest areas of improvement. Unlike C, there are no header files and each source file is part of a defined package. Other files can then use the package by simply importing it, much like Python. This is somewhat cleaner than the C/C++ header system.

Finally, I'll just mention that hash tables and strings are part of the language rather than being provided by libraries. This is nice as they are both extremely useful data structures, and having them built in allows syntactic simplification.

It'll be interesting to see how fast the language is taken up by the open source community.

2009-11-14: Sparse matrix libraries

I've been looking into sparse matrix libraries as I need to do some work with some sparse data. Initially I tried J, but unfortunately the sparse matrix implementation in J602 seems to be rather slow. Consequently, I've been searching for a library packages for other languages. So far, c++ seems to be the most promising language in terms of sparse linear algebra.

The first toolkit I tried was the sparse matrices present in uBLAS (part of Boost). Unfortunately, this implementation proved to be quite slow, and was not thread safe. Interestingly, g++ was FAR slower than the sun studio compiler.

I then found a blog entry which pointed to a couple of libraries: FLENS and Eigen. I decided to try the latter due to the praise in the blog entry, and also as it's well maintained as part of the KDE project. The API for accessing matrices is also quite good, with functions for applying things column/row/component wise. Unfortunately, these are only present for the dense matrix objects, and for sparse arrays one has to do manual looping. There are iterators of course so looping over non-zero entries is trivial.

Eigen seems to work very well, and is MUCH faster than boost. Unfortunately, the DynamicSparseMatrix object does not seem to compile using the sun compilers. I have thus been using only GCC to compile, but I wonder if the sun compiler would produce much more efficient code. I may try and patch the object to make it compile under Sun studio, but first I should get my actual experiment code working before trying to speed it up :)

Finally, the last sparse matrix library I came across is Getfem++. This library also looks interesting, but I have not tried it out in any way. It is obviously very powerful and does a lot more than I need. It's not a header based library, so installation isn't as easy as with uBLAS and Eigen.

2009-09-23: The birthday problem

Some time ago I solved the birthday problem in J. While not difficult, I thought I'd post the code here so it's recorded somewhere and I won't lose it.

The first solution I came up with was to use the Poisson distribution:

  f =: !@]%~^*^@-@[
  g =: 1:-2&!(f&0)@% 365"_
  2+1 i.~ 0.5<g 2+i.30
23

f is the Poisson PMF with $λ$ as the left argument and $k$ as the right. g calculates the probability of a pair out of $n$ people sharing a birthdate. The last line simply calculates the probability out to 30 people and finds the first exceeding 50%.

This is very cute mathematically, but there's another way to solve the problem and that is to use simulation:

  require'stats'
  dates =: 10000?@#365
  P =: (2+i.30) ([:mean(#@~.<#)\)"0 _ dates
  2+1 i.~ 0.5<P
23

The above code uniformly samples from the 365 possible birthdays (leap years ignored) a pool of 10000 people. It then uses this pool to form groups and estimates the expected number of people sharing a birthdate (the line defining P). The last line is essentially the same as in the Poisson case.

2009-09-17: Karpov, Kasparov, and Anand

There's some interesting matches coming up between Karpov and Anand and also Karpov and Kasparov. The former is scheduled to start on the 31/Oct and the latter on the 21/9.

2009-09-15: Piano repairs

My digital piano recently developed a rather annoying problem whereby certain keys started to stick. They would feel quite heavy when being depressed, and either return very slowly or get stuck at the bottom of the travel. I called a few people whom I located through Yamaha's service centre search on their website, but nobody wanted to come out to my apartment to fix the piano. Examining the keys revealed they were not sitting straight and were descending at an angle causing them to rub against the neighbouring key. As I couldn't shift it that easily, and I was reluctant to pay a technician to do so, I thought I'd just fix it myself. I suspected this was probably due to broken hinges and that it could be fixed easily by replacing the keys.

The piano itself is very well manufactured. Though I didn't have a service manual, it was trivial to work out how to disassemble the piano. There's plenty of space inside and things are easily accessible. All the wires are also neatly bundled and taped out of the way.

After a bit of work I had the keyboard module out, leaving the interior rather empty:

Empty interior

In this picture you can see the left mid-range driver attached to the bass reflex box. I have removed the right mid-range driver as it was actually damaged and not electrically connected (still sounds ok amazingly). I will need to source this part from Yamaha. The PCB in the middle with the electrolytic caps is the amplifier module. I think these are just chip amps, which are not necessarily bad (gainclone for example), but I didn't examine them much as I was more concerned with fixing the keys rather than upgrading the audio circuitry.

The keyboard module looks like this:

Keyboard module

I have removed the two problem keys, which you can see behind the module. Here's a close up:

Closeup keyboard module

You can see a thin metal spring protruding from the hole. This is used to return the key to the top position. The keys pivot on an axis at the very back of the keys. Just for interest, the back of the keyboard module looks like this:

Rear keyboard module

Each of the keys has a metal bar which provides the weight behind the key. When the key is depressed, the metal bar travels upwards and impacts on a bar of felt. Unlike a real piano, there is no escapement mechanism (and no need for one) so the bar remains in contact with the felt until the key is released.

Here's a closeup of the broken keys:

Broken keys

On each key, one side of the hinge point is broken which was causing the diagonal descent I noted earlier. Presumably one can order these specific keys from Yamaha to replace them. In the meantime I simply repaired these keys so I have something to play on while sourcing parts.

Repaired keys

Finally, here's a partially assemble picture of the piano with the top off but the keyboard module replaced:

Partially assembled interior

It turns out this is a known problem for some models of Clavinovas, and that Yamaha is replacing the entire keyboard module on affected models. I have yet to confirm this with Yamaha directly, but it would certainly be nice. Pity the mid-range drivers aren't being replaced too.

The temporary repairs to the keyboard module seem to have been quite successful. The keys no-longer stick and descend without any additional friction. To finish off the repairs, I'm going to try and source either a new keyboard module or replacement keys and the driver from Yamaha.

2006-09-12: Filing efficiency

Recently I purchased a pack of hanging folders for my filing cabinet. I like to have a general set of folders rather than subject specific folders, i.e., I prefer to have a set of folders related to the first letter of a keyword so I can rapidly file away items, and usually find them reasonably quickly. For example, I would file stuff related to my car under 'S' for "Subaru."

Unfortunately, the pack only contained 20 folders, though for some reason they did provide all 26 A-Z labels. Thus it's necessary to combine some letters, and I wanted to do this efficiently somehow. One approach is to figure out the probability of each letter occurring and then arranging groups to maximise the uniformity. The solution presented here is a greedy optimiser, so it may not find the global optimal.

First I load in some data to calculate letter frequencies. I don't have a corpus of keywords so I used Tolstoy's "War and Peace," which is probably a very bad choice.

require'files'
a =: a.{~97+i.26
corpus =: (#~e.&a) {.@>(' ',TAB,LF) cutopen dltbs tolower freads'2600.txt'

2600.txt contains the book and was obtained from Project Gutenberg. a is simply the 26 letters in lowercase. The contents of 2006.txt is read in, converted to lower case, and the first letter of every word extracted.

We can now calculate the word probabilities, and represent this in a boxed array.

require'stats'
P =: <"0 mean"1=corpus
T =: P,.<"0 a

Now, the idea is to continuously merge the two groups with the smallest probabilities. To merge two groups, the probabilities are added and the characters tracked.

M =: +&(>@{.) ; (/:~@,)&(>@{:)

This dyadic verb adds the scores and concatenates the reference strings. The concatenated string is sorted into ascending string order for readability later on.

The two groups with lowest probability can be found by simply sorting the structure T. As I placed the probabilities in the first column, this can be done simply using the verb /:~. We can then pop of the top two elements, merge them using M, and add the result back to the list. This is then repeated until we have 20 elements.

  n =: 20 
  /:~{:"1 (2&}.,[:M/2&{.)@:(/:~)^:(#>n"_)^:_ T 
+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+--+------+
|a|b|c|d|e|f|g|h|i|j|k|m|n|o|q|r|t|v|lw|psuxyz|
+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+--+------+

Looking at the groups, I think the choice of "War and Peace" was particularly bad. I should also potentially consider filing using the first two letters instead of the first letter only (with 40 folders). This would reduce the seek time when looking for items.

2009-09-05: An Atom feed

I'm proud to announce that this blog now has an Atom feed! It is unfortunately powered by J, not bash, so it is different from the rest of the CMS system I use to run this site.

I have not bugtested it much, so please let me know if you encounter any issues. In particular, I have not verified the MathML is working (assuming it's supported by the reader).

2009-08-30: Critical regions

Yesterday I spent the day writing a database for a project of mine. I initially started to use JDB, but I found it too buggy and the database kept getting corrupted. My database is a simple column database stored in files that are memory mapped. Most of the tables are populated only once and then only read from, but the table that stores results is frequently written to. Indeed, this was the trickiest point as there are many concurrent processes which may write to the table. Of course, one must guarantee that writes are atomic operations, so some form of locking is needed.

It turns out this was rather difficult. My first attempt was to use a lockfile and flock() to lock the file while writing. This worked perfectly in tests, but during actual runs the inserts were not atomic and I was losing many inserts. After a long time I discovered that flock() acts locally only on Linux NFS clients prior to 2.6.12.

The trick to file locking over NFS (at least with Linux) is to creat() an individual lockfile on the same filesystem (e.g. using the PID) and then link() to the main lockfile. If link() is successful, the process has the lock and can proceed into the critical region. Locks are freed by simply unlink()ing them.

2009-08-29: Interview with Arthur Whitney

While researching programming languages, I came across an interesting interview of Arthur Whitney. A lot of what he says I agree with, such as what makes elegant code, self describing code rather than comments, and how it's nice to start from scratch. I have not had much experience with Arthur's languages, but I have recently begun reading the book Q for Mortals. There's some interesting differences between Q and J and the integrated KDB+ seems very useful, and it's quite amazing that the interpreter can fit into the L1 cache!

2009-08-25: The Fibonacci sequence

The Fibonacci sequence is well known and is given by the simple recurrence relation $F_{n} = F_{n - 1} + F_{n - 2}$ . Some of the Project Euler questions involve the Fibonacci sequence, in particular large Fibonacci numbers. Thus, the question of how to calculate it efficiently arises.

There is a closed form solution known as Binet's formula, which is

$F (n) = \frac{ϕ^{n} - {(1 - ϕ)}^{n}}{\sqrt{5}},$

where $ϕ = \frac{1 + \sqrt{5}}{2}$ is the golden ratio. Calculating any particular number of the sequence therefore looks relatively straight forward, except when one wishes to calculate it exactly and for large $n$ .

To be able to calculate the number for large $n$ , we need to use arbitrary precision integers. These allow us to represent any $n \in ℚ$ exactly, but of course not $\sqrt{5}$ as it is irrational. The solution is to do calculations on the ring extension $ℚ [\sqrt{5}]$ .

We can do this in J by representing numbers in $ℚ [\sqrt{5}]$ by a vector $[a, b]$ that represents $a + b \sqrt{5}$ . Multiplying two such numbers can be done using convolution

ppr =: [:+//.*/

Applying such a verb to two vectors $[a, b]$ and $[c, d]$ gives $[a c, b c + d a, b d]$ representing the number $a c + (b c + d a) \cdot \sqrt{5} + b d \cdot 5$ . Of course, the last term can be combined with the first term, giving the vector $[a c + b d \cdot 5, b c + d a]$ . Defining this explicitly:

mul =: (0 1 0)"_ +//. (1 1 5)"_ * ppr

We need to take powers efficiently, which we can do through repeated squaring. Reusing a verb form defined in an earlier post gives us

pow =: 4 : 'mul/ mul~^:(I.|.#:y) x'

The golden ratio is now simply $ϕ = [\frac{1}{2}, \frac{1}{2}]$ , and $1 - ϕ = [\frac{1}{2}, - \frac{1}{2}]$ . Putting it all together gives

fib =: (1r2 1r2)&pow {:@:- (1r2 _1r2)&pow

Note that the final division is done by simply taking the $\sqrt{5}$ term after calculating the numerator, as the first term is necessarily zero. This function allows us to calculate large Fibonacci numbers quickly and exactly:

   (6!:2) 'echo #":fib 1e5'
20899
1.10964

So, $F_{10^{5}}$ has 20899 digits, and it took 1.10964s to compute.

2009-08-21: The Cooley--Tukey FFT

The discrete Fourier transform (DFT) is a well known algorithm from digital signal processing. It is used for transforming from the time domain signal representation to the frequency domain representation. This has several very important applications ranging from spectral filtering to polynomial multiplication (convolution), the latter being the motivation behind this post. Indeed, this post was inspired by a friend of mine who asked me some questions about how to do convolution quickly. Using the DFT, cyclic convolution can be carried out in $O (n log n)$ , whereas a direct calculation in the time domain takes $O (n^{2})$ .

Theorem 1

Let $F_{k} [x] = \sum_{n = 0}^{N - 1} x_{n} e^{- \frac{2 π i}{N} n k}$ and $F_{n}^{- 1} [X] = \frac{1}{N} \sum_{k = 0}^{N - 1} X_{k} e^{\frac{2 π i}{N} n k}$ denote the Fourier and inverse Fourier transform functions respectively. Furthermore, let $x$ and $y$ be two signals of length $N$ in the time domain with Fourier representations $X_{k} := F_{k} [x]$ and $Y_{k} := F_{k} [y]$ . Assume $x_{k} = 0$ and $y_{k} = 0$ for all $k < 0$ and $k \geq N$ . Then, $F_{n}^{- 1} [X \cdot Y] = {(x ⋆ y_{N})}_{n}$ with ${(y_{N})}_{n} := \sum_{p = - \infty}^{\infty} y_{n - p N}$ , where $⋆$ denotes convolution defined by ${(x ⋆ y_{N})}_{n} := \sum_{k = 0}^{N - 1} x_{k} \cdot {(y_{N})}_{n - k}$ .

Proof:

$F_{n}^{- 1} [X \cdot Y] = \frac{1}{N} \sum_{k = 0}^{N - 1} X_{k} \cdot Y_{k} \cdot e^{\frac{2 π i}{N} n k}$

$= \frac{1}{N} \sum_{k = 0}^{N - 1} (\sum_{l = 0}^{N - 1} x_{l} e^{- \frac{2 π i}{N} l k}) \cdot (\sum_{p = 0}^{N - 1} y_{p} e^{- \frac{2 π i}{N} p k}) e^{\frac{2 π i}{N} n k} by definition of X, Y$

$= \sum_{l = 0}^{N - 1} x_{l} \sum_{p = 0}^{N - 1} y_{p} \cdot (\frac{1}{N} \sum_{k = 0}^{N - 1} e^{\frac{2 π i}{N} k (n - p - l)})$

By Euler's identity, the last term is 1 only when $p = n - l (mod N)$ and zero otherwise. Thus, it can be replaced by the infinite sum $\sum_{m = - \infty}^{\infty} [n - l - m N = p]$ , where $[\cdot]$ denotes the Iverson bracket.

$= \sum_{l = 0}^{N - 1} x_{l} \sum_{p = 0}^{N - 1} y_{p} \cdot \sum_{m = - \infty}^{\infty} [n - l - m N = p]$

The summation over $p$ can be extended to $\pm \infty$ as $y_{l}$ is zero for $m \notin [0, N - 1]$ by assumption. Continuing:

$= \sum_{l = 0}^{N - 1} x_{l} \sum_{p = - \infty}^{\infty} \sum_{m = - \infty}^{\infty} y_{p} \cdot [n - l - m N = p]$

$= \sum_{l = 0}^{N - 1} x_{l} \sum_{m = - \infty}^{\infty} y_{n - l - m N}$

$= {(x ⋆ y_{N})}_{l}$

Q.E.D.

The consequence of this is that cyclic convolution between $x$ and $y$ can be computed in $O (n log n)$ time by $F^{- 1} [F [x] F [y]]$ , and discrete convolution can be calculated by appropriately padding the sequences. As a side remark, note that as Fourier transforms can be used to calculate convolutions, convolutions can be used to calculate Fourier transforms. This idea gives rise to Rader's and Bluestein's FFT algorithms, which we shan't discuss further.

Now, we'll consider how to do Fourier Transforms in J. The initial verbs presented here were initially presented by John Randall in the previous link, but we'll develop some extensions later on.

Recall the Fourier transform $F_{k} [x] = \sum_{n = 0}^{N - 1} x_{n} e^{- \frac{2 π i}{N} n k}$ . This is just a polynomial with coefficients $x_{n}$ evaluated at the roots of unity. We may therefore take advantage of the polynomial primitive p.

rou =: [: r. 2p1 * i. % ]       NB. Roots of unity
F   =: (] p. (+ @: rou @: #))"1 NB. Fourier transform

The first monadic verb simply calculates the roots of unity. The second monadic verb is the actual transform function, which evaluates the polynomial at the (conjugated) roots of unity. The inverse transform function $F_{n}^{- 1} [X] = \frac{1}{N} \sum_{k = 0}^{N - 1} X_{k} e^{\frac{2 π i}{N} n k}$ can be implemented much the same way:

FI  =:  ((# %~ ] p. (rou @: #)))"1 NB. Inverse Fourier transform

The two differences here are that the roots of unity are not conjugated, and the result is divided by the length. The forward and inverse transforms can be combined into one verb with a defined obverse as

FT  =: F :. FI NB. Combined Fourier transform verb

Polynomial multiplication can now be implemented by padding the signals, and multiplying in the Fourier domain:

pad =: ,. 0: $~ $
fpr =: (*/&.:FT)@:pad@:,: NB. Polynomial multiplication by FT

We can compare the results of this against direct polynomial multiplication:

require'numeric'
ppr =: +//.@(*/)              NB. Polynomial product
trim=:(#~ 1: 0} +./\.@(0&~:)) NB. Trims zero coefficents
  (ppr-:trim@clean@fpr)~i.5 
1

The clean function is needed here as there are some numerical errors that arise when transforming between the domains, and the trim function is required to removed the superfluous higher order zero components. So, there we have the Fourier transform and it's inverse defined in three lines of code.

Lets explore the running time of the FT:

time =: 3 : '(6!:2) ''FT i.y'''"0
  ]T =: time ]n=: 2^>:i.13
4.37927e_5 2.50244e_5 2.97546e_5 5.40161e_5 0.000141144 0.000489197 0.00143997 0.004702 0.0166071 0.0613451 0.241803 0.963937 3.8649
   ]B =: T %. n^/ i.3 NB. Fit a polynomial
0.000519551 _6.20351e_7 5.76555e_8
   B (+/ .*)~ n^/i.3  NB. View the polynomial curve
0.000518541 0.000517992 0.000518279 0.000524386 0.000558739 0.000716006 0.00138477 0.00413925 0.015316 0.0603405 0.241074 0.965277 3.86463

This appears to be scaling with $O (n^{2})$ as expected. This can be improved using the divide and conquer paradigm.

Theorem 2:

Let $x$ be a sequence of length $N$ such that $N = 0 (mod 2)$ . Then

$F_{k} [x] = F_{k} [x'] + e^{- \frac{2 π i}{N} \cdot k} F_{k} [x'']$

and

$F_{n}^{- 1} [x] = \frac{1}{2} (F_{n}^{- 1} [x'] + e^{\frac{2 π i}{N} \cdot n} F_{n}^{- 1} [x''])$

where $x_{k}' := x_{2 \cdot k}$ and $x_{k}'' := x_{2 \cdot k + 1}$ are two subsequences of $x$ consisting of the even and odd indices respectively.

Proof:

$F_{k} [x] = \sum_{n = 0}^{N - 1} x_{n} e^{- \frac{2 π i}{N} \cdot n \cdot k}$

$= \sum_{l = 0}^{\frac{N}{2} - 1} x_{l} e^{- \frac{2 π i}{N} \cdot (2 l) \cdot k} + \sum_{m = 0}^{\frac{N}{2} - 1} x_{m} e^{- \frac{2 π i}{N} \cdot (2 m + 1) \cdot k}$

$= \sum_{l = 0}^{\frac{N}{2} - 1} x_{l} e^{- \frac{2 π i}{N} \cdot (2 l) \cdot k} + e^{- \frac{2 π i}{N} \cdot k} \sum_{m = 0}^{\frac{N}{2} - 1} x_{m} e^{- \frac{2 π i}{N} \cdot (2 m) \cdot k}$

$= F_{k} [x'] + e^{- \frac{2 π i}{N} \cdot k} F_{k} [x''] .$

Similarly,

$F_{n}^{- 1} [x] = \frac{1}{N} \sum_{k = 0}^{N - 1} x_{n} e^{\frac{2 π i}{N} \cdot n \cdot k}$

$= \frac{1}{N} \sum_{l = 0}^{\frac{N}{2} - 1} x_{l} e^{\frac{2 π i}{N} \cdot (2 l) \cdot n} + \frac{1}{N} \sum_{m = 0}^{\frac{N}{2} - 1} x_{m} e^{\frac{2 π i}{N} \cdot (2 m + 1) \cdot n}$

$= \frac{1}{2} (\frac{2}{N} \sum_{l = 0}^{\frac{N}{2} - 1} x_{l} e^{- \frac{2 π i}{N} \cdot (2 l) \cdot n} + \frac{2}{N} e^{- \frac{2 π i}{N} \cdot n} \sum_{m = 0}^{\frac{N}{2} - 1} x_{m} e^{- \frac{2 π i}{N} \cdot (2 m) \cdot n})$

$= \frac{1}{2} (F_{n}^{- 1} [x'] + e^{\frac{2 π i}{N} \cdot n} F_{n}^{- 1} [x'']) .$

Q.E.D.

Thus, the transform has been decomposed into two smaller transforms, for which the same process can be repeated if it the smaller sequences are of even length. The $e^{- \frac{2 π i}{N} k}$ factors are called the twiddle factors. This is the radix-2 decimation in time form of the Cooley--Tukey Fast Fourier Transform (FFT) algorithm.

One final thing to note is that $F_{k} [x']$ for $k \geq \frac{N}{2}$ is equal to $F_{k - \frac{N}{2}} [x']$ and similarly for $x''$ , so these sums need only be calculated for $\frac{N}{2}$ values.

To implement this, we'll use the prior implementation of the FT as a fallback for when $N \neq 0 (mod 2)$ . Note that if $N$ is a power of 2, then this process can be repeated down to a subsequence of length 1, for which the Fourier transform is simply the identity function. Thus, applying the radix-2 FFT to such a sequence is very efficient. So lets write the FFT verbs:

NB. Fast Fourier transform (Cooley--Tukey)
FF  =: (([:+/ ((#$1:),:#$+@:rou@:#) * (2:|i.@#) ,~@:$:/.])`F @. (2:|#))"1 
FFI =: (([:-:@:+/ ((#$1:),:#$rou@:#) * (2:|i.@#) ,~@:$:/.])`FI @. (2:|#))"1
FFT =: FF :. FFI

The verbs FF and FFI do the forward and inverse transforms respectively. Both of these first check if $N \neq 1 (mod 2)$ and revert to the FT if it does. If not, the key primitive /. is used to recursively apply itself to the even and odd components. The twiddle factors (which are roots of unity again, conjugated in the forward case) are then multiplied in and the final summation calculated.

Some sanity checks comparing the FT and FFT verbs:

   +/*: (FT - FFT) i.2^13
 7.07245e_10j2.37437e_10
   +/*: (FT inv - FFT inv) i.2^13
1.05388e_17j_3.53809e_18

and some timing results:

  (6!:2) 'FT i.2^13'
3.88696
  (6!:2) 'FFT i.2^13'
0.276201

It's clear that the FFT in this case is significantly faster than the FT, but of course this is the best case scenario as the length of the input was a power of 2. Of course, one could always pad sequences out to a power of 2 boundary if appropriate.

So there we have it, the FT in two lines of code and the FFT in three.

2009-08-17: The AKS primality test

Recently, Terence Tao wrote about the Agrawal-Kayal-Saxena (AKS) primality test. The AKS test was discovered in 2002 and is the first deterministic polynomial-time algorithm to determine primality of a given number. The concept and proof of the AKS test is quite beautifully simple, but I shall not discuss the proof here as Terry has a nice outline on his blog. Instead, I thought it'd be instructive to implement the algorithm in J, which is a language I'm currently learning.

So to begin, some methods for handling polynomials are needed. In J, polynomials are most easily represented by a vector of coefficients, and there are some simple primitives for evaluating polynomials in such a representation. However, there are no primitives for addition, multiplication, or – most importantly – exponentiation. Thus we begin by defining some simple polynomial operations

sub  =: -/@,:
add  =: +/@,:
trim =: (#~ 1: 0} +./\.@(0&~:)) NB. Trims zeroes
ppr  =: +//.@(*/) NB. Polynomial product
ppow =: 4 : 'ppr/ ppr~^:(I.|.#:y) x' NB. Polynomial exponentiation by squaring

The first two diadic verbs simply add and subtract two given polynomials. The third verb is monadic and trims the trailing zeros (i.e., the higher order terms with a coefficient of 0). The second diadic verb ppow exponentiates a polynomial using repeated squaring. It's actually quite amazing how easy it is to define using just a few primitives. There is an essay on repeated squaring in J if more details on the repeated squaring method are required.

Of course, we need to do polynomial exponentiation on the ring $F_{p} [X] / X^{r} - 1$ , where $F_{p} [X]$ is the ring of polynomials of variable $X$ over the finite field $F_{p}$ of $p$ elements, so some more verbs will be needed to do this efficiently. First, we'll need to be able to reduce polynomials to the ring $(ℤ / n ℤ) [X] / (X^{r} - 1)$ :

red =: {.@[ | (({:@[|i.@:#@:]) +//. ]) NB. Polynomial reduction

This dyad takes a polynomial $f \in F_{p} [X]$ as it's left argument, and a tuple $(n, r) \in ℤ \times ℤ$ as the right argument and reduces it to the ring $(ℤ / n ℤ) [X] / (X^{r} - 1)$ . Using this reduction, polynomial exponentiation on the ring can be defined:

pmodpow =: 4 : 'y red (y&red@[ ppr y&red@])/ ([:ppr~ y&red)^:(I.|.#:{.y) x'

This dyad takes a polynomial $f$ as it's left argument, and a tuple $(n, r)$ as it's right argument, and calculates $f^{n} (mod n, X^{r} - 1)$ . We can also easily define integer exponentiation on the ring $ℤ / n ℤ$ much the same way:

modpow =: 4 : '({:y) | */ ({:y) |*~^:(I.|.#:{.y) x'

Here, modpow takes $x \in ℤ$ as it's left argument and a tuple $(n, r)$ on the right and returns $x^{n} (mod r)$ . We now have all the tools needed to implement the AKS algorithm.

aks =: monad define
  NB. Check if perfect power
  k =. i.@:>:&.(p:inv) <. 2^. y
  roots =. k %: y
  if. +./ (=<.) roots do. 0 return. end.

  NB. Calculate smallest r such that O_r(n)>log^2 n
  r =. 2x
  found =. 0
  while. -.found do.
    found =. 1
    for_k. i.&.(-&2x) <.*: ^. y do.
      if. 1=y modpow k,r do.
        r =. r + 1
        found =. 0
        break.
      end.
    end.
  end.

  if. *./ 1~:y +. i.&.<: r do. 0 return. end.

  for_a. x:i.&.<: <. (%:r) * ^.y do.
    lhs =. trim (a,1x) pmodpow y,r
    rhs =. trim (,:a) add (=i.&>:) r|y
    p =. y|lhs sub rhs
    if. +./ p  do. 0 return. end.
  end.
  1
)

Here's an example of using AKS to filter the first 100 integers:

  (#~aks"0) >:i.100
2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97

2009-08-11: Convex hulls

I recently solved Project Euler Problem 252 which is an interesting geometrical challenge involving convex polygons. While I won't describe my solution here, which incidentally scales with $O (n^{3}),$ I will discuss the related problem of finding the convex hull of a set of points.

Of course, this problem has been well studied and there exists a $O (n)$ online algorithm for finding the convex hull of a simple polychain. This algorithm is called Melkman's Algorithm.

The algorithm first begins by ordering the first three points of the chain into a convex anticlockwise chain. These three points are then added to a deque, where the head and tail of the deque share the same point. I.e., if the three points are labelled $a$ , $b$ , and $c$ , and $a b c$ forms a convex anticlockwise chain, then the initial deque is $a b c a .$ The remaining points are then handled iteratively:

Let $v$ be the next point.
Ignore $v$ if it falls within the convex hull. As we are dealing with a simple polychains, this can be checked by simply verifying it lies in the convex cone formed by the top two and bottom two points of the deque.
If $v$ lies outside the cone, pop points from the top and bottom of the deque until $v$ lies within the new cone, then add $v$ to the top and bottom of the deque.

Once all points have been processed, the deque contains the convex hull. Note that though this algorithm scales with $O (n),$ if a sort is required to ensure the given points are arranged as a simple polychain then it scales with $O (n log n) .$

Here's a J implementation:

det =. (-/ . *)
left =: (0: < [:det"2[:|:"2-"_ 1) f. NB. Tests if a point is left of a line
right =: -.@left

init =: monad def '(,~{:) ]`(({.,|.@}.)`(|.@}:,{:)@. ((0 1 2)"_ -: /:))@. ({. right }.) y'

melkman =: monad define
  y =. /:~y NB. Ensure we have a simple polychain by sorting
  D =. init 3{.y NB. Create initial deque

  for_v. 3}.y do.
    NB. Check if v is not inside convex cone
    if. -.(v left 2{.D) *. (v left _2{.D)  do. 
      while. v right _2{. D do. D =. }: D end. NB. Pop points from bottom of deque until v lies to the left
      D =. D,v NB. Add v to the bottom of the deque

      while. v right 2{.D do. D =. }. D end. NB. Pop points from the top until v lies to the left
      D =. v,D NB. Add v to the top of the deque
    end.
  end.

  D
)

and an example which generates some random points and plots the convex hull:

require'plot'
pline =: dyad define
pd'reset'
pd'type line'
pd<"1 |:x
pd'type marker'
pd'markers circle'
pd<"1|:y
pd'show'
)

P =: 100 2 ?@$ 0
(pline~melkman) P

2009-08-05: The start of a new blog

I have once again been inspired to keep a blog, although this time it will be different from my last one. To start with, it'll be in English most of the time. Keeping a blog in French is an interesting way to learn French, but my French isn't that good anyway so it was probably more embarrassing then helpful.

Secondly, I plan to include my research interests in this blog. My last blog was purely a personal blog, and as a result I didn't really have that much to blog about. This time, this blog will be as much for me as for anybody else, and I'll use it to remember things I've read, seen, or interesting ideas.

For this first blog entry I'll spend some time on some videos I have been watching recently. The first of these was an interesting talk by Gregory Chaitin given at CMU. His lecture briefly touches upon problems with infinite sets in Cantor's theory, David Hilbert's famous problems, Russell's paradoxes, and the works of others such as Godel and Turing. He is a talented lecturer and his presentation is fun to watch.

The second set of videos which I have found are two series of lectures given by Richard Feynman. The first is the Character of Physical Law lecture series that was given by Feynman in 1964 at Cornell. I have only watched the first episode, which discussed the law of gravitation. Feynman has a very good lecturing style which is engaging, and his use of the blackboard during his talks is inspiring given the prolific nature of powerpoint presentations these days.

The second series is the Douglas Robb Memorial Lectures given by Feynman at the University of Auckland in 2007. I have not yet watched these videos, but no doubt Feynman's style shines as it did in 1964.

Finally, the last video I shall mention in this entry is the presentation Phase transition phenomenon in Compressed Sensing given by Jared Tanner at the recent SMLS'09 conference. The presentation discusses compressed sensing, and in particular studies sparsity recovery using L-1 regularisation for Gaussian random matrices. Lower bounds on the phase transitions for Gaussian random matrices satisfying the restricted isometry property are given, but also other alternative bounds based on geometry are given. The geometrical bounds are obtained by counting faces of random projections of the L-1 ball, which is a convex polytope. The paper by David Donoho and Jared Tanner Counting Faces of Randomly-Projected Polytopes when the Projection Radically Lowers Dimension covers this face counting approach, but I'm yet to read it.