Writing systems

Over the past few days my background process has been drawing writing systems to my sketchbook from Wikipedia. Got Phoenician, Aramaic, Greek, Etruscan, Latin, Futharks, Hebrew, Arabic, Urdu, Old South Arabian, Ge'ez, Japanese kana, Hangul, Brahmi, Devanagari, Gujarati. Planning to do the uni/bi/triliteral Egyptian hieroglyphs & Demotic & Hieratic, Gothic, Syrian, Meroitic, Ogham, Manchu, Georgian, Armenian, Glacolitic. If I go nuts in the head, maybe 1000 hanzi/simplified & Egyptian hieroglyphs.

Today's amusement

Cristão - kirishitan - 吉利支丹 - 切支丹 - 切死丹 - 鬼理死丹 -- 鬼理死魔巣. 煮翻...

A new set of magnitude numbers

An erratic thought for the morning.

Following the SI-prefixes:
M-illion = 1e6
G-illion = 1e9
Trillion = 1e12
P-illion = 1e15
Exillion = 1e18
Z-illion = 1e21
Yrillion = 1e24

Or maybe the powers of thousand:
Twollion = 1000^2
Threellion = 1000^3
Fourllion = 1000^4
Fivellion (pente-peta) = 1000^5
Sixllion (hexa-exa) = 1000^6
Sevenllion = 1000^7
Eightllion =1000^8

Or the powers of ten?
Sixion = 1e6
Nineion = 1e9
Twelveion = 1e12
Fifteenion = 1e15
Eighteenion = 1e18
Twentyoneion = 1e21
Twentyfourion = 1e24

Would make operating with the mid-powers easier. Sevenion divided by sixion is an onion. Pffft.

Could extend it downwards to negative ions or nins.
Threenin = 1e-3, etc.

Sevenion divided by twelveion is a fivenin. Fivenin times threenin is eightnin. Sevenion times threenin is fourion. (E-notation is better for math though, 1e7 / 1e12 = 1e-5; 1e-5 x 1e-3 = 1e-8; 1e7 x 1e-3 = 1e4.)

Memory aid for multiplying and dividing powers: rotate the operator by 45 degrees clockwise.

Tomtegebra concept art

Concept art from last summer for Tomtegebra, see also this for more background. I ended up scaling it back a little. And by a little I mean by everything.

Gazing at the CPU-GPU crystal ball

Reading AnandTech's Core i7 980X review got me thinking. CPU single-thread performance has roughly doubled over the past four years. And we have six cores instead of just two, for a total speedup in the 5-7x range. In the last two years, GPU performance has quadrupled.

The current top-of-the-line CPU (Core i7 980X) does around 100 GFLOPS at double-precision. That's for parallelized and vectorized code, mind you. Single-threaded scalar code fares far worse. Now, even the 100 GFLOPS number is close to a rounding error compared to today's top-of-the-line GPU (Radeon HD 5970) with its 928 GFLOPS at double-precision and 4640 GFLOPS at single-precision. Comparing GFLOPS per dollar, the Core i7 980X costs $999 and gets roughly 0.1 GFLOPS/$, whereas the HD 5970 costs $599 and gets 1.5 GFLOPS/$ at double precision and 7.7 GFLOPS/$ at single precision.

The GFLOPS numbers are a bit quirky. The HD 5970 number is based on 160 processors running at 725 MHz where each processor can execute five instructions per cycle. And if each is executing a four-element multiply-add, that's 8 FLOPS per instruction. To put it all together, 160 processors * 0.725 GHz * 5 instructions * 8 FLOPS = 4640 GFLOPS. For doubles, each processor can do four scalar multiply-adds per cycle: 160 * 0.725 * 4 * 2 = 928 GFLOPS. If you're not doing multiply-adds, halve the numbers.

The Core i7 GFLOPS number is apparently based on 5 ops per cycle, as 6 cores * 3.33 GHz * 5 FLOPS = 100 GFLOPS. I don't know how you achieve five double ops per cycle. To get four double ops per cycle, you can issue a 2-wide SSE mul and add and they are executed in parallel. Maybe there's a third scalar op executed in parallel with that? If that's the case, you could maybe get 180 GFLOPS with floats: two 4-wide SSE ops and a scalar op for 9 FLOPS per cycle. The GFLOPS/$ for 180 GFLOPS would be 0.18. For a non-multiply-add workload, the numbers are halved to 40 GFLOPS for doubles and 80 GFLOPS for floats.

If we look at normal software (i.e. single-threaded, not vectorized, gets no ops in parallel), the Core i7 does 3.33 GHz * 1 ops = 3.3 GFLOPS. That's a good 30x worse than peak performance, so you better optimize your code. If you're silly enough to run single-threaded scalar programs on the GPU, the HD 5970 would do 0.725 GHz * 2 op mul-add = 1.45 GFLOPS. Again, halve the number for non-mul-add workloads.

Anyhow, looking at number-crunching price-performance, the HD 5970 is 15x better value for doubles and 43x better value for floats compared to the 100 GFLOPS and 180 GFLOPS numbers. If you want dramatic performance numbers to wow your boss with, port some single-threaded non-vectorized 3D math to the GPU: the difference in speed should be around 700x. If you've also strategically written the code in, say, Ruby, a performance boost of four orders of magnitude is not a dream!

Add in the performance growth numbers from the top and you arrive at 1.6x yearly growth for CPUs and 2x yearly growth for GPUs. Also consider that Nvidia's Fermi architecture is reducing the cost of double math to a CPU-style halving of performance. Assuming that these growth trends hold for the next two years and that GPUs move to CPU-style double performance, you'll be seeing 250 GFLOPS CPUs going against 9200 GFLOPS GPUs. The three and four year extrapolations are 409/18600 and 655/37100, respectively. The GPU/CPU performance ratios would be 37x, 45x and 57x for the two-to-four-year scenarios. The corresponding price-performance ratios would be 62x, 75x and 95x.

With regard to performance-per-watt, the Core i7 980x uses 100W under load, compared to the 300W load consumption of the HD 5970. The 980x gets 1 GFLOPS/W for doubles and 1.8 GFLOPS/W for floats. The HD 5970 does 3.1 GFLOPS/W for doubles and 15.5 GFLOPS/W for floats.

If number crunching performance was all that mattered in a CPU price, top-of-the-line CPUs would be priced at $50 today. And $10 in four years...

The CPU-GPU performance difference creates an interesting dynamic. For Intel, the horror story is there being no perceivable difference between a $20 CPU with a $100 GPU vs. a $500 CPU with a $100 GPU. That would blow away the discrete CPU market and you'd end up with cheap CPUs integrated on the motherboard with the GPU as the differentiating factor. Much like the Atom boxes, come think of it. The plain Atom can't play video or play 3D games. An Atom with ION can [edit] play videos that have GPU-accelerated codecs and light 3D games [/edit].

Intel strategy?

To slow down the migration towards GPUs, you need to make targeting GPUs an unattractive proposition to software developers. A software developer does a cost-benefit analysis based on the size of the market and the cost of entering it. To keep them away from a market, reduce the its size and increase the cost. When it comes to retarding GPU computing, the goal is to make fewer computers ship with a capable GPU and to make developing for GPUs harder. If you're a GPU manufacturer, your goals are the exact opposite.

To make fewer computers ship with a capable GPU, they must ship with a graphics-only GPU. Competition in the GPU market is based on performance, so vendors are unlikely to buy a slow GPU. To make vendors buy a slow GPU, you must integrate it to the package you're selling, and thus lower the probability that they will buy a superfluous extra GPU. You also want to make it harder for vendors to buy a motherboard with a capable GPU integrated. These actions should restrict the market of capable GPUs to enthusiast machines.

This is still only a stop-gap measure, as it's pretty difficult to beat a 50x performance gap. In the long term, you want to switch the graphics-only integrated GPU to a home-built capable GPU and then proceed with commoditizing the CPU market. Or if that fails, maybe go for a monopoly on the chipset market: use IP to keep out competition, hoist prices.

Of course, this kind of maneuvering tends to lead to lawsuits.

4x4 float matrix multiplication using SSE intrinsics

Here's a bit of fun from earlier this week: 4x4 matrix multiplication using SSE. Benchmarked it to be 4x faster than the scalar version (on a Pentium M, using GCC 4.4.1 with flags -O3 -mfpmath=sse -msse2 -std=c99).

To refresh your memory, here's the definition of matrix multiplication: A x B = R, where R[i][j] = dot(row(A, i), column(B, j)) and dot(u, v) = sum(i=0 to N, u[i] * v[i]). Only defined when A's width equals B's height. Matrix multiplication is used for transformations in computer graphics, 4x4 matrix multiplication specifically in 3D graphics.

First the scalar version. Multiply A with B and store the result in R. The matrices are stored in column-major order in a flat float array (so that the 4x4 transform matrix translation column is the last four values in the array).

void 
mmul(const float *a, const float *b, float *r)
{
  for (int i=0; i<16; i+=4)
    for (int j=0; j<4; j++)
      r[i+j] = b[i]*a[j] + b[i+1]*a[j+4] + b[i+2]*a[j+8] + b[i+3]*a[j+12];
}

Now, a naive slow SSE version. Load a row from A and a column from B, multiply them together and add the values of the result vector together.

Note that we need to align all the matrices to 16-byte boundaries. You can use memalign(16, 16*sizeof(float)) or posix_memalign or compiler pragmas to do that. If the matrices are not aligned to 16-byte boundaries, the SSE code will segfault. I'm using the SSE intrinsics, so you won't see any cool inline assembler tricks here.

void
mmul_sse_dot(const float * a, const float * b, float * r)
{
  __m128 a_line, b_line, r_line;
  float mc[16] __attribute__((aligned(16)));  // 16-byte aligned temp array
  for (int i=0; i<16; i+=4) {
    b_line = _mm_load_ps(&b[i]);              // b_line = vec4(column(b, i))
                                              // remember that we are column-major
    for (int j=0; j<4; j++) {
      a_line = _mm_set_ps(a[j+12], a[j+8], a[j+4], a[j]); 
                                              // a_line = vec4(row(a, j))
                                              // note that SSE is little-endian
      r_line = _mm_mul_ps(a_line, b_line);    // r_line = a_line * b_line
      _mm_store_ps(mc, r_line);               // copy r_line to memory
      r[i+j] = mc[0]+mc[1]+mc[2]+mc[3];       // r[i][j] = sum(r_line)
                                              //         = dot(a_line, b_line)
    }
  }
}

Sadly, this way is slow in practice (30% slower than the scalar version, in fact) as it stores only one value at a time and there's no operation to sum the elements of an SSE vector.

What we must do instead is execute four dots in parallel. Multiply the first value of the row A[I] with the first row of B and store the result to row R[I] (R[I] = A[I][0] * B[0]). Then multiply the second value of row A[I] with the second row of B, add to row R[I] (R[I] += A[I][1] * B[1]). Repeat for all the remaining values in row A[I].

Now the first value of row R[I] is the sum of products of the row A[I] and the first values of the rows of B. I.e. the dot product of row A[I] and the first column of B.

Correspondingly, the second value of R[I] is the sum of products of row A[I] and the second values of the rows of B. I.e. the dot of row A[I] and the second column of B. And so forth.

In pseudo-code:

R[I] = A[I][0] * B[0] + 
       A[I][1] * B[1] + 
       A[I][2] * B[2] + 
       A[I][3] * B[3]

But as our matrices are column-major, we need to flip A and B, so that we're multiplying whole columns with the values of a row. Ok, time for code:

void
mmul_sse(const float * a, const float * b, float * r)
{
  __m128 a_line, b_line, r_line;
  for (int i=0; i<16; i+=4) {
    // unroll the first step of the loop to avoid having to initialize r_line to zero
    a_line = _mm_load_ps(a);         // a_line = vec4(column(a, 0))
    b_line = _mm_set1_ps(b[i]);      // b_line = vec4(b[i][0])
    r_line = _mm_mul_ps(a_line, b_line); // r_line = a_line * b_line
    for (int j=1; j<4; j++) {
      a_line = _mm_load_ps(&a[j*4]); // a_line = vec4(column(a, j))
      b_line = _mm_set1_ps(b[i+j]);  // b_line = vec4(b[i][j])
                                     // r_line += a_line * b_line
      r_line = _mm_add_ps(_mm_mul_ps(a_line, b_line), r_line);
    }
    _mm_store_ps(&r[i], r_line);     // r[i] = r_line
  }
}

Now we're cooking. 4x faster than the scalar version. And it even produces the same results.

You can get the full code here: mmul.c.

[Edit]
Wrote a C++ version with a vec4 struct to make the code super short and simple:

void mmul_vec4
(const float * a, const float * b, float * r)
{
  for (int i=0; i<16; i+=4) {
    vec4 rl = vec4(a) * vec4(b[i]);
    for (int j=1; j<4; j++)
      rl += vec4(&a[j*4]) * vec4(b[i+j]);
    rl >> &r[i];
  }
}

The best thing? Vec4 is actually a wrapper on __m128 SSE operations, so mmul_vec4 runs as fast as mmul_sse. And the wrapper was dirt simple to write too, see below.

struct vec4
{
  __m128 xmm;

  vec4 (__m128 v) : xmm (v) {}

  vec4 (float v) { xmm = _mm_set1_ps(v); }

  vec4 (float x, float y, float z, float w)
  { xmm = _mm_set_ps(w,z,y,x); }

  vec4 (const float *v) { xmm = _mm_load_ps(v); }

  vec4 operator* (const vec4 &v) const
  { return vec4(_mm_mul_ps(xmm, v.xmm)); }

  vec4 operator+ (const vec4 &v) const
  { return vec4(_mm_add_ps(xmm, v.xmm)); }

  vec4 operator- (const vec4 &v) const
  { return vec4(_mm_sub_ps(xmm, v.xmm)); }

  vec4 operator/ (const vec4 &v) const
  { return vec4(_mm_div_ps(xmm, v.xmm)); }

  void operator*= (const vec4 &v)
  { xmm = _mm_mul_ps(xmm, v.xmm); }

  void operator+= (const vec4 &v)
  { xmm = _mm_add_ps(xmm, v.xmm); }

  void operator-= (const vec4 &v)
  { xmm = _mm_sub_ps(xmm, v.xmm); }

  void operator/= (const vec4 &v)
  { xmm = _mm_div_ps(xmm, v.xmm); }

  void operator>> (float *v)
  { _mm_store_ps(v, xmm); }

};

Here's the full code for the C++ version: mmul.cpp.

3D models and parsing binary data with JavaScript

Background

A few days ago, I was fiddling around with loading 3D models for use with WebGL (here demos). I had an OBJ parser and some test models lying around from an earlier project, and the parser was only a two hundred lines, so I ported it over to JavaScript. But the OBJ files were pretty huge, a 40-thousand quad model without normals or texcoords was 2.4 megs. Not really something you want to send over the internet.

The OBJ file format stores coordinates and face indexes as text. While it makes the model files human-readable and easy to parse, it also increases the file size. Uncompressed, a single coord takes 9 bytes, plus an optional sign byte.

My first instinct was to just use gzip compression on server-side. And it did manage to bring the file size down to 880k for the above-mentioned 2.4 meg OBJ. Still. Maybe there would be more gains to be had by storing the coords and face indexes as binary data.

I went for [0,1]-normalized 16-bit fixed-point numbers for the coordinates and 16-bit unsigned ints for the face indexes. I also had to store the vertex and face counts, for which I used 32-bit unsigned ints (which, come think of it, is a bit silly, as using 16-bit face indexes limits the vert count to 65536). Finally, I added two three-float vectors for coord scaling and translation, in order to inflate the normalized coords back to something close to their original values.

Result: 560k for uncompressed binary file, 470k gzipped. Nice. Though it could be brought down further by using a triangle strip instead of separate quads and tris. And there are fancy 3D mesh compression papers that talk about compressing down to less than 2 bits per triangle (for the example model: 40kquads -> 80ktris = 160kb = 20kB). But enough of that for now. Next up, how to parse the binary data.

Parsing binary data

To parse a binary string with JavaScript, you use charCodeAt and the usual bit munging operators (&, |, <<, >>). First, you need the byte value at some point in the string. As charCodeAt operates on characters and not bytes, you need to tell the browser to get the data as an unparsed string by calling overrideMimeType("text/plain; charset=x-user-defined") on your XMLHttpRequest. Then mask away the high byte of each char code to get the low byte value. The following works at least on Firefox and WebKit, I don't know about IE and Opera.

  var byteValue = str.charCodeAt(index) & 0xFF; 

Loading unsigned ints is simple enough, just load the bytes for the int and shift them to their places:

  // read big-endian (network byte order) unsigned 32-bit int from data, at offset
  readUInt32 = function(data, offset) {
    return ((data.charCodeAt(offset) & 0xFF) << 24) +
           ((data.charCodeAt(offset+1) & 0xFF) << 16) +
           ((data.charCodeAt(offset+2) & 0xFF) << 8) +
           (data.charCodeAt(offset+3) & 0xFF);
  }
  // read big-endian (network byte order) unsigned 16-bit int from data, at offset
  readUInt16 = function(data, offset) {
    return ((data.charCodeAt(offset) & 0xFF) << 8) +
           (data.charCodeAt(offset+1) & 0xFF);
  }

Reading in normalized fixed-point numbers as floats isn't much harder either, as they're generated with uint16(normalized_float * 65535):

  readNormalizedUFixedPoint16 = function(data, offset) {
    return readUInt16(data, offset) / 65535.0;
  }

Floats, on the other hand, are a bit more involved. Writing the float parser was educational, now I have a better understanding of the buggers! The following code doesn't handle the special numbers and likely loses precision, is slow, etc., so if you have a better way of parsing binary floats, please share.

  // read big-endian (network byte order) 32-bit float
  readFloat32 = function(data, offset) {
    var b1 = data.charCodeAt(offset) & 0xFF,
        b2 = data.charCodeAt(offset+1) & 0xFF,
        b3 = data.charCodeAt(offset+2) & 0xFF,
        b4 = data.charCodeAt(offset+3) & 0xFF;
    var sign = 1 - (2*(b1 >> 7));                     // sign = bit 0
    var exp = (((b1 << 1) & 0xff) | (b2 >> 7)) - 127; // exponent = bits 1..8
    var sig = ((b2 & 0x7f) << 16) | (b3 << 8) | b4;   // significand = bits 9..31
    if (sig == 0 && exp == -127)
      return 0.0;
    return sign * (1 + sig * Math.pow(2, -23)) * Math.pow(2, exp);
  }

And there you have it. I haven't gotten around to parsing signed integers, though I guess subtracting 1 << (bits-1) from the unsigned value would do the trick.