GCD

On comp.lang.asm.x86, Jon Kirwan asked for compiler output for the "greatest common divisor" function from the following C implementation:

Here's my favorite C compiler with /os/s options (optimize for size, use no stack checking; other optimization options did not help much) versus a human (i.e., me) implementation:

A particular weakness of this compiler is that it is not sure what registers to start with. But regardless it is missing a huge number of optimization opportunities when compared with the human implementation. There is no comparison is speed or size. Of course, as a human, I have an immense advantage since I can see mathematical relationships that are missed by all C compilers.

Although xchg is not a particularly fast Intel instruction, it does help make it very tight, and probably not more than a cycle off of optimal for performance. The main loop of the routine exists entirely within an instruction pre-fetch buffer (16 bytes.)

Much as I would like, it is pointless to try and describe the exact performance of the algorithms given above. They are performance limited almost entirely by branch mis-prediction penalties and branch taken clocks.

Update: With the advent of modern CPUs such as the Alpha, Pentium-II or Athlon which have deep pipelines and correspondingly bad branch misprediction penalties, it makes sense to remove branches through "predication". In this case we can use min() and max() functions to remove the inner most if() statement in the original loop:

But on a 32 bit machine:

So putting in the appropriate substitutions and simplifying as appropriate we arrive at:

Once again, we take the Pepsi challenge and we get:

So I am still ahead, but in theory a good enough compiler could have equaled my results. Unlike the  first loop, the performance of this loop is characterizable.

Update: Michael Abrash's famous book "The Zen of Code Optimization" actually says that the pure subtraction method is too slow (this contradicts actual testing that I performed.) I have also gotten feedback from people insisting that divide is better for certain special cases or larger numbers. And the final icing on the cake is that in "Handbook of Applied Cryptography" (ISBN 0-8493-8523-7) they give a "divide out powers of two" algorithm that appears to be a good compromise (but still suffers from the branch mis-prediction penalties.) I will need to investigate this a little deeper, and may need to revise these results. Outside analysis is welcome.

Update: Just as a question of completeness, I am presenting what I have discovered so far. Using what apparently is a classical approach, rather than either dividing or subtracting, we subtract powers of 2 times the smaller value from the larger. Furthermore, rather than pursuing the gcd algorithm all the way, once we have small enough values for a and b, we just perform a table lookup (table filling not shown here.)

Obviously some cases of (a,b) input need to be checked (a=0, b=0 sends the program into an infinite loop.)

I believe that it might be possible to improve upon this by performing table look ups on several bits of each operand at once which could give prime linear factors which could be used to reduce the operands by several steps in one shot. So I am refraining from performing assembly language analysis until I resolve this.

Update: There have been a number of submissions that have come to my attention. The first is from Pat D:

The second submission is from Ville Miettinen who works (or worked) for hybrid graphics using the cmovCC instructions:

Scan array of packed structure for non-zero key

Suppose we want to scan an array of structures just to determine if one has a non-zero value as a particular entry. That is to say, suppose we wish to implement the following function:

As it stands, this code suffers from too many problems to analyze at the assembly level immediately. Applying some ingenuity and some standard tricks as described on my optimization page I arrive at:

While it may look a lot longer and seem to add unnecessary complexity, it in fact helps the compiler by practically giving it a variable to register map. Now many may comment that using goto's or labels is bad and causes compilers to turn off optimizations or whatever. Well, my compiler doesn't seem to have a problem with it, which is why I've used it.

Assuming we have:

struct foobar {
    unsigned int key;
    void * content;
};

here's the assembly output of the compiler versus my hand coding:

The compiler wins the size contest, but mine is a heck of a lot faster because I avoid an extra test and branch in the inner loop. As the comments show, my loop saves a clock (on the Pentium.)

Impressed? Well, don't be too impressed. In a sense all of this should be academic. The real problem here is that we are accessing spaced out elements out of an array of packed structures. If the array is substantial in size then we end up losing for all the extra data loaded (the cache can only load contiguous lines of data.) Splitting the structure up into parallel arrays of data which have locality of reference (i.e., all elements are read/written in the relevant inner loops) will tend to do more for performance than these sorts of optimizations. But even in that situation I don't think most compilers will use a repz scasd.

Nevertheless, in situations where you are adhering to a bad array structure design (like the vertex buffer array used in Microsoft's Direct 3D API) the above code does apply.

63 bit LFSR

The following is a typical cyclic redundancy checkcomputing algorithm. It is basically a deterministic iterating step of a function that takes 63 bits of input scrambles the bits in some uniformly pseudo-random way, then returns those 63 bits. Such algorithms can be used for high performance (but low security) encryption, pseudo-random numbers, testing for data integrity or a hash function.

I have implemented this by using 63 of 64 bits in a pair of 32 bit call by reference variables. (My compiler does not support a native 64 bit integer type.)

Here the cycle count is a close one on the Pentium MMX (*in fact pre-MMX CPUs brings the two to a tie since the first shld takes an extra clock to decode.) This is mostly because of the surprisingly slow shld instructions. On Pentium Pro/II, K6 and 6x86MX processors, shld is significantly improved in performance making this a more compelling solution based on size and speed.

The following is optimized purely for performance on Pentium based CPUs.

Quick and dirty bubble sort

Andrew Howe, from Core Designs (makers of Tomb Raider), sent me an extremely short sort loop. Originally written for WATCOM C/C++, I have stripped off the cruft so that you can see it here in its most optimal form.

Notice the use of xchg, stosd, and loop, something you just wont see from C compilers. I did not bother working out the timings of this routine. Remember that for some applications, bubble sort will out perform some of the more sophisticated sorting algorithms. For example, 3D applications that need Z-sorting, but can rely on temporal locality to always maintain at least a mostly sorted list. (For another interesting departure along these lines see <<Sorting with less than all the facts>>)

A fast implementation of strlen()

Recently, someone wrote to me with the comment that strlen() is a very commonly called function, and as such was interested in possible performance improvements for it. At first, without thinking too hard about it, I didn't see how there was any opportunity to fundamentally improve the algorithm. I was right, but as far as low level algorithmic scrutiny is concerned, there is plenty of opportunity. Basically, the algorithm is byte scan based, and as such the typical thing that the C version will do wrong is miss the opportunity to reduce load redundancy.

Here, I've sacrificed size for performance, by essentially unrolling the loop 4 times. If the input strings are fairly long (which is when performance will matter) on a Pentium, the asm code will execute at a rate of 1.5 clocks per byte, while the C compiler takes 3 clocks per byte. If the strings are not long enough, branch mispredictions may make this solution worse than the straight forward one.

Update!

While discussing sprite data copying (see next example) I realized that there is a significant improvement for 32-bit x86's that have slow branching (P-IIs and Athlon.)

I think this code will perform better in general for all 32 bit x86s due to less branching. 16 bit x86's obviously can use a similar idea, but it should be pretty clear that it would be at least twice as slow. (I'm really starting to like this bit mask trick! See next example) Thanks to Zi Bin Yang for pointing out an error in my code.

Here's the C version of the same thing:

It compiles to substantially the same thing (on a P6 it should have the same performance) but it has a few problems with it. Its not really portable to different sizes of int (the numbers 0x7f7f7f7f and 0x80808080 would have to be truncated or expanded, and casting p from a pointer to a scalar is not guaranteed to work) and its not any easier to understand than the assembly.

Note that although alignment on reads is not necessary on modern x86's it is performed here, as a means of avoiding a memory read failure. This can happen as a result of reading past the end of the string, but under the assumption that your memory allocator aligns to at least DWORD boundaries.

Update!

Ryan Mack has sent in an MMX version of this loop:

The individual loop iterations are not dependent on each other, so on a modern out of order processor this is just instruction issue limited. On a P-II we have 6 ALU microops + 1 load microop, which should translate to a 3 clock per 8 byte throughput. On AMD processors there are 7 riscops which is 3 clocks on an Athlon and 4 clocks on a K6 per 8 byte throughput.

Performance measurements on an Athlon show that (ignoring branch   misprediction) for long strings this 64 bit method is upwards of 55% faster than the previous 32 bit method (which in turn is 22% faster than the WATCOM clib implementation.) However, for very short strings, this method is as much as 3 times slower then the previous method (which is about the same same speed as as the WATCOM clib implementation.)  But   keep in mind that if your strings are shorter, the branch mispredicts (which were not modeled in my tests) will take a comparatively larger percentage of the time.  In any event, if you intend to use one of these alternatives, you should performance measure them to see which is better for your application.

Update!

Norbert Juffa has written in to inform me of an old trick which is a significant improvement over the device I discovered:

I did a little digging and indeed you can find an old thread discussing this on google groups. So here's an implementation in C:

Compared to other solutions on an Athlon, this is indeed the fastest solution for strings up to about 32 characters (the MMX solution is faster for longer strings). This function greatly outperforms the strlen of most libraries.

Sprite data copying

A common question that comes up for programmers that are interested in rendering sprite graphics (i.e., game programmers) is how to write a high performance sprite blitter. In a recent discussion on rec.games.programmer, the questioner had a requirement that an array of pixels from one array be copied to a destination array unless the pixel value was zero (representing a transparency.) The starting code looked something like:

Of course this code has several problems with it: (1) It has a branch prediction in the inner loop that is likely to fail very often and (2) The memory is being read one byte at a time which wastes memory bandwidth and increases loop overhead.

That said, this solution is actually good if the memory unit has "write combine" and "write mask" capabilities and that memory write access to the destination is more of a bottleneck than branch prediction, loop overhead or source reads. (This is true if the destination is VRAM, or the CPU is a highly clocked P-II.)

However, if reading the destination is not a problem (if the destination is a system memory buffer), then reading 4 bytes at a time and building a mask is not a bad solution, provided you can build one without using conditional branches (otherwise there would be nothing to gain.) Russ Williams came up with a straight forward approach using the carry flag (it is optimized for the Pentium):

The above loop works great for Pentiums and is out of reach for C compilers because of the use of the carry flag. However the P-II does not like the partial register stall (a non-overlappable 7 clock penalty.) After some discussion motivated by a comment by Christer Ericson, the following trick was observed:

The above expression will generate a bitmask in the high bits of each byte; 00h indicating a zero byte, and 80h indicating a non-zero byte. With some straight forward trickery you can turn this into 00h for zero and FFh for non-zero, giving us the mask we desire. Here's the code in assembly language:

Of course the above loop is not especially suited for Pentiums since practically the whole loop has data dependencies, and there is an AGI stall. These issues can be dealt with by unrolling the loop once and using ordinary addition instead of lea (The edx and ebp registers are available for this).

However, P-II's will attempt to "naturally" unroll this loop, and doesn't suffer from AGI stalls and therefore would stand to benefit much less from hand scheduling of this loop (though it probably would not hurt.) However, since the loop is 22 micro-ops (which is 2 more than the P-II reorder buffer limit) it cannot totally unroll it.

A K6 can also only partially "unroll" this loop because it contains more than 6 decode clocks.

Actually, with the right source, C compilers should have little difficulty coming up with the source to the above (or something very similar.) But I'm not sure how you are supposed to convince your C compiler that it needs to be unrolled once.

Update:After some thought, its obvious that the above algorithm can and should be performed using MMX and some of the new SSE instructions (movntq in particular). I may implement this at a later date.

Another idea, which apparently comes from Terje Mathisen, is to use the P-II's write combining features in conjunction with minimizing branch misprediction penalties. The idea is to use a mask to select between two destination addresses for the write bytes:

Of course, temp is a reused system array meant just as a trash bin for transparency writes. The outer loop gets complicated, but in ideal situations, this will make better use of the P-II chipset features without suffering from branch mispredict penalties. The carry flag trick still eludes modern C compilers. (Of course this loop should be unrolled a bit to improve performance.)

However, like the key scan problem listed above, all this may be academic. For most sprites that are rendered once and used over and over again, transparencies and solid source pixels are usually in long runs. So using an alternate encoding which give transparency versus non-transparency runs will usually be more efficient. Blitting transparency would come down to adding a constant to the pointers, and the blitting solid colors would be immediate without required examination of the source.

Update: Paolo Bonzini has written in with a more efficient inner loop. The idea is apply the mask to the xor difference of the source and destination.

Paolo also submitted a trivial MMX loop:

Now its time for MMX

On sandpile DS asked how to perform the following:

I thought of two ways:

or

The first has fewer instructions but higher latency which is ideal for a  CPU like the Athlon.  The second has more instructions but shorter latency,  which is beneficial to CPUs like the K6, Cyrix or P55C.  I am not sure which is the better loop for P6 architectures.

Of course once we enter the realm of MMX, the poor C compiler really has nothing to offer.

Generalized shift

Norbert Juffa was studying the Compaq Visual Fortran compiler and thought that generalized shift could be improved. In the Fortran language a general shift is always a left shift, but it may take a negative shift value (which is a right shift) and shift values beyond the word size (in either direction) in bits results in 0. Here's the algorithm I came up with:

The ANSI C standard, requires the we use unsigned int when using the right shift operation due to an unspecified "implementation defined" status of right shifting on signed numbers. If you are wondering why I mention it -- its because this is not well known, and can lead to incorrect resuts. This time I have included Microsoft Visual C++:

The comments indicate the clock on which an ideal 3-way out-of-order x86 machine (i.e., an Athlon) could execute these instructions. Notice the issue stalls for each of these code sequences. This means that there may be a case for x86 CPUs to go with 4-way instruction decoding and issue in the future. In any event, MSVC++ seems to be doing quite well in terms of the final performance, but it is unable to simplify the back to back neg reg0 instructions. The classic software register renaming technique for avoiding a dependency stall (after mov esi, eax the role of the two registers can and should be swapped -- the two registers will be identical after this instruction except that one of them is available one clock later than the other) is also not being used. WATCOM C/C++ is having difficulty with register size mismatching, and is unaware of the standard sbb reb, reg trick.

64bit BCD add

Norbert Juffa has given me some devilishly clever code for performing a 64 bit BCD addition.

The use of carries and rcr instructions opposite their C equivalent (sort of) shows the massive simplifications possible in assembly that are not offered in C.

You can find more on BCD arithmetic here at Douglas W. Jones' website

Pixel bashing

"Gordy" has been working on a software based graphics library and asked on USENET for optimized paths for many common pixel bashing routines. Here are some of the solutions that I presented.

The first is an MMX routine for compressing 4 bit per pixel to 1 bit per pixel in a specific ordering. The pixel translation simply takes the bottom bit from each nibble. The order of the original nibbles in a given qword is: 21436587 (don't ask me -- its Gordy's library.) We can assume that the nibble bits are premasked (i.e., the other bits are 0.)

MagicConst dd 48124812h,48124812h ;....     movq     mm7, MagicConst     pcmpeqb  mm6, mm6        ; 1111111111111111     psllw    mm6, 12         ; 1111000000000000 @quads:     movq     mm0, [esi]      ; 0006000500080007     movq     mm1, mm6        ; XXXX000000000000     pmullw   mm0, mm7        ; 8765705800870070     pand     mm1, mm0        ; 8765000000000000     psrld    mm0, 28         ; 0000000000004321     psrlw    mm1, 8          ; 0000000087650000     por      mm0, mm1        ; 0000000087654321     packuswb mm0, mm0        ; [?B?A?B?A]     packuswb mm0, mm0        ; [BABABABA]     movd     eax, mm0     mov      [edi],ax     add      esi, 8     add      edi, 2     dec      ecx     jnz      @quads

Since the bits are isolated, a sequence of shift and ors can be replaced by a sequence of shifts and adds. But a sequence of shifts and adds can be made into a single multiply. From a performance point of view this is a good trade off since pmullw is a fully pipelined instruction and in a modern x86, many internal instances of this loop should be able to execute in parallel. This solution consumes 64 source bytes in just a handful of clocks -- its not likely that any C based formulation could even come close.

The second was for averaging streams of bytes together (and rounding down.) Gordy was looking for a solution for pre-K6-2 MMX/3DNow! (which did not have pavgusb or pavb instructions) as well as a straight integer solution. Here are the solutions I proposed:

; MMX Solution by Paul Hsieh LTop: movq    mm0, [esi] movq    mm1, [edx] movq    mm2, mm0 pand    mm0, mm1 pxor    mm1, mm2 psrlb   mm1, 1 paddb   mm0, mm1 movq    [edi], mm0 dec     ecx jnz     LTop

These solutions use the follow observation:

While the original (A + B) addition may overflow, the addition (A and B) + (A xor B)/2 will not.

The third routine was one for performing a saturating add. Using MMX its just a trivial application of the paddusb instruction (as "Gordy" pointed out), but what of the integer fallback routine? I came up with the following integer routine:

This time we are seperating the bottom bit to make sure we can perform parallel adds with only 8 bits of work area at a time.

The carry over into the 9th bit is translated into a 0x00 or 0x7F mask. By or-ing in the mask into the destination then shifted left by one, this achieves the same effect as if the mask contained 0x00 or 0xFF values.

This time, since the carry flag is captured and used on the first add, the C compiler cannot duplicate this routine.

Update: Oliver Weichhold posted an interesting question in comp.lang.asm.x86:

While one might instictively look at clever ways of performing the multiple shifts and masking that might seem to be required, its really not neccesary. Consider that what you really want is a mapping of 16 input bits that you just want placed in some output positions, where each input bit is independent from the other bits.

This just screams "Look up tables" and indeed it leads to some pretty simple code:

The tables can be set up with some kind of simplistic C code at initialization time.

Converting Uppercase

Conventional wisdom say that to quickly convert an ASCII string to uppercase, that one should use a table to convert each character. The problem is that the address mode of x86 processors require 32 (or 16) bit registers, so some kind of byte -> dword (or word) conversion also happens. This leads to partial register conversion which is penalized on AMD CPUs and is devastating to Intel CPUs.

So what is our alternative? Well, if we can amortize our performance using SIMD, then what we are really trying to do is translate the range [61h, 7Ah] to [41h, 5Ah]. The approach I came up with was to exploit the fact that ASCII is only well defined in the range [00h, 7Fh] by rotating the range then masking, then rotating and capturing the middle range into the 5th bit and subtracting that from the 4 input bytes:

This routine is within reach of most modern C compilers. (Stay tuned for the C version.) A cursory glance, however, reveals that this routine is easily converted to MMX.

Update: The need to be able to support non-ASCII (i.e., include the range [80h, FFh]) has come up. One of the bonuses of being able to do the whole range is that UTF-8 encoded unicode can be directly translated this way. Anyhow, its just a minor change, that unfortunately costs a cycle on the critical path:

Update: Norbert Juffa wrote in to inform me that as part of Sun's opening up of the Solaris source code they have made some of his sources used for some C string functions, including caseless comparisons and tolower available to the public. You can view the source here.

Update: I have designed an alternative function that deals with the entire UTF-8 range, but is also more parallel and thus should execute somewhat faster:

; Integer SIMD uppercase(string) on UTF-8 data v2 by Paul Hsieh mov  eax, src[esi] mov  ebx, 0x80808080 mov  edx, eax          ; src or   ebx, eax          ; (0x80808080 | src) mov  ecx, ebx sub  ebx, 0x7b7b7b7b   ; (0x80808080 | src) - 0x7b7b7b7b sub  ecx, 0x61616161   ; c = (0x80808080 | src) - 0x61616161 not  ebx               ; d = ~((0x80808080 | src) - 0x7b7b7b7b) not  edx               ; ~src and  ebx, ecx          ; c & d and  edx, 0x80808080   ; ~src & 0x80808080 and  ebx, edx          ; e = (c & d) & (~src & 0x80808080) shr  ebx, 2            ; e >> 2 sub  eax, ebx          ; src - (e >> 2);

The C code is fairly comparable:

; Integer SIMD uppercase(string) on UTF-8 data v2 by Paul Hsieh uint32_t upperCaseSIMD (uint32_t x) { uint32_t b = 0x80808080ul | x; uint32_t c =   b - 0x61616161ul; uint32_t d = ~(b - 0x7b7b7b7bul); uint32_t e = (c & d) & (~x & 0x80808080ul);         return x - (e >> 2); }