J two - ohhh...

Jun 26, 2007 13:27

I've fixed the J Mandelbrot set program. Needless to say, the fix took precisely three characters :-)

[By the way, did anyone read that post? Does anyone else think J is interesting? Or did you just see a wall of punctuation and get scared off?]

The problem was that the bailout "absolute value greater than 2" test is performed after iterating the "plus-square" function 400 times. Now, iterates of plus-square can get big pretty fast: in fact, even within the small region under consideration, some points were overflowing after only eleven iterations. The first time it overflowed, it would return a complex number one of whose components was infinite: thereafter, it would return complex numbers whose real and/or imaginary parts were not-a-number (NaN). Now, the absolute value of NaN + i NaN is NaN: this is of course perfectly sensible, and exactly what NaN is for. Once an intermediate value becomes NaN, you want it to thread through your calculation, somewhat in the manner of Haskell's Maybe monad, not bothering you but getting on with the rest of the calculation (which is usually some large piece of array manipulation). A few data points that produce infinity or NaN can usually be ignored if the rest of the calculation goes OK.

Where things get odd is NaN's behaviour with respect to comparison operators. NaN is not equal to anything, even itself. NaN is not less than anything or greater than anything either. The upshot of which is that 2 < NaN returns False (or rather 0, since J uses Iverson's Convention). Arguably, 2 < NaN should return NaN (or maybe "Mu"?), but that wouldn't work in languages in which booleans and floats are a separate type. Hence, points which have overflowed are found to give absolute values less than 2, and thus (incorrectly) to be part of the Mandelbrot set.

So, the obvious fix would be to perform the bailout test earlier: to stop iterating as soon as |z| went over 2, and return true. Unfortunately, my J-fu is not yet strong enough for me to be able to do this without re-writing a big chunk of the code using an explicit while loop, which is not very Jish (explicit iteration is discouraged, not least because it's much slower). But there's an easier way. As well as not being greater than 2, NaN is also not less than 2. So, by changing the sense of the test round, and re-ordering the display characters (so points for which the test is true are displayed as '#', and those which test false are displayed as '.' rather than the other way round), we get the following:
{&'.#' @ (2:>|) @ ((+*:)^:200 0:) (18 %~ i:_20) j.~/ (28 %~ _59 + i.75)
...........................................................................
...........................................................................
...........................................................#...............
...........................................................................
........................................................#..................
......................................................###..................
.....................................................######................
.....................................................######................
......................................................####.................
............................................#....###############...........
............................................#####################.###......
.............................................#######################.......
.........................................###########################.......
..........................................###########################......
........................................###############################....
...............................#........##############################.....
...........................########....################################....
..........................###########..###############################.....
.........................#############.###############################.....
......................#..#############.##############################......
....###############################################################........
......................#..#############.##############################......
.........................#############.###############################.....
..........................###########..###############################.....
...........................########....################################....
...............................#........##############################.....
........................................###############################....
..........................................###########################......
.........................................###########################.......
.............................................#######################.......
............................................#####################.###......
............................................#....###############...........
......................................................####.................
.....................................................######................
.....................................................######................
......................................................###..................
........................................................#..................
...........................................................................
...........................................................#...............
...........................................................................
...........................................................................
Hurrah!
Total changes: one character to change '<' to '>', two to change '.#' to '#.' :-)

Now here's the interesting bit: the Haskell version seems only to work because of what I can only assume is a bug in GHC's handling of IEEE 754 special values. The Haskell version behaves exactly like the J version until it gets to the stage of taking the magnitude of NaN + i NaN, at which point (in Hugs) it dies with an "arithmetic overflow!" error (grrrr: that's what Infinity and NaN are for), or (in GHC) it returns Infinity (which is greater than 2, of course). G'huh? How does Infinity make sense here? It's not a number, it doesn't have a magnitude, much less an infinite one. Even more weirdly, sqrt (nan*nan + (nan*nan)) returns NaN as expected.

Yet again, I suffer pain because someone didn't read IEEE 754 properly. Which is the IEEE's fault for charging heavily for their standard documents, rather than making them freely available like they would have done if they actually cared about wide dissemination and interoperability. Grrrrr.

j, programming, chaos, haskell, computers, ieee754, maths, beware the geek

Previous post Next post
Up