Fixed point maths in x86 assembly

Fixed point stores a number with a fractional part in an ordinary integer, with the point at a fixed place in it. It's an approximation to floating point with a limited range, and it's typically used for fast 3D calculations and inside polygon shading loops, where the integer instructions are quicker than the FPU, or there isn't an FPU at all.

The usual format on the PC is 32 bits, with a 16-bit integer part and a 16-bit fractional part (fixed 16:16), which fits nicely into a 32-bit register. The code here is 386 assembly, with C where it helps, and I'll give some of the explanations in the more familiar base 10 so that it's easy to see what's going on. Let's start with the format.

The format

A fixed point number is an integer that's really the number multiplied by a constant. In floating point the position of the point can move; here it stays where it is. I'm sticking to powers of two, in particular 2¹⁶ (65536), so 00010000h represents 1.0 and 00038000h represents 3.5, because 8000h is half of 10000h, one:

The decimal equivalent, which I'll call fixed 2:2, takes 100 to mean 1, so 150 is 1.5. To convert any real number to this format you just multiply by 100 and round to the nearest integer.

In C a fixed point number is a 32-bit signed integer. With a 32-bit compiler that's an int or a long; with a 16-bit one it has to be a long. Converting from floating point means multiplying by 65536, and storing the result in an integer drops whatever is left after the point. Dividing by 65536.0 converts it back:

typedef long fixed;                  //  32 bits, 16:16

double r = 3.5;
fixed  f = (fixed)(r * 65536.0);     //  00038000h
double s = f / 65536.0;              //  3.5 again

Fixed or floating point

An important question is whether fixed point has any advantage over floating point for your application. These are the main points for and against each:

Usually a mixture of fixed and floating point turns out to be best. In C, fixed point for simple arithmetic probably isn't worthwhile, since the FPU is pretty fast and the shifts that multiplies and divides need make them not much faster than floating point, and much less versatile. In assembler, though, fixed point has a clear speed advantage, and groups of fixed point instructions can usually be scheduled together quite well. The places where it pays off even in a high-level language are square roots and trig, so those get the most space below.

Adding and subtracting

Fixed point addition is a single 32-bit add instruction, and subtraction a single sub:

_result dd ?

start:
  mov eax, 00010000h    ;  1.0
  add eax, 00024000h    ;  + 2.25
  sub eax, 00058000h    ;  - 5.5
  mov _result, eax      ;  = -2.25

The negative result is taken care of: _result is FFFDC000h, which is -2.25 in two's complement. It's the same in C:

fixed result = 0x10000 + 0x24000 - 0x58000;   //  1.0 + 2.25 - 5.5 = -2.25

These work because for addition and subtraction it doesn't matter where the point is, as long as it's in the same place in both numbers. If the numbers you are adding are 10000h times too large, then the result will be too.

Negation and absolute value

Negation is the normal 32-bit neg instruction (or f = -f in C). The absolute value can be done with a compare and a jump, but it's faster without the jump, using a common trick for things that would otherwise need one:

abs_eax:
  cdq                   ;  EDX = 00000000h if EAX is positive, FFFFFFFFh if negative
  xor eax, edx          ;  Flip every bit if it was negative...
  sub eax, edx          ;  ...and add one, which makes it positive
  ret

cdq copies the sign bit of EAX into every bit of EDX. For -2.25, FFFDC000h, that makes EDX FFFFFFFFh, the xor flips every bit to give 00023FFFh, and subtracting -1 adds one to give 00024000h, 2.25. For a positive number EDX is 0 and neither instruction changes anything. cdq only works on EAX, so for another register you xchg it into EAX first and back again afterwards.

In C there's a version of the same trick:

f = (f ^ (f >> 31)) - (f >> 31);     //  f >> 31 is 0, or -1 if f is negative

It relies on two things: f is a signed integer, and the compiler sign extends when it shifts a signed integer right, which the PC compilers do (the C standard leaves it up to the compiler). The f >> 31 becomes your cdq. The 31 would have to be 15 or 7 for 16 and 8-bit integers. Thanks to Karl for this one :)

Multiplication

When you multiply two fixed point numbers, the result has twice as many fractional digits. In decimal 2:2, 1.50 × 3.20 = 4.80, but 150 × 320 = 48000. That's 4.8 in fixed 4:4, and it needs dividing by 100 to turn it back into 2:2.

The intermediate 4:4 result needs room. 09.99 × 09.99 fits back into 2:2 as 9980, 99.80, but on the way it's 998001 in 4:4, which is twice as many digits. The largest value 2:2 can hold is 99.99, and that can't be multiplied by anything greater than 1.0 without overflowing.

This is where the 386 helps. imul multiplies EAX by a 32-bit number and leaves the whole 64-bit product in EDX:EAX, which for two 16:16 numbers is fixed 32:32: the integer part is in EDX and the fraction is in EAX. The 16:16 answer is the middle 32 bits:

_num1   dd  00018000h   ;  1.5
_num2   dd -00020000h   ;  -2.0
_result dd ?

start:
  mov eax, _num1
  imul _num2            ;  EDX:EAX = the product, in fixed 32:32
  mov ax, dx            ;  Put the integer part under the top of the fraction...
  rol eax, 16           ;  ...and swap the halves, to leave fixed 16:16
  mov _result, eax      ;  = -3.0

Negatives are handled correctly, and _result ends up as FFFD0000h, -3.0. The top half of EDX is the overflow from the multiplication: if you need more range, EDX and the top half of EAX together are a 48-bit result in fixed 32:16. The bottom half of EAX is extra fraction, which is of doubtful accuracy because the numbers going in were only accurate to 16 bits.

It's the intermediate result that causes problems in C. It has to be a type twice the size of the fixed point format, so for 16:16 it has to be 64 bits, but in C the result of a multiply is the same type as the numbers you were multiplying. 0x10000 * 0x10000 should give 0x100000000, which shifted right by 16 bits gives the correct 0x10000, but only the bottom 32 bits are kept, so the answer is 0, which is terribly wrong. In 2:2 terms, it's like working out 09.99 × 09.99 and keeping only the bottom four digits of 998001, 8001, which leaves 00.80 after the divide.

One way round it is to shift each number right by 8 bits before the multiplication. That loses a lot of accuracy, since the fractional part of the result is only accurate to 8 bits, but at least the top of the number survives:

fixed fmul8(fixed a, fixed b)
{
    return (a >> 8) * (b >> 8);      //  Only 8 fractional bits of each are used
}

It's much easier if your compiler has a 64-bit integer type, such as __int64. Be careful though, because some compilers don't do 64-bit arithmetic inline: they call library functions for it, which are much slower.

fixed fmul(fixed a, fixed b)
{
    return (fixed)(((__int64)a * b) >> 16);
}

That's not a very satisfactory choice, so I would suggest that any fixed point maths in C is written in inline assembly.

Division

Division goes the other way. idiv divides the 64-bit number in EDX:EAX by a 32-bit one, so the number to be divided has to be in fixed 32:32 form, like the result straight after a multiply, and the divisor is stored the normal way. In decimal 2:2, 1.5 / 2 = 0.75 is 150 / 200, with the 150 multiplied up by 100 first to make 15000, so that 15000 / 200 = 75, which is 0.75 in 2:2.

_num1   dd  00018000h   ;  1.5
_num2   dd -00020000h   ;  -2.0
_result dd ?

start:
  mov eax, _num1
  rol eax, 16           ;  Fraction to the top half of EAX, integer part to the bottom
  movsx edx, ax         ;  EDX = the integer part, sign extended
  xor ax, ax            ;  EAX = the fraction alone, so EDX:EAX is 1.5 in fixed 32:32
  idiv _num2            ;  EAX = EDX:EAX / -2.0, in fixed 16:16
  mov _result, eax      ;  = -0.75

_result comes out as FFFF4000h, -0.75. The C version needs the 64-bit type again:

fixed fdiv(fixed a, fixed b)
{
    return (fixed)(((__int64)a << 16) / b);
}

This can cause all sorts of undesirable divide errors if it isn't used carefully. 00080000h / 00000001h, 8.0 divided by 1/65536, is 524288, which is too big for 16:16, and the processor stops the program with a divide error rather than give a wrong answer. The same happens when dividing by zero.

A divide that can't overflow

When dividing by something that can get very small, like the depth in a perspective projection, it's handy to have a divide that can't fail. This one returns the largest number of the right sign instead, which works as plus or minus infinity:

;  Inputs:
;    EAX = dividend (fixed 16:16)
;    ECX = divisor (fixed 16:16)
;  Returns:
;    EAX = EAX / ECX (fixed 16:16), or if that won't fit,
;          the largest value of the right sign

fdiv_safe:
  push ebx ecx edx
  mov ebx, eax
  xor ebx, ecx          ;  The sign bit of EBX is now the sign of the answer
  test eax, eax         ;  Make both numbers positive
  jns short fds_1
  neg eax
fds_1:
  test ecx, ecx
  jns short fds_2
  neg ecx
fds_2:
  mov edx, eax          ;  EDX:EAX = the dividend in fixed 32:32
  shr edx, 16
  shl eax, 16
  cmp edx, ecx          ;  The answer only fits in 32 bits if EDX < ECX,
  jae short fds_big     ;  which also catches dividing by zero
  div ecx
  test eax, eax         ;  It has to fit in 31 bits, to leave room for the sign
  js short fds_big
  test ebx, ebx         ;  Put the sign back
  jns short fds_done
  neg eax
fds_done:
  pop edx ecx ebx
  ret
fds_big:
  mov eax, 7FFFFFFFh    ;  The largest positive number...
  test ebx, ebx
  jns short fds_done
  inc eax               ;  ...or 80000000h, the most negative one
  jmp short fds_done

It divides the two sizes as unsigned numbers and puts the sign back afterwards, which makes the overflow test simple. An unsigned div only fails when the answer needs more than 32 bits, and that happens exactly when the top half of the dividend, EDX, is at least the divisor, so one compare catches it, including a divisor of zero. The answer then has to fit in 31 bits as well, since a signed 16:16 number only has 31 bits for its size. 32767.0 / 0.5 gets through the first test but comes out as FFFE0000h, and without the second test that would be returned as -2.0.

Multiply then divide

Working out a × b / c takes advantage of the fact that the 32:32 result from imul can be plugged straight into idiv, with no shifting at all, and the whole 64-bit intermediate keeps every bit of accuracy:

_num1   dd  00018000h   ;  1.5
_num2   dd  00030000h   ;  3.0
_num3   dd  00040000h   ;  4.0
_result dd ?

start:
  mov eax, _num1
  imul _num2            ;  EDX:EAX = 1.5 * 3.0, in fixed 32:32
  idiv _num3            ;  EAX = EDX:EAX / 4.0, in fixed 16:16
  mov _result, eax      ;  = 1.125

As long as b / c is between -1 and 1, the answer is no bigger than a, so there's no risk of a divide error. That's often the case with this calculation: stepping along a polygon edge, for instance, the x at a row is x1 + (x2 - x1) × (y - y1) / (y2 - y1), where the fraction is never more than one.

Sine and cosine tables

I use a sine and cosine lookup table with 1024 entries for 360 degrees, each one a fixed 16:16 value. The cosine is the same curve as the sine a quarter of a turn on, so the cosine table can overlap the sine table quite considerably, saving space: one table of 1024 + 256 entries holds both.

.DATA

_sin    dd  256 dup (?)
_cos    dd 1024 dup (?)

_sin[i] for i up to 1023 reads on into _cos, and _cos[i] is _sin[i + 256], the sine 90 degrees later, which is the cosine.

The table is filled in at the start of the program. The first 512 entries, half a turn, come from the Taylor series for sine, x - x³/3! + x⁵/5! - …, and the rest is filled in by symmetry, since sin(x) = -sin(x - 180°). It's calculated in fixed 8:24, so that the 16:16 result that's stored is accurate to 16 fractional bits:

_sinval   dd ?
_sinrad   dd ?
_sinpower dd ?

calc_sin:
  xor ebx, ebx          ;  EBX = the table entry, 0 to 511
  xor di, di            ;  DI = its offset in the table
  mov ebp, 512

calc_sinloop:

    ;  Convert to radians (= EBX * PI / 512), in fixed 8:24.

  mov ecx, 3243F6Bh     ;  PI * 1000000h
  mov eax, ebx
  mul ecx
  div ebp
  mov _sinrad, eax
  mov _sinpower, eax

    ;  Add up the series x - x^3/3! + x^5/5! - ...

  xor ecx, ecx
  mov _sinval, ecx
  inc cx

calc_sincalc:

  mov eax, _sinpower    ;  Add this term
  add _sinval, eax

 REPT 2
  inc cx                ;  Next term = this one * x / CX
  mov eax, _sinpower
  mul _sinrad
  rol edx, 8            ;  Take the 8:24 part of the 16:48 product
  rol eax, 8
  mov dl, al
  mov eax, edx
  xor edx, edx
  div ecx
  mov _sinpower, eax
 ENDM

  mov eax, _sinpower    ;  Subtract the next term
  sub _sinval, eax

 REPT 2
  inc cx
  mov eax, _sinpower
  mul _sinrad
  rol edx, 8
  rol eax, 8
  mov dl, al
  mov eax, edx
  xor edx, edx
  div ecx
  mov _sinpower, eax
 ENDM

  cmp cx, 256           ;  Plenty of terms, long after they reach zero
  jng calc_sincalc

  mov eax, _sinval      ;  Store it in fixed 16:16
  sar eax, 8
  mov _sin[di], eax
  add di, 4

  inc bx
  cmp bx, 512
  jne calc_sinloop

    ;  Make the rest of the table using sin(x) = -sin(x - 180 degrees).

calc_sinrest:
  mov eax, _sin[di-512*4]
  neg eax
  mov _sin[di], eax
  add di, 4
  cmp di, (1024+256)*4
  jne calc_sinrest

  ret

For each entry the angle in radians is the entry × π / 512, with π as 3243F6Bh, 3.14159265 × 2²⁴. Entry 256 comes out as 01921FB5h, 1.5707963, and its series adds up to exactly 1.0, 00010000h.

The series doesn't need its powers and factorials worked out separately: each term is the last one multiplied by x twice and divided by the next two counts, which is what the two REPT 2 blocks do between an add and a subtract. Multiplying two 8:24 numbers gives a 16:48 product in EDX:EAX, and the two rols and the mov dl, al take the 8:24 part out of the middle of it. The terms are always positive and the add and subtract alternate them, so the running total can dip below zero along the way without any harm.

It's probably not as efficient as it could be, but I don't care: it runs once at the start, and every entry comes out within 1/65536 of the true sine.

Using the tables

I usually store angles as 16 bits, so a value of 10000h is 360 degrees. That's the fractional part of a fixed 16:16 number, the fraction of the way around a complete circle the angle represents, and it means angles wrap round by themselves: adding 90 degrees (4000h) to 300 degrees overflows into 30 degrees, with no checks needed.

The tables are most useful for 3D rotation, where a point is turned by multiplying it by the sine and cosine of the angle. Here is a basic example of how to use them:

;  Inputs:
;    BX  = angle (0-65535 for a full turn)
;    ECX = radius (fixed 16:16)
;  Returns:
;    ESI = ECX * cos(BX)
;    EDI = ECX * sin(BX)
;    BX is changed

rsincos:
  push eax edx

    ;  Turn the angle into an offset in the tables.
  shr bx, 6
  shl bx, 2
    ;  Cosine multiply.
  mov eax, _cos[bx]
  imul ecx
  mov ax, dx
  rol eax, 16
  mov esi, eax
    ;  Sine multiply.
  mov eax, _sin[bx]
  imul ecx
  mov ax, dx
  rol eax, 16
  mov edi, eax

  pop edx eax
  ret

The top 10 bits of the angle pick the entry and the bottom 6 are dropped, so angles go in steps of 0.35 degrees. 30 degrees is 1555h, which is entry 85, 29.88 degrees, and with a radius of 10.0 that gives 8.6703 and 4.9821. When that isn't close enough, the bottom 6 bits can be used to interpolate between one entry and the next.

CORDIC

CORDIC works out the sine and cosine (both at once) with nothing but shifts and adds, to whatever accuracy the numbers can hold, with no table of results at all. It can also find multiples of sines and cosines, the angle and length of a vector, and exponentials.

It starts from the standard two-dimensional rotation matrix, abbreviating cos(θ) to c, sin(θ) to s and tan(θ) to t, and taking the cosine out as a scalar:

[ c  -s ]         [ 1  -t ]
[ s   c ]  =  c * [ t   1 ]

The idea is that a rotation by any angle can be made up from a fixed, predetermined set of smaller and smaller rotations, each one taken either clockwise or anticlockwise, until we reach the right angle. Two choices make that cheap. First, the preset rotations are the angles whose tangent is a power of two, 1, 1/2, 1/4 and so on, which are 45°, 26.57°, 14.04°, 7.13°… Then multiplying by t is just a shift, and each rotation is two shifts, an add and a subtract. Second, the number of rotations is fixed in advance. cos(θ) = cos(-θ), so it doesn't matter which way each rotation goes: all the c factors multiply together into one constant, 0.60725 for sixteen steps.

So to find the sine and cosine of an angle, start with x = 0.60725 and y = 0, which cancels that constant out in advance, and turn towards the angle: anticlockwise while there's angle left to turn, clockwise after overshooting. At the end, x is the cosine and y the sine.

In C each step looks like this, where i is the step and also the shift:

for (i = 0; i < 16; i++) {
    dx = x >> i;
    dy = y >> i;
    if (a >= 0) {                    //  Angle left to turn: anticlockwise
        x -= dy;
        y += dx;
        a -= cordic_angle[i];
    } else {                         //  Overshot: clockwise
        x += dy;
        y -= dx;
        a += cordic_angle[i];
    }
}

Here are the steps for 30 degrees. The angle left to turn swaps sign as it overshoots one way and then the other, and roughly halves each time, so each step gives about one more bit of accuracy:

(The angle is 1555h, which is 29.998 degrees, so the true answers are 0.86604 and 0.49997.) In assembler it's a short loop:

;  Inputs:
;    BX = angle (0-65535 for a full turn)
;  Returns:
;    ESI = cos(BX) (fixed 16:16)
;    EDI = sin(BX) (fixed 16:16)

    ;  atan(2^-i) for each step, in 1/1000000h of a turn
_cordic_angle dd 200000h, 12E405h, 9FB38h, 51112h, 28B0Dh, 145D8h, 0A2F6h, 517Ch
              dd 28BEh, 145Fh, 0A30h, 518h, 28Ch, 146h, 0A3h, 51h

cordic:
  push eax ebx ecx edx ebp

    ;  CORDIC only reaches 99 degrees each way, so turn anything past
    ;  90 degrees round by 180, and negate the answer at the end.

  xor ebp, ebp
  mov ax, bx
  add ax, 4000h         ;  Positive if BX is within 90 degrees of 0
  jns short cordic_ok
  add bx, 8000h
  not ebp               ;  EBP = FFFFFFFFh: negate the answer
cordic_ok:

  movsx ebx, bx         ;  EBX = the angle left to turn, in 1/1000000h of a turn
  shl ebx, 8
  mov esi, 26DD3B6Ah    ;  x = 0.60725 in fixed 2:30, to cancel the growth
  xor edi, edi          ;  y = 0
  xor ecx, ecx          ;  ECX = i

cordic_loop:
  mov eax, edi          ;  EAX = y >> i
  sar eax, cl
  mov edx, esi          ;  EDX = x >> i
  sar edx, cl
  test ebx, ebx
  js short cordic_back
  sub esi, eax          ;  Turn anticlockwise
  add edi, edx
  sub ebx, _cordic_angle[ecx*4]
  jmp short cordic_next
cordic_back:
  add esi, eax          ;  Turn clockwise
  sub edi, edx
  add ebx, _cordic_angle[ecx*4]
cordic_next:
  inc ecx
  cmp ecx, 16
  jne cordic_loop

  sar esi, 14           ;  Fixed 2:30 back to 16:16
  sar edi, 14
  xor esi, ebp          ;  Negate them both if the angle was turned round,
  sub esi, ebp          ;  with the same trick as abs
  xor edi, ebp
  sub edi, ebp

  pop ebp edx ecx ebx eax
  ret

The rotations add up to 99.88 degrees, so on its own CORDIC can only reach angles within about 99 degrees of zero. Anything more than 90 degrees away is turned round by 180 degrees first, by adding 8000h to the angle, and both answers are negated at the end with the same xor and sub trick as the absolute value.

The fixed point formats are chosen for accuracy. x and y never get bigger than 1, so they're kept in fixed 2:30, with 14 more fractional bits than 16:16, to soak up the rounding of all those shifts, and the angle is kept 8 bits finer than the 16-bit angle that comes in. With 16:16 and the plain 16-bit angle the answers are only good to about 20/65536, because the angle runs out of bits long before x and y do; like this, every angle round the circle comes out within 3/65536.

A few notes on the method:

• The table of rotations is the only table there is, one entry per step, and the number of steps sets the accuracy: roughly one bit per entry.
• The loop can be unrolled, which turns every shift into a constant, and on the Pentium some of the shifts and adds can be paired.
• Without the starting 0.60725, CORDIC rotates any vector (x, y) by the angle and multiplies its length by 1.64676, which can be useful by itself.
• Run backwards, it turns (x, y) into a length and an angle: start with the angle at zero and choose each rotation by the sign of y instead, turning the vector down onto the x axis. When y reaches 0, x is the length (times 1.64676) and the angle is how far it turned.
• With the hyperbolic functions cosh and sinh in place of cos and sin (and the rotation matrix [ ch sh ; sh ch ]), it gives exponentials, since exp(x) = cosh(x) + sinh(x), and run backwards it gives logarithms, although the tables needed are larger.

Square roots

There are three ways of working out a square root here: an iterative one using Newton's method, which is the one I use for unit vectors, a bitwise one that only needs shifts and adds, and a recursive one.

Newton's method

Square roots can be calculated by means of an iterative formula, x := (x² + n) / 2x, which can be optimised to x := (x + n/x) / 2, which is efficient: one divide per step. Starting from any guess above the root, each step gets closer, and the number of correct bits roughly doubles each time once it's near. Here it is finding the root of 2.0, from a first guess of just under 2:

The condition for stopping is that the new value isn't any smaller than the old one. With whole numbers, many values toggle between two answers one apart instead of settling on one, which stops a simple "has it changed" test from catching them, but stopping as soon as it stops going down always lands on the lower of the two, which is the root rounded down.

;  Inputs:
;    EDX:EAX = number (fixed 32:32, positive)
;  Returns:
;    EAX = square root (fixed 16:16)

sqrt:
  push ebx ecx edx esi edi
  and edx, 7FFFFFFFh    ;  Clear the sign bit
  mov esi, eax          ;  Keep the number in EDI:ESI
  mov edi, edx

    ;  First guess: 2^k - 1, where 2^k is just above the root

  mov ebx, -1
  bsr ecx, edx          ;  Find the number's top bit
  jz short sqrt_low
  shr ecx, 1            ;  It's in EDX, so k = bit/2 + 17
  neg ecx
  add ecx, 15           ;  ECX = 32 - k
  jmp short sqrt_guess
sqrt_low:
  bsr ecx, eax
  jz short sqrt_zero    ;  The root of zero is zero, and EAX is already 0
  shr ecx, 1            ;  It's in EAX, so k = bit/2 + 1
  neg ecx
  add ecx, 31           ;  ECX = 32 - k
sqrt_guess:
  shr ebx, cl

    ;  Newton's method: x = (x + n/x) / 2

sqrt_loop:
  mov eax, esi
  mov edx, edi
  div ebx               ;  EAX = n / x
  add eax, ebx          ;  EAX = (x + n/x) / 2, keeping the carry
  rcr eax, 1
  cmp eax, ebx          ;  Stop as soon as it stops getting smaller
  jae short sqrt_done
  mov ebx, eax
  jmp sqrt_loop
sqrt_done:
  mov eax, ebx
sqrt_zero:
  pop edi esi edx ecx ebx
  ret

The function takes in a 64-bit integer and produces a 32-bit result. If the 64-bit input is fixed 32:32, the output can be interpreted as fixed 16:16, which is handy because a lot of the time a square root is needed straight after a multiplication, and the result of the imul is already in the right form to be passed straight into sqrt.

The first guess makes a big difference. Starting from FFFFFFFFh every time, the guess only halves on each step until it gets near the root, and it takes around 20 divides; starting from a power of two just above the root, found from the number's top bit with bsr, it usually takes five, the last of which only confirms the answer. The guess and every step after it stay above the root, and the root of a 32:32 number is always bigger than its top half in EDX, so the div can never overflow.

Bit by bit

A bitwise square root works out the root a bit at a time, deciding in turn whether each bit in the answer should be set. The square root of a 32-bit number is no more than 0xFFFF, so it starts at the top bit and works down towards the bottom. The simple way to test a bit is to take the result so far, set the bit and square it: if the square is higher than the number, the bit shouldn't be set. This version does the same thing using only additions and shifts:

long sqrt_bit(long r)
{
    long t;                          //  temp
    long b;                          //  bit
    long c = 0;                      //  result

    for (b = 0x40000000; b != 0; b >>= 2) {
        t = c + b;
        c >>= 1;
        if (t <= r) {
            r -= t;
            c += b;
        }
    }
    return c;
}

It's the long division method for square roots, done in binary. c holds the root found so far, kept shifted up by the bit being tried, and r holds what's left of the number once the square of the root so far is taken away. For the root of 200, the first steps with b above 64 change nothing, and then:

The root is 14, and the 4 left in r is what's over: 200 = 14² + 4.

This could be used as a fast way to calculate the square root of a fixed point number. Passing in a fixed 16:16, the function returns the result in fixed 8:8, so the result has to be shifted left by 8 places, and it's only accurate to 8 fractional bits. This assembler version does that, and since the root and the bit never overlap, the additions can be ors:

;  Inputs:
;    EAX = number (fixed 16:16)
;  Returns:
;    EAX = square root (fixed 16:16, to 8 fractional bits)

sqrt16:
  push ebx ecx edx
    ;  EAX = r, EBX = b, ECX = t, EDX = c
  xor edx, edx
  mov ebx, 40000000h
sqrt16_loop:
  mov ecx, edx          ;  t = c + b
  or ecx, ebx
  shr edx, 1            ;  c >>= 1
  cmp eax, ecx          ;  if (t <= r) {
  jb short sqrt16_nobit
  sub eax, ecx          ;      r -= t
  or edx, ebx           ;      c += b
sqrt16_nobit:           ;  }
  shr ebx, 2
  jnz sqrt16_loop
  mov eax, edx          ;  Turn the 8:8 root into 16:16
  shl eax, 8
  pop edx ecx ebx
  ret

The loop could be unrolled, and the jump inside it can be got rid of with the carry flag, which makes it more Pentium-friendly. sbb of a register from itself gives -1 if the carry is set and 0 if it isn't, a mask that chooses between t and nothing:

sqrt16b_loop:
  mov ecx, edx          ;  t = c + b
  or ecx, ebx
  shr edx, 1            ;  c >>= 1
  cmp eax, ecx          ;  Carry is set if r < t...
  sbb esi, esi          ;  ...which makes ESI FFFFFFFFh, or 0 if t <= r
  not esi               ;  ESI = FFFFFFFFh if t <= r
  and ecx, esi          ;  ECX = t, or 0
  and esi, ebx          ;  ESI = b, or 0
  sub eax, ecx          ;  r -= t, or 0
  or edx, esi           ;  c += b, or 0
  shr ebx, 2
  jnz sqrt16b_loop

Because this doesn't use many registers, the function can be made to calculate two square roots at once in one loop, with the two sets of instructions taking advantage of Pentium instruction pairing.

Recursively

The recursive method uses the fact that dividing a number by four halves its square root. So it finds the square root of n >> 2, multiplies that by two and decides whether or not it needs to add one to the result. The only advantage it has over the iterative method is that it uses a multiply instead of a divide, and it returns the integer just below (or equal to) the square root:

long sqrt_rec(long n)
{
    long k;

    if (n == 0) return 0;
    k = sqrt_rec(n >> 2) << 1;      //  The root of n/4, doubled
    if ((k + 1) * (k + 1) <= n) k++;  //  One more if it still fits
    return k;
}

For 200 it goes down through 50, 12, 3 and 0, and comes back up with roots of 0, 1, 3, 7 and 14. It works for any number up to 46340²; above that, (k + 1)² no longer fits in 32 bits.

Vectors

These pieces of code are useful for 3D maths. Each takes a vector in ESI, EDI and EBP, with a second one in memory where it needs two.

Dot product

Given a pair of vectors (ax, ay, az) and (bx, by, bz), the dot product of the two is ax·bx + ay·by + az·bz. This is equal to the cosine of the angle between the vectors multiplied by their lengths, and it's useful for light sourcing calculations amongst other things: if the two vectors have unit length, as a surface normal and the direction of a light usually do, the dot product is just the cosine of the angle between them, which is how bright the surface is.

;  Inputs:
;    ESI, EDI, EBP = vector a
;    _vect_x, _vect_y, _vect_z = vector b
;  Returns:
;    EAX = a . b (fixed 16:16)

_vect_x  dd 00010000h
_vect_y  dd 00010000h
_vect_z  dd 00010000h

vect_dot:
  push ebx ecx edx
    ;  Multiply the X terms
  mov eax, _vect_x
  imul esi
  mov ebx, eax
  mov ecx, edx
    ;  Add the Y terms
  mov eax, _vect_y
  imul edi
  add ebx, eax
  adc ecx, edx
    ;  Add the Z terms
  mov eax, _vect_z
  imul ebp
  add ebx, eax
  adc ecx, edx
    ;  Take the middle 32 bits of ECX:EBX as fixed 16:16
  mov eax, ebx
  mov ax, cx
  rol eax, 16
  pop edx ecx ebx
  ret

The three products are added up as 64-bit numbers with add and adc, and only cut back to 16:16 at the end, so nothing is lost in between. (1, 2, 3) · (4, 5, 6) comes out as 00200000h, 32.

Cross product

The cross product of two vectors is a vector at right angles to both of them, (ay·bz - az·by, az·bx - ax·bz, ax·by - ay·bx). Crossing two edges of a polygon gives its normal, for the lighting above, and the sign of the normal's z says whether the polygon faces the camera or away from it.

;  Inputs:
;    ESI, EDI, EBP = vector a
;    _vect_x, _vect_y, _vect_z = vector b
;  Returns:
;    _cross_x, _cross_y, _cross_z = a x b

_cross_x dd ?
_cross_y dd ?
_cross_z dd ?

vect_cross:
  push eax ebx ecx edx
    ;  X = ay*bz - az*by
  mov eax, edi
  imul _vect_z
  mov ebx, eax
  mov ecx, edx
  mov eax, ebp
  imul _vect_y
  sub ebx, eax
  sbb ecx, edx
  mov eax, ebx
  mov ax, cx
  rol eax, 16
  mov _cross_x, eax
    ;  Y = az*bx - ax*bz
  mov eax, ebp
  imul _vect_x
  mov ebx, eax
  mov ecx, edx
  mov eax, esi
  imul _vect_z
  sub ebx, eax
  sbb ecx, edx
  mov eax, ebx
  mov ax, cx
  rol eax, 16
  mov _cross_y, eax
    ;  Z = ax*by - ay*bx
  mov eax, esi
  imul _vect_y
  mov ebx, eax
  mov ecx, edx
  mov eax, edi
  imul _vect_x
  sub ebx, eax
  sbb ecx, edx
  mov eax, ebx
  mov ax, cx
  rol eax, 16
  mov _cross_z, eax
  pop edx ecx ebx eax
  ret

Each difference is taken between the full 64-bit products with sub and sbb, like the sums in the dot product. (1, 0, 0) × (0, 1, 0) is (0, 0, 1), and (1, 2, 3) × (4, 5, 6) is (-3, 6, -3).

Length of a vector

The length of the vector (x, y, z) is √(x² + y² + z²). The squares add up in 32:32, which is exactly what the Newton square root takes:

;  Inputs:
;    ESI, EDI, EBP = vector
;  Returns:
;    EAX = its length (fixed 16:16)

vect_length:
  push ebx ecx edx

  mov eax, esi          ;  x*x
  imul esi
  mov ebx, eax
  mov ecx, edx

  mov eax, edi          ;  + y*y
  imul edi
  add ebx, eax
  adc ecx, edx

  mov eax, ebp          ;  + z*z
  imul ebp
  add eax, ebx
  adc edx, ecx

  call sqrt             ;  The sum is fixed 32:32, the root comes back 16:16

  pop edx ecx ebx
  ret

The length of (3, 4, 0) is exactly 5.0, 00050000h, and (1, 1, 1) gives 1.73204, the square root of 3 to the last bit.

Unit vector

Dividing each component of a vector by the length of the vector makes its length 1. Each component is put into 32:32 the same way as in the division above, and divided by the length:

;  Inputs:
;    ESI, EDI, EBP = vector
;  Returns:
;    ESI, EDI, EBP = the same direction, with length 1

unitize:
  push eax ebx edx

  call vect_length
  mov ebx, eax

    ;  X component
  mov eax, esi
  rol eax, 16           ;  EDX:EAX = x in fixed 32:32
  movsx edx, ax
  xor ax, ax
  idiv ebx
  mov esi, eax
    ;  Y component
  mov eax, edi
  rol eax, 16
  movsx edx, ax
  xor ax, ax
  idiv ebx
  mov edi, eax
    ;  Z component
  mov eax, ebp
  rol eax, 16
  movsx edx, ax
  xor ax, ax
  idiv ebx
  mov ebp, eax

  pop edx ebx eax
  ret

No component can be longer than the whole vector, so the divides never overflow. (3, -4, 0) becomes (00009999h, FFFF3334h, 0), which is 0.6 and -0.8 to within 1/65536.

References