Archive for August, 2009

Improved geohashing algorithm

Friday, August 14th, 2009

Geohashing is an awesome idea (though I've never actually gone to any of the hash points). However, it does suffer from some problems - often the chosen location will be too far away, or in an inaccessible area. It is also slightly biased towards the poles (though because few people live in the polar regions this isn't a particularly big deal for most people). It would be nice to have an improved algorithm which solves these problems.

We'll start the same way, by hashing the date and that date's DOW opening to get a set of 128 random bits R.

Given any area A we want to find a function (which we'll call f(A)) mapping a set of bits R onto a point x, (f(A))(R) = x. We would also like to have the property that if B is a subset of A containing x and if (f(A))(R) = x then (f(B))(R) = x. So if the New York State meetup happened to be in New York City, all the people who were just going to the New York City meetup would be there too. Another desirable property would be if the points (f(A))(R) were evenly spread throughout A.

Here is an algorithm that I think fulfills both properties.

First, use R to pick a random point on the globe, taking 64-bit random floating point fractions u and v in [0,1) from R in the usual way and then transforming them into longitude and latitude as:
\displaystyle long = 360frac(u+a(N))-180
\displaystyle lat = \frac{180}{\pi}cos^{-1}(2frac(v+b(N))-1)
Where frac takes the fractional part, leaving a number between 0 and 1 (one could also use XOR instead of add and frac) and N is the smallest integer > 0 that causes the resulting coordinates to end up inside A.

a and b are functions that find points that are furthest away from any tried so far. They have binary expansions as follows:

N  a(N)  b(N)  c(N)     d(N)
0  0.000 0.000 0.000000 0.000000
1  0.100 0.100 0.110000 0.100000
2  0.100 0.000 0.010000 0.010000
3  0.000 0.100 0.100000 0.110000
4  0.010 0.010 0.001100 0.001000
5  0.110 0.110 0.111100 0.101000
6  0.110 0.010 0.011100 0.011000
7  0.010 0.110 0.101100 0.111000
8  0.010 0.000 0.000100 0.000100
9  0.110 0.100 0.110100 0.100100
10 0.110 0.000 0.010100 0.010100
11 0.010 0.100 0.100100 0.110100
12 0.000 0.010 0.001000 0.001100
13 0.100 0.110 0.111000 0.101100
14 0.100 0.010 0.011000 0.011100
15 0.000 0.110 0.101000 0.111100
16 0.001 0.001 0.000011 0.000010

And so on, to as many binary places as you need. If you interleave the bits of a(N) and b(N) you get c(N), which looks like d(N) but with each even bit XORed with the bit to its left. d(N) is the binary expansion of the integers but with the bits in the reverse order (and flipped to the other side of the binary point).

The sequences a and b together are related to the Bayer matrix which describe a 2D ordered dither. If you want 65 shades of grey but only have a grid of 8x8 pixels, each of which can only be black or white, the most regular patterns with n black bits are described by those that have 8a(N) and 8b(N) black if N

Component pattern

Thursday, August 13th, 2009

Sometimes, your object hierarchy is determined by the kinds of data you're throwing around - if you have a pile of pieces of information about each customer in a database, it's obvious that you should have a Customer object. But if this is the only rule you use for deciding when to make a class, you'll often end up with enormous classes. This tends to happen for the object representing the program itself. A program tends to have a lot of things that there are only one of, and if you stick all these things in one big object then its unwieldy.

So you break up your Singleton object into components. You can do this in whatever way makes most logical sense for your program - perhaps a set of routines and data for allocating memory go in one object, and a set of routines and data for dealing with the program's main window go in another. Then you need a master object responsible for wiring all these components together. The construction of this object may be somewhat non-trivial:

In the master constructor, each of the components is constructed. However, we don't tell the components about the master yet because it's constructor is incomplete.

Once the constructor is complete, we need to loop through all the components again and "site them" (give them a pointer to the master that they can use to obtain pointers to other components they need).

Depending on how complex the interdependencies between components are, more initialization steps may be needed. The order of the components may also make a difference.

This doesn't have to be used for just the singleton - any sufficiently complex object can benefit from being broken up this way.

I looked for this pattern in the literature but didn't see it anywhere. However, it seems to be a reasonably common thing to do so I'm recording it here.

The programmer cooperative

Wednesday, August 12th, 2009

If I suddenly found myself with lots of money to spare, one thing I'd like to do is hire some talented programmers to work on whatever they wanted to work on. I suspect that there are lots of talented programmers out there who have ambitious projects they'd like to take on but which they just don't have time for because they have full time jobs. Most of those projects probably have no commercial application, but it's very difficult to tell in advance what is likely to make money and what isn't. If you get enough of these projects together, at least some of them are bound to make money.

Any profits made could be shared with the entire cooperative, or used to hire more programmers. Sort of like a business and employment cooperative except for programmers rather than entrepreneurs. Ideally once it got big enough it would become self-sustaining and then undergo exponential growth. Part of the nature of software is that it requires very little capital investment (really just the programmers' salary). Another part is that because it can be copied indefinitely and exhibit network effects, once a piece of software is written it can have great influence and importance. So a successful software project can make much more money than it costs, and one only needs a small percentage of projects to be financially successful in order to sustain the entire system.

I imagine that a very large number of programmers would like to join such a cooperative - what could be better than getting paid and not having anyone tell you what to work on? So one would need some method for choosing between a very large number of applications. Let's suppose for a moment we have some infallible way of weeding out people who, given the choice, would rather go skiing or watch TV all day than write computer programs of any sort. Suppose also we have a way of weeding out those programmers who are actually working on something for someone else and want to get paid for it twice.

One way to choose amongst the remainder would be for them each to submit a few ideas of some things they would like to work on, and how much they want to be paid. The cooperative would then pick those ideas that have the most "bang per buck" (i.e. are most interesting, not too expensive and have some chance of commercial success). The member might change their mind about what they work on, but the initial idea would give some idea about what sort of things they're interested in.

Every few months, each member would give a demo of what they've been working on to the rest of the company. Or, alternatively, make regular blog posts about their progress. This gives some cross-pollination of ideas, and may lead to collaborations within the cooperative. Such collaborations can be very flexible and informal, ending as soon as they stop being fun for any of the parties involved.

CGA Hydra

Tuesday, August 11th, 2009

A while ago, Trixter challenged me to figure out if it was possible for a CGA card with both composite and RGB monitors attached to it to display a different image on each display. At first I thought this was impossible because the composite output is just a transformation of the RGB output - the RGB output contains all the information that the composite output contains.

But that reasoning only works if you're close up. If you stand back sufficiently far from the screens, adjacent pixels will blur into each other so this is no longer necessarily true. Suppose we have a pattern that repeats every 4 high-resolution pixels (or half an 80-column character, or 1/160th of the screen width, or one colour carrier cycle) and we stand sufficiently far back that this looks like a solid colour. On the RGB monitor this will just be an average of the 4 colours making up the pattern. So, for example, black-black-white-black and white-black-black-black will look the same on the RGB monitor, but they will look different on the composite monitor because these two patterns have different phases with respect to the color carrier, so they will have different hues.

That explains how we can get details on the composite monitor but not on the RGB monitor, but what about the other way around? This is a bit more complicated, because it requires knowing some more details about how the CGA generates (non-artifact) colour on the composite output. For each of the 8 basic colours (black, blue, green, cyan, red, magenta, yellow and white) there is a different waveform generated on the card. The waveform for the current beam colour is sent to the composite output. The waveforms for black and white are just constant high and low pulses, but the waveforms for the 6 saturated colours are all square waves of the colour carrier frequency at different phases. The green and magenta lines switch between high and low on pixel boundaries, the other 4 at half-pixel boundaries (determined by the colour adjust trimpot on the motherboard).

What this means is that if you're displaying a green and black or magenta and black image, the pixels are essentially ANDed with this square wave. The pixels corresponding to the low parts of these waves have no effect on the composite output. So you can use these pixels to make the image on the RGB monitor lighter or darker whilst having no effect on the composite image.

Here's what the finished result is supposed to look like (with another image on an MDA display as well):

Note that I've allowed the composite image to show through on the RGB monitor a little in order to improve contrast.

Floating point is not evil

Monday, August 10th, 2009

In response to this:

Floating point is not evil and it is deterministic, but you need to know what the values you're working with actually are. Basically, floating point is like scientific notation, except with 2s instead of 10s. In other words, instead storing numbers like 1.234*104 it stores numbers like 1.00110100102*210. It's actually stored as a pair of numbers, a mantissa and an exponent. The number of bits in the mantissa is fixed (it's 23 bits for single precision and 52 for double) and the leading "1" is implied. Each number representable as an IEEE floating point constant has exactly one representation except for zero (+0 and -0 have different bit patterns for complicated reasons). There are complications (denormals, NaNs and infinities) which can usually be ignored and which I won't go into.

Floating point numbers are handy for lots of purposes but they do have a couple of problems.

The first problem is that floating point numbers are inefficient. They are very quick with today's hardware but consider what you'd do if neither floating point hardware nor floating point libraries were available. For most applications, you'd use fixed point numbers - you'd store an integer and it would be implied by the type of number you're working with that the actual numeric value is obtained by dividing this integer value by 2n for some n. For most purposes you probably wouldn't store that n with each integer - all your numbers have the same number of significant digits. For example, if you're writing a graphics program you might decide that units of 1/256 of a pixel width are always enough, so n would always be -8. When writing floating-point programs, most programmers don't do this calculation to figure out what the precision needs to be, they just use single precision floating point or switch to double if that isn't precise enough. While constant precision is preferable for a general purpose calculator, most actual applications are better served by constant resolution.

The other problem is that sooner or later you'll run out of precision. If you're plotting Mandelbrot sets, sooner or later you'll zoom in far enough that adjacent pixels have complex numbers with the same floating-point representation. If you're using FFTs to multiply big integers, sooner or later you'll want to multiply integers so large that floating-point numbers won't have sufficient precision. If you're using hardware floating point, this is quite difficult to solve (you need to find or write a big-float library) and will cause a big speed hit, so most people will give up at that point. However, if you're already using fixed point bignums, it's just a question of adding another digit.

Memory handling for fractal program

Sunday, August 9th, 2009

I have been poking again at my fractal plotting program recently. The latest version is available here (source). I've added multi-threading (so all available cores will be used for calculation) and Direct3D rendering (which makes zooming much smoother and allows us to rotate in real time as well with the Z and X keys - something I've not seen any other fractal program do). There is still some jerkiness when zooming in or out (especially past a power-of-two boundary) but it's greatly reduced (especially when there's lots of detail) and with a little more work (breaking up the global matrix operations) it should be possible to eliminate it altogether.

The linked grid data structure in the last version I posted worked reasonably well, but it used too much memory. At double-precision, every block needed 72 bytes even after it had been completed. This means that an image that covered a 1600x1200 screen with detail at 16 texels per pixel for anti-aliasing would use more than 2Gb. I could use a different block type for completed points but there would still need to be at least 4 pointers (one for each direction) at each point, which could take 500Mb by themselves.

So I decided to switch to a data structure based on a quadtree, but with some extra optimizations - at each node, one can have not just 2x2 children but 2n by 2n for any positive integer n. When a node contains 4 similar subtrees (each having the same size and subnode type) it is consolidated into a 2n+1 by 2n+1 grid. If we need to change the type of one of the subnodes of this grid, we do the reverse operation (a split). In this way, if we have a large area of completed detail it will use only 4 bytes per point (storing the final number of iterations). This works much better - memory usage can still be quite high during processing (but seems to top out at ~700Mb rather than ~2Gb for the linked grid data structure) and goes right back down to ~200Mb when the image is completed (only about 35Mb of which is the data matrix, the rest is Direct3D textures).

One disadvantage of the grid tree method is that it doesn't naturally enforce the rule preventing 4 or more blocks being along one side of one other block. That rule is important for ensuring that we don't lose details that protrude through a block edge. However, this turns out to be pretty easy to explicitly add back in. This gives a more efficient solid-guessing method than I've seen in any other fractal program (by which I mean it calculates the fewest possible points in areas of solid colour).

This data structure was surprisingly difficult to get right. The reason for this is that with all the splitting and consolidating of grids, one must avoid keeping grid pointers around (since they can get invalidated when a split or consolidation happens in a nearby part of the tree). That means that many algorithms that would be natural to write recursively need to be written iteratively instead (since the recursive algorithm would keep grid pointers on the call stack).

Another tricky thing to get right was the memory allocation. My first attempt just used malloc() and free() to allocate and deallocate memory for the grids. However, the way grids are used is that we allocate a whole bunch of them and then deallocate most of them. However, they are quite small and the ones that don't get deallocated are spread throughout memory. This means that very little memory can be reclaimed by the operating system, since nearly every memory page has at least one remaining block (or part of one) on it.

The problem here is that the malloc()/free() interface assumes that a block, once allocated, never moves. However, grids turn out to be quite easy to move compared to general data structures. This is because the grid tree is only modified by one thread at once (it must be, since modifying it can invalidate any part of it) and also because we know where all the grid pointers are (we have to for split and consolidate operations). Only a dozen or so different grid sizes are currently used in practice, so it's quite practical to have a separate heap for each possible grid size. If all the objects in a heap are the same size, there's an extremely simple compaction algorithm - whenever we delete an item, just move the highest item in the heap down into the newly free space. So a heap consists of a singly-linked list of items. These are allocated in chunks of about 1Mb. Each chunk is contiguous region of memory, but chunks within a heap don't have to be contiguous. This allows the heaps to interleave and make the best use of scarce (on 32-bit machines) address space.

Sounds simple enough but again there turned out to be complications. Now it's not just split and consolidate that can invalidate grid pointers - now deleting a grid can also invalidate grid pointers. This means more recursive algorithms that need to be made iterative (such as deleting an entire subtree) and certain operations need to be done in a very specific order to avoid keeping more grid pointers around than are necessary. I now understand why such techniques are so rarely used - they're really difficult to get right. I do think that a general garbage-collected heap would be slower and use more memory than this special purpose allocator, but I'm not sure if difference would be large enough to make all this effort worth it. Perhaps someone will port my code to a garbage-collected language so the performance can be compared.

Sound for fractal program

Saturday, August 8th, 2009

I'd like to be able to add some kind of sound output to my fractal program. What I'd like to do is generate a click every time a point is completed. So, when it's working on a difficult area it would emit a clicking sound but when it's finishing a lot of points it would generate white noise. The principle is very similar to that of a geiger counter's audio output.

This is trivial to do under DOS, but in Windows it is surprisingly difficult. Ideally we want to generate a waveform at 44.1KHz. When we finish a point, we flip the current sample in the audio buffer between its maximum and minimum positions. The problem is determining what the current sample should be (or, alternatively, moving a pointer to the next sample 44,100 times per second). What we need is a high-resolution timer, but none of the high-resolution timers seem to be particularly reliable. timeGetTime() is too low-resolution. QueryPerformanceCounter() has bugs on some systems which cause it to skip backwards and forwards a few seconds on occasion. The RDTSC CPU instruction isn't guaranteed to be synchronized between CPUs, and its count rate may vary due to power saving features on some CPUs. Another possibility is reading the playback cursor position from DirectSound, but that is unreliable on OSes before Vista.

What I may end up doing is faking it - measure the average number of blocks completed over some small period of time and use that to generate the clicks via a Poisson distribution.

Public medicine and medical research

Friday, August 7th, 2009

The US pays twice as much for healthcare per person as the average western country, and supplies a lower quality of care. So where does all that extra money go to?

Some of it goes to administrative overhead (it costs a lot of money to pay the people at your insurance company who are looking for excuses to deny you coverage). But a larger fraction goes towards the profits of the insurance companies. People invest in these companies (often as part of a 401k or other retirement plan) and get back a profit (at the expense of those who pay insurance premiums and consume healthcare).

But not all of a company's profits go immediately back to its investors - any reasonable company will spend some portion of their profits on research into ways to make even more money in the future. In terms of medical insurance companies, this generally means medical research - research into new techniques, medicines and gadgets that improve the state of the art and make us all healthier in the long run.

So if the US switches to public medicine, what happens to the state of medical research? If the massive profits go away, will the research money dry up too? I have heard concerns that socialized medicine will cause medical research to grind to a halt, leaving the entire world with 2009-level medical technology for the forseeable future. I don't think it will be quite that bad - there are other sources of medical research money than US insurance companies (the US government and other governments).

I hope that in the process of transforming the country to socialized medicine, the government continues the investment in medical research and changes its focus a bit. I have the impression that the direction of medical research is set at least somewhat by the treatments that will make the most money rather than the treatments that are medically best. Drug companies pursue research on drugs that they can patent, not new uses for existing drugs on which the patents have expired. Cures are pursued at the expense of prevention and finding causes. With a change in focus we can hopefully get more effective research for less money overall.

Why insurance is the wrong model for healthcare

Thursday, August 6th, 2009

Insurance is a good model for things like car insurance and home insurance because the duration of a claim is very short - one never has to change one's insurance provider for reasons beyond one's control during a car crash or an earthquake. But a medical condition can last for many years or even a lifetime, during which there can be many individual claims.

Once you have a medical condition that is likely to take many years to resolve and be extremely expensive, you're beholden to your insurance provider, since you are a liability for them. Should you be allowed to change insurance providers? Okay, most of us can't change providers anyway without spending a lot more (to avoid paying tax on our health premiums we're stuck with our employer-provided health plans) but suppose that problem was fixed?

Let's suppose the answer is no, you shouldn't be allowed to change (at least not without taking a cripplingly large increase in premiums turning it from insurance into a payment plan) - no business should be required to take on a customer that is guaranteed to lose them money. Then the "free market" aspect of health insurance is denied to those people who need it most - it's not a free market when you can't make a choice.

So let's suppose the answer is that you should - that insurers should not able to deny access to insurance based on pre-existing conditions (which really ought to include genetic predispositions). Then insurers have to factor in the costs of such patients to all their customers' premiums. That doesn't sound so bad at first glance but consider the consequences - suppose that insurance company A had low premiums but covered very little and insurance company B had high premiums and covered everything. Then all the healthy people get insurance from A and switch to B when they get sick. B's prices then rise to cover the cost of their coverage and what you essentially end up with is equivalent to not having insurance at all. This means that government has to regulate what coverage is provided as well as who it is provided to. All insurance companies will then provide that level of coverage but none will provide any more. We've gone from having no choice to having no real choice. Either way we fail to have a free market for healthcare.

Socialized medicine aka single-payer healthcare aka politically controlled medicine has problems too. In the UK, different districts decide which treatments they will offer, resulting in a "postcode lottery" for some treatments. Waiting lists of years are not unheard of for treatments that are not imminently life threatening, and doctors and nurses are overworked and underpaid compared to their counterparts in the US. However, the average Brit pays about half what the average US resident pays in healthcare costs and gets better coverage on average. There is also no worry about medically-induced bankruptcy, and no worries about one's insurance company refusing coverage (due to a "pre-existing condition") once a serious medical condition is found.

I disagree that single-payer healthcare gives patients an incentive to over-consume. It should not be the insurance company that gets to decide that - doctors should be ones to decide which treatments are medically necessary. And if some doctors are over-prescribing compared to other doctors in their geographical area or speciality, that is something that should be investigated.

Another criticism of single-payer healthcare is that providers have no incentive to compete on price - the government will have to pay whatever they ask for critical treatments. But there's nothing special about single-payer healthcare there - there's nothing stopping the government from acting the same way as insurance companies to do to keep a check on prices. In fact the government has more power than an insurance company there - the insurance company only has the power to drop an overcharging provider from their preferred providers list while the government has the ability to take away their medical license.

Yet another criticism is that if the government gets to decide what treatments are on offer, people don't have any choice. Well, most of us don't have any choice anyway under the insurance model but at least with the government you can vote them out if they aren't doing a good job with healthcare. Yes, socialized medicine means a government bureaucrat between the patient and the doctor, but better an elected bureaucrat than an insurance-company bureaucrat that one has no control over.

378Kb disk hack for Fractint 15.1

Wednesday, August 5th, 2009

I first heard of Fractint soon after I got interested in Fractals. I knew it was very fast and had lots of fractal types and features and I wanted to run it on the PC1512 (the only PC I had regular access to at the time). I bought a "Computer shopper" magazine which had a copy of Fractint 15.1 on the cover disk. Only the 3.5" disk was available, so I had to get a friend who had both drives to copy it onto 5.25" disks for me. The disk contained a 323Kb self-extracting ZIP file which, when ran, produced (amongst other files) FRACTINT.EXE which was a 384Kb file. This would not have been a problem except that the machine only had 360Kb disks - FRACTINT.EXE wouldn't fit!

Somewhere (probably on another magazine disk) I found a program that was able to format disks to non-standard formats by calling BIOS functions (or possibly by programming the drive controller directly, I don't remember). It turns out that the drives in the PC1512 were actually able to format 42 tracks before the head bumped into the stop, not just the standard 40. This gave me a 378Kb disk, still not quite enough. I discovered that if I requested the program to format a 44-track disk, it would fail but record in the boot sector that the disk was 396Kb. With this format, I was able to persuade the ZIP program to try extracting FRACTINT.EXE. Of course, it failed to write the last few Kb but it turned out that because Fractint used overlays, it still worked - only one or two of the overlays were inaccessible (I think printing and shelling to DOS were broken, and I didn't miss them).

I made a Mandelbrot zoom video with this setup. It was on the point +i (zooming into a more detailed area would have taken too much time and disk space). I played the video back by converting the images to PCX files (which are much quicker to decode than GIFs) copying them to a ramdrive and then using PICEM to play them back at about 5fps.