Rudiments of computer math

The inner core of floating point
Charles Douglas Wehner
30 December 2003

These notes are for guidance only.
No liability accepted.

The programs can be assembled with
the free NASM assembler from GNU.

The art of good programming is to lose the plot. That is to say, if a program is written for ease of understanding by the programmers it will probably contain few ingenious tricks.

If, however, several things are taking place at once, the CPU will work at breathtaking speed towards its goal. It is not necessary for the programmer to continue to remember what he wrote - it is quite sufficient that these details are encapsulated in the program.

However, good documentation is essential for the maintenance of the system.

An example of a good programming trick arose when a friend known as Stavros showed the author his university assignment. He had used a constant and a slope to find the first guess for the Newton-Raphson square-root.

You MULTIPLY by the slope, and ADD the constant. This gives accuracy sufficient for IEEE Real 32 in just three turns of the loop.

The author observed that it took perhaps TWO-AND-A-HALF turns. This cannot be capitalised upon because two would be insufficient whilst three is a little too good.

The author decided to SHIFT THE PACKED R32 ONE PLACE RIGHT.

Bizarre! And yet it works.

Floating point consists of a number often called the MANTISSA (defined as something of little magnitude or importance) and an EXPONENTIATING number (often wrongly called the EXPONENT).

The "EXPONENT" is usually binary. That is to say, it defines how many times the Mantissa must be shifted left or right in binary for it to become the correct number.

For example, if E=4 and M=1.5 the number is 24. Each left shift doubles the Mantissa. Thus, we have 1.5 shifted to 3 6 12 and 24.

The observant will notice that from this definition, the "Exponent" is actually the INTEGER of the LOGARITHM BASE 2 of the number.

To cover the special case of the number being zero, it is usual for the Exponent to be set to zero. So we add an offset to the Exponent. As a result, zero binary shifts will not be confused with the number zero. In Real 32, this offset is 127. Notice how even Exponents have become odd and vice-versa due to the odd offset.

It is usual for the Mantissa to be left-justified in binary. That means that the leftmost bit is ALWAYS set. So 1.5 has a Mantissa of 11000000 00000000 00000000.

The leftmost bit means "1", the next means "0.5" and so on, halving with each step to the right.

Because the leftmost bit is ALWAYS set, it is wasted. So once we have reached this stage (having packed the shifts with offset into the Exponent and having a left-justified Mantissa) we may as well RESET this bit when the number is positive, or LEAVE IT SET for negative numbers. This is the USUAL format of floating-point.

NOT SO with IEEE math! With Real 32, the ENTIRE Exponent byte is shifted right until the bottom bit overlays the top bit of the Mantissa. The Signum (sign of the number) is specified by resetting the TOP of the Exponent for positive, or setting it for negative.

You might now understand why the author dislikes IEEE math. You cannot reclaim your Exponent without bit-shifting.

However, 1.5 with an odd number of binary shifts will have a Mantissa 11000000 00000000 00000000 whilst 1.5 with an even number of shifts has 01000000 00000000 00000000.

Shift these right and you get 01100000 00000000 00000000 and 00100000 00000000 00000000 respectively.

The "VIRTUAL BIT" cannot be shifted. Whatever happens, the leading bit always represents 1. So we have 11100000 00000000 00000000 and 10100000 00000000 00000000 for the two first-guesses.

These numbers are 1.75 and 1.25 respectively.

So we automatically generate a first guess for root 3 of 1.75. The correct root of 3 is about 1.732. Similarly, for root 1.5 we have 1.25. The correct root is about 1.2247.

In addition, the powers of two are also halved. So 3 times 2 to the power 10 gives 1.75 times 2 to the power 5. 3 times 2 to the power 9 gives 1.25 times 2 to the power 4.

These are not as good as the Stavros guesses - good for two and a half turns. These guesses are good for perhaps two-and-three-quarters.

But we CANNOT go round the loop a non-integer number of times. We go three times around the loop as with Stavros. But we have replaced his MULTIPLY and ADD by a single 32-bit SHIFT.

That is the name of the game.

Lost the plot? The PROGRAM has not. I used this trick in many versions, and the computers always delivered a square root in record time. And if the program needs maintenance, one can always think it through again.

It is very convenient for the lay user to have floating point for the generation of sines and cosines. However, the writers of commercial graphics programs usually construct their line-drawing routines from INTEGER mathematics. If they do not, the program runs too slowly.

The inner core of floating point runs at the integer level. A central processor does not understand floating point. Registers contain binary switches.

However, pressure from lazy programmers forced Intel to incorporate floating point routines into the Pentium range. This has the commercial advantage that users of the facilities find themselves locked into the Intel system, with programs that are not portable to other processors.

In the design of the inner core of floating point, we generate routines that can equally well be applied to integer graphics systems. However, graphics require scaled-down versions because 32-bit precision is excessive.

Cosine

The cosine of a negative number is the same as that of the equivalent positive number. Clear the signum (topmost bit in IEEE Real, topmost of the Mantissa in other systems).

Now divide by Pi. Look at the lowest bit of the quotient. If it is set, it represents half-cycles 2, 4, 8 etc. These half-cycles are the mirror-image of the others. Use CMA to complement the accumulator, or "XOR eax,-1" or similar. The numbers should now run backwards for these half-cycles.

The Mantissa should represent PAIRS of quarter-cycles (quadrants). The first quadrant is positive. The second quadrant gives a NEGATIVE result, as well as running backwards.

So you shift left into carry, or "ADD eax,eax" to convert the Mantissa into a single quadrant. Also, PRESERVE the carry. If it was set, you will CURRENTLY complement the accumulator ("XOR eax,-1") and LATER negate the result.

The numbers represent 0 to 0.99999999976716935634613037109375 of a quadrant. They will have to be multiplied by 3373259426. This is the fixed-point representation of Pi/4. The product will be in EDX:EAX on a Pentium. So we move the EDX number into EBX for the following routine:

Remember to remove the test loadings (MOV EBX,1686629713 &c.). The program is 32 bit IBM-compatible only.

start: mov ebx,1686629713 ;; 22-5 degrees mov ebx,2248839617 ;; 30 degrees mov ebx,3373259426 ;; 45 degrees mov eax,ebx mul ebx mov ebx,edx mov eax,9 mul ebx mov eax,1184 sub eax,edx mul ebx mov eax,106522 sub eax,edx mul ebx mov eax,5965232 sub eax,edx mul ebx mov eax,178956971 sub eax,edx mul ebx mov eax,2147483648 sub eax,edx mul ebx mov eax,-1 sub eax,edx

Amended 3 December 2006

display: push eax mov dx,leader+256 mov ah,9 int 33 pop eax mov cx,32 thirtytwo: push cx xor dl,dl mov ebx,eax add eax,eax adc dl,dl add eax,eax adc dl,dl add eax,ebx adc dl,0 add eax,eax adc dl,dl push eax add dl,48 mov ah,2 int 33 pop eax pop cx loop thirtytwo int 32 leader: db "0.$"

The display routine shows everything in the register. For full precision, a 32-bit number must be displayed as 32 figures in any number system other than powers of two. This does not mean that 32 decimal places are available. There are NINE decimal places, suitable for an EIGHT-FIGURE display.

But you will ask "WHAT ARE THE MAGIC NUMBERS"?

First, we have Horner`s rule. Isaac Newton used to evaluate polynomials backwards:

______________________________________
______________________________
______________________
______________
______
1 / 6 + 1 / 5 + 1 / 4 + 1 / 3 + 1 / 2 + 1

This leads to the magic number e. The author did the same kind of thing, but according to Professor Donald E Knuth, Horner published first.

Little Jack Horner
Sat in a corner
Eating his Christmas π
He put in his thumb
And pulled out an algorithm
And said "What a good boy am I".

Given Y=a plus bX plus cX-squared plus dX-cubed

Take d and multiply by X
Add to c and multiply by X
Add to b and multiply by X
Add to a and STOP

Given Y=a minus bX plus cX-squared minus dX-cubed

Take d and multiply by X
Subtract from c and multiply by X
Subtract from b and multiply by X
Subtract from a and STOP

So Horner`s rule says RUN THE POYNOMIAL BACKWARDS. Here we have the Cosine equation from Colin McLaurin and Brook Taylor, where

Cos(X)= 1 - X-squared/2 + X-fourth/24 - X-sixth/ 720 etc.

So we prepare X-squared by "mov eax,ebx ; mul ebx ; mov ebx,edx" and then use integers representing 1/2, 1/24 and so on - but in reverse order.

They were made by taking 4294967296 and dividing it by 2. That result was divided by 3 and then by 4. This next result was divided by 5 and by 6. And so on. Only the integer part was kept.

So if we multiply a 32-bit integer in EBX by the number 178956971 in EAX, for example, the product appears in EDX:EAX. Taking EDX down to EAX (transfering the number) has the effect of dividing by 4294967296. So we get 178956971/4294967296 or one 24th.

But the design of a computer program depends on the parameters of the CPU. Does it waste time to load the instruction "mov eax,178956971"? Is it quicker just to load 178956971 by means of "LODSD"?

Here is that same program as a loop. It takes up less space.

start: mov ebx,1686629713 ;; 45 degrees mov ebx,2248839617 ;; 60 degrees mov ebx,3373259426 ;; 90 degrees mov eax,ebx mul ebx mov ebx,edx mov si,coefficients+256 mov cx,7 xor edx,edx cld jmp hopin cosloop: mul ebx hopin: lodsd sub eax,edx loop cosloop ret

display program as before. leader: db "0.$" coefficients: dd 9 dd 1184 dd 106522 dd 5965232 dd 178956971 dd 2147483648 dd -1

If it is desired not to divide by Pi/4, this can be done by altering the coefficients. Here the coefficients are increased in size by Pi-squared/16 once for the second coefficient, twice for the third, three times for the fourth and so on.. However, it takes one more turn of the loop, one more coefficient.

start: mov ebx,2147483648 ;; 45 degrees mov ebx,2863311531 ;; 60 degrees mov ebx,-1 ;; 90 degrees mov eax,ebx mul ebx mov ebx,edx mov eax,27 mul ebx mov eax,2023 sub eax,edx mul ebx mov eax,108241 sub eax,edx mul ebx mov eax,3948192 sub eax,edx mul ebx mov eax,89607967 sub eax,edx mul ebx mov eax,1089502240 sub eax,edx mul ebx mov eax,1003736220 sub eax,edx mul ebx mov eax,-1 sub eax,ebx sub eax,edx

display program as before.

The Sine can be created by means of
Sin X = X - X-cubed/3! + X-fifth/5! &c.
or by entering into the cosine program. For example, subtract Pi/4 radians from the argument and treat as a cosine.

Unfortunately, accuracy is lost for very small arguments. As a result, one should first check that the argument is greater than that which gives a different result to its argument. These very small arguments are returned unchanged, and only the larger arguments resolved via the polynomial.

There is a great deal that could be said on the subject, but as the title says, these are the RUDIMENTS of floating point.

Here are a couple of small sinusoid programs (sine and cosine) suitable for graphics or for the fast fourier transform of JPEG:

(COSINE) start: mov bx,32768 mov ax,bx mul bx mov bx,dx mov ax,91 mul bx mov ax,2730 sub ax,dx mul bx mov ax,32768 sub ax,dx mul bx xor dx,-1


(SINE)

start:
   mov bx,32768
   mov cx,bx

   mov ax,bx
   mul bx
   mov bx,dx

   mov ax,13
   mul bx
   mov ax,546
   sub ax,dx
   mul bx
   mov ax,10922
   sub ax,dx
   mul bx
   mov ax,-1
   sub ax,dx
   mul cx

display: push dx mov dx,leader+256 mov ah,9 int 33 pop ax mov cx,16 sixteen: push cx xor dl,dl mov bx,ax add ax,ax adc dl,dl add ax,ax adc dl,dl add ax,bx adc dl,0 add ax,ax adc dl,dl push ax add dl,48 mov ah,2 int 33 pop ax pop cx loop sixteen int 32 leader: db "0.$"

It will be seen that the sine was generated by evaluating
1 - X-squared/3! + X-fourth/5! - X-sixth/7!
which for small arguments approximates to 1

Thereafter, the result is multiplied by X. For a Fourier or line-drawing program, it suffices to multiply by the SHIFTED version of X. For example, if 10000000 0000000 0000000 represents the Mantissa of an EIGHTH, 00010000 00000000 00000000 represents the shifted number.

The right-shift will have caused bits to fall off the register - and these are lost forever. Accordingly, in floating point we would keep the 10000000 00000000 00000000 for the final multiply, whilst turning it into 00010000 00000000 00000000 to square to 00000010 00000000 00000000 for the routine.

When the result is multiplied by the unshifted Mantissa, the old Exponent becomes the Exponent of the result. However, it is always possible that 10000000 00000000 00000000 or whatever, when multiplied in, gives a leading zero. Floating point requires that the Mantissa is kept left-justified.

So in the 32-bit routines, we would do the following:

evaluate: push ebx SHIFT EBX ACCORDING TO EXPONENT mov eax,ebx mul ebx mov ebx,edx RUN THE HORNER`S RULE ROUTINE pop ebx mul ebx and edx,edx js continue add eax,eax adc edx,edx DECREMENT THE EXPONENT CREATE FP-NIL IF IT DOES NOT DECREMENT continue:

Arctangent

The arctangent can be created by means of
ATN X = X - X-cubed/3 + X-fifth/5 &c.

This has all the appearances of the sine, except for one small difference - there is no exclamation mark!!!!!!!!!!!!!!!!!!

So the factorial convergence of the sinusoids has gone. Instead of dividing by 1, by 6, by 120, by 5040 and 362880, we use 1, 3, 5, 7 and 9. This is a disaster.

The only solution is to use very small arguments. So if the number is about 1/2, we will be multiplying by about 1/4 at every "turn" of the loop. This delivers exponential convergence instead of factorial convergence.

To begin with, the tangent of 45 degrees is unity. Numbers greater than unity can be turned into their RECIPROCALS, and the routine can be CALLED rather than jumped to. On return, we have the chance to subtract the result from 90 degrees to obtain the correct result.

Numbers from zero to 1 (the tangents of zero to 45 degrees) have to be reduced in size. Fortunately, it can be shown that the tangent of 22.5 degrees is ROOT-2 MINUS 1.

"IT CAN BE SHOWN" is the lazy mathematician`s cliche. OK. Here goes.

The tangent of (A minus B) is the DIFFERENCE of Tan(A) and Tan(B), divided by Tan(A) TIMES Tan(B) plus 1.

So if A is 45 degrees, and B is 22.5 we have
Tan(22.5) = 1 - Tan(22.5) / 1 + Tan(22.5)
Tan(22.5) + (Tan(22.5))squared = 1 - Tan (22.5)
(Tan(22.5))squared + 2Tan(22.5) - 1 = 0

There is a formula for resolving equations of the form
aX-squared + bX + c = 0
It is:
X = (-b +/- Root(b-squared - 4ac) )/ 2a

This gives -2 +/- ROOT(4 + 4) / 2

Which is -1 +/- ROOT(2)

So Tan(22.5) is 0.414213562, which when multiplied by 4294967296 gives 1779033704 as a magic number for our program.

If the argument is less than this, we JUMP to the program. If it is more, we have to reduce the argument and CALL. After the call we will add 22.5 degrees to the result.

Here is a piece of that code:

mov eax,1779033704 (Tan 22.5) sub ecx,eax (DIFFERENCE) mul ebx (PRODUCT OF ARGUMENT AND 0.414) stc (SET THE CARRY) rcr edx,1 (CREATES 1+PRODUCT /2) mov ebx,edx (SAVE IT) mov edx,ecx (FETCH PRODUCT OF ARGUMENT AND 0.414) div ebx (DIVIDE BY HALF 1+PRODUCT) clc (RESULT IS DOUBLE SIZE) rcr eax,1 (SO HALVE THE RESULT)

This then will deliver a reduced argument, the tangent of the "ANGLE LESS 22.5", without the angle being known. On return, the result has 1686629713 added on. This is 22.5 degrees in radians, as a fraction of 4294967296.

start: mov ebx,1073741824 ;; QUARTER mov ebx,2147483648 ;; HALF mov ebx,-1 ;; Tan 45 call arctangents call display int 32

arctangents: mov eax,1779033704 sub eax,ebx jnc evaluate mov ecx,ebx mov eax,1779033704 sub ecx,eax mul ebx stc rcr edx,1 mov ebx,edx mov edx,ecx div ebx clc rcr eax,1 mov ebx,eax call evaluate add edx,1686629713 ret

evaluate: push ebx mov eax,ebx mul ebx mov ebx,edx mov eax,204522252 mul ebx mov eax,226050910 sub eax,edx mul ebx mov eax,252645135 sub eax,edx mul ebx mov eax,286331153 sub eax,edx mul ebx mov eax,330382100 sub eax,edx mul ebx mov eax,390451572 sub eax,edx mul ebx mov eax,477218588 sub eax,edx mul ebx mov eax,613566757 sub eax,edx mul ebx mov eax,858993459 sub eax,edx mul ebx mov eax,1431655765 sub eax,edx mul ebx mov eax,-1 sub eax,edx pop ebx mul ebx ret

When creating a giant laboratory mathematical system, for example, the constants of floating-point will not be known. It is the job of the programmer to conjure them up.

The author wrote first a 32-decimal-place system, and then a 65-place one.

You begin by creating a four-function system, just like the cheapest pocket calculator - but more accurate. Then you use it to create ROOT 2, by starting with 1.5 as a first guess. It does not matter how slow the system is whilst you are developing it. You just want the tangent of 22.5 degrees to the full precision.

The author also generated the tangent of ELEVEN-AND-A-QUARTER degrees for use with the 32-place mathematical system, and FIVE-AND-FIVE-EIGHTH for the 65-place system

Here is the algorithm for finding the tangent of half the angle that you have the tangent of:

1. Think of a number, the tangent.
2. Square it.
3. Add 1.
4. Find the square root.
5. Subtract 1.
6. Divide by the number you first thought of.
eg:
1. 1
2. 1
3. 2
4. 1.414213562
5. 0.414213562
6. 0.414213562
(Tangent of 22.5 degrees)
again:
1. 0.414213562
2. 0.171572875
3. 1.171572875
4. 1.0823922
5. 0.0823922
6. 0.198912367
(Tangent of 11.25 degrees)
again:
1. 0.198912367
2. 0.039566129
3. 1.039566129
4. 1.019591158
5. 0.019591158
6. 0.098491403
(Tangent of 5.625 degrees)
This box was added 9 Nov 2005.

Notice how the tangent of five and five-eighth degrees is just under one tenth. This means that it will converge by one decimal place of itself. As the arctangent polynomial uses odd powers of the argument - advancing in twos, that it - the polynomial will deliver two decimal places per coefficient when five and five-eighths degrees or less are used.

The snippet

arctangents:
   mov eax,1779033704
   sub eax,ebx
jnc evaluate
   mov ecx,ebx
   mov eax,1779033704
   sub ecx,eax
   mul ebx
   stc
   rcr edx,1
   mov ebx,edx
   mov edx,ecx
   div ebx
   clc
   rcr eax,1
   mov ebx,eax
   call evaluate
   add edx,1686629713
   ret

can be rewritten for the higher precision, with new MAGIC NUMBERS. These routines can be NESTED. So Pi/8, Pi/16 and Pi/32 are added to the angle as required by the math.

Notice that the Mantissa of Pi/8, Pi/16 or Pi/32 is always that of Pi. It is slower but more accurate to use less nested argument-reduction routines. So when evaluating Pi for laboratory use, one might drive the big math from Basic, and use ROOT-TWO MINUS 1 only. When the big math has been built, the deeper levels of argument-reduction are added in order to speed it up.

Logarithm

It was pointed out that the "Exponent" byte is the INTEGER of the logarithm to base 2 of a floating-point number. This can be reclaimed, and we only have the fractional part to find.

The fractional part will be the LOG base 2 of the Mantissa. Put them both together, float the result and you have the LOG base 2.

What is it good for? Well, logarithms are proportional to one another. Multiply any logarithm base 2 by 0.69314718 - the natural log of 2 - and you have converted that logarithm.

The ANTILOG base 2 of 0.5 is 1.414213562, root 2. So if you multiply the Mantissa by itself and it exceeds 2, it will have contained root 2. If it takes TWO multiplications to reach 2, it contained the FOURTH root of 2, which is 0.25 in base-2 logarithms.

So by repeated squaring of the Mantissa we can generate the LOG base 2 bit-by-bit:

start: mov edx,3373259426 xor ebx,ebx mov cx,32 mulloop: mov eax,edx mul edx mov eax,edx add eax,eax jc both add edx,edx both: adc ebx,ebx loop mulloop

display: push ebx mov dx,leader+256 mov ah,9 int 33 pop eax mov cx,32 thirtytwo: push cx xor dl,dl mov ebx,eax add eax,eax adc dl,dl add eax,eax adc dl,dl add eax,ebx adc dl,0 add eax,eax adc dl,dl push eax add dl,48 mov ah,2 int 33 pop eax pop cx loop thirtytwo int 32 leader: db "0.$"

If the processor has a very fast multiply routine, this is quite adequate. The concept is also good when building a giant math, because speed is not always important.

The reader will have noticed that much of computer transcendental math consists of MAGIC NUMBERS. By means of the McLaurin-Taylor equation for antilogarithm one can obtain e to any accuracy. It will be 2.718281828459045 or thereabouts. The LOG base 2 of this gives the RECIPROCAL of the LOG base e of 2.

Log base e of 2 was required to convert between bases, as described. So slow technology that needs no MAGIC NUMBERS becomes ENABLING TECHNOLOGY to find the magic numbers for the faster math.

Supposing we take the logarithm of 1.5, of 1.25, of 1.125, of 1.0625 and so on. Suppose also that we multiply these by 4294967296. We get MAGIC NUMBERS.

These represent the sizes of numbers in the Mantissa of a logarithm equivalent to multiplying by 1.5, by 1.25 and so on.

In the following, 0.69314718 is represented by 2977044472. So we can expect a Mantissa of 1.999999999 to deliver this result.

However, are there any GAPS? We have assumed that the Mantissa bytes are FFFFFFFF, but if we shift and add the Mantissa without causing an overflow, we have found a gap.

The LOG base e of that shift-and-add must then be SUBTRACTED from the result. We leave the shift-and-add added into the argument because we do not want to rediscover the same gap. However, the Mantissa with its gap plugged needs to be re-shifted. Hence the "each" loop.

To use this routine, you use the Mantissa of the number. You remove the offset from the binary Exponent, and multiply by 0.69314718. To this you add the result of this routine. Then you float the lot.

start: mov ebx,2699853634 ;; 1.5 times 1.25 mov ebx,3205727920 mov ebx,2147483648 mov ebx,3221225472 mov ebx,-1 mov edx,2977044472 push edx mov edx,ebx mov si,coefficients+256 mov ecx,32 cld logloop: shr edx,1 lodsd add ebx,edx jc undo pop edx sub edx,eax push edx mov edx,ebx push cx neg cx add cx,33 each: shr edx,1 loop each pop cx nojob: loop logloop pop edx

display: push edx mov dx,leader+256 mov ah,9 int 33 pop eax mov cx,32 thirtytwo: push cx xor dl,dl mov ebx,eax add eax,eax adc dl,dl add eax,eax adc dl,dl add eax,ebx adc dl,0 add eax,eax adc dl,dl push eax add dl,48 mov ah,2 int 33 pop eax pop cx loop thirtytwo int 32 undo: sub ebx,edx jmp nojob leader: db "0.$" coefficients: dd 1741459379 dd 958394255 dd 505874286 dd 260380768 dd 132163268 dd 66589974 dd 33424039 dd 16744533 dd 8380427 dd 4192257 dd 2096640 dd 1048448 dd 524256 dd 262136 dd 131070 dd 65536 dd 32768 dd 16384 dd 8192 dd 4096 dd 2048 dd 1024 dd 512 dd 256 dd 128 dd 64 dd 32 dd 16 dd 8 dd 4 dd 2 dd 1

Will you please stop looking up there - I`m down here. Yes, I know what has distracted you. It`s those numbers 1, 2, 4, 8.... They REALLY are powers of two. It REALLY is bizarre.

When a moiety is added to a logarithm, it is equivalent to multiplying a number by 1 PLUS a moiety. But it is specifically the NATURAL logarithm where for small moieties they are the SAME. So as the number has more and more zeroes after the 1, as in 1.00000000xyz, the logarithm tends to linearity (0.00000000xyz).

We can use this fact to make the routine run many times as fast. All the gaps in the top two bytes will have been filled in after just 16 loops, and the remaining ones can simply be transferred across from the argument with add edx,ebx at the end of the routine. Note the truncated MAGIC NUMBER series.

start: mov ebx,-1 mov ebx,2699853634 ;; 1.5 times 1.25 mov ebx,3205727920 mov ebx,2147483648 mov ebx,3221225472 mov edx,2977044472 push edx mov edx,ebx mov si,coefficients+256 mov ecx,16 cld logloop: shr edx,1 lodsd add ebx,edx jc undo pop edx sub edx,eax push edx mov edx,ebx push cx neg cx add cx,17 each: shr edx,1 loop each pop cx nojob: loop logloop pop edx add edx,ebx

The importance of truncating the sequence may not be immediately obvious until it is pointed out that the shift is often an N-SQUARED function.

There exists the shift instruction

    shr edc,ecx

on some Intel CPUs. This will presumably wire the input at edx directly to the output at edx. It might just take one or two ticks of the system clock. Without any reference books to guide him, the author tried it and it failed to work.

In order to make progress, the author wrote

 each:
   shr edx,1 
 loop each

which is GUARANTEED to take increasing time for increasing contents of the loop-counter ECX.

So the number of shifts, if all possible shifts took place, would be (N)(N+1)/2, which for 32 is 528 but for 16 is just 136. This means that the second program is almost four times faster.

Antilogarithm

We have seen how to generate a cosine. This uses the McLaurin-Taylor equation. There is an exact equivalent for the natural antilogarithm:
Exp(X) = 1/0! + X/1! + X-squared/2! + X-cubed/3! + ....
Amended 3 December 2006

We would prepare the magic numbers by dividing 4294967296 by the factorial of 12, of 11, of 10, and so on down to 2. The cosine program itself can be used, but with the code "mov eax,ebx, mul ebx, mov ebx,edx" removed. The number of multiplications is extended to twelve.

A reversal of the above shift-and-add method is more elegant. One cannot give hard and fast rules as to which is the BETTER, because much depends upon the parameters of the processor.

In the following program, we do not fill in the gaps in the logarithm but subtract the logarithms of 1.5, 1.25, 1.125 from it if we can. Each time we succeed, say with log 1.125, we have to multiply the result, say by one and an eighth.

start: mov ebx,1073741824 mov ebx,2699853634 ;; 1.5 times 1.25 mov ebx,3205727920 mov ebx,2147483648 mov edx,1073741824 mov si,coefficients+256 mov ecx,32 cld exploop: lodsd cmp ebx,eax jc nojob sub ebx,eax push cx neg cx add cx,33 mov eax,edx each: shr eax,1 loop each pop cx add edx,eax nojob: loop exploop

display: push edx mov dx,leader+256 mov ah,9 int 33 pop eax mov cx,32 thirtytwo: push cx xor dl,dl mov ebx,eax add eax,eax adc dl,dl add eax,eax adc dl,dl add eax,ebx adc dl,0 add eax,eax adc dl,dl push eax add dl,48 mov ah,2 int 33 pop eax pop cx loop thirtytwo int 32 leader: db "0.$" coefficients: dd 1741459379 dd 958394255 dd 505874286 dd 260380768 dd 132163268 dd 66589974 dd 33424039 dd 16744533 dd 8380427 dd 4192257 dd 2096640 dd 1048448 dd 524256 dd 262136 dd 131070 dd 65536 dd 32768 dd 16384 dd 8192 dd 4096 dd 2048 dd 1024 dd 512 dd 256 dd 128 dd 64 dd 32 dd 16 dd 8 dd 4 dd 2 dd 1

Whereas the logarithm tends to 1 PLUS UNITY-times the argument as it approaches zero, the antilogarithm tends to its own SLOPE.

In the jargon, d/dx(exp(x)) tends to exp(x).

So we do not just add the sixteen linear bits onto the result, as with the logarithm. Instead, we multiply them by the result obtained so far (by the SLOPE), and then add on.

start: mov ebx,1073741824 mov ebx,2699853634 ;; 1.5 times 1.25 mov ebx,3205727920 mov ebx,2147483648 mov edx,1073741824 mov si,coefficients+256 mov ecx,16 cld exploop: lodsd cmp ebx,eax jc nojob sub ebx,eax push cx neg cx add cx,17 mov eax,edx each: shr eax,1 loop each pop cx add edx,eax nojob: loop exploop mov ecx,edx mov eax,edx mul ebx add edx,ecx

display: push edx mov dx,leader+256 mov ah,9 int 33 pop eax mov cx,32 thirtytwo: push cx xor dl,dl mov ebx,eax add eax,eax adc dl,dl add eax,eax adc dl,dl add eax,ebx adc dl,0 add eax,eax adc dl,dl push eax add dl,48 mov ah,2 int 33 pop eax pop cx loop thirtytwo int 32 leader: db "0.$" coefficients: dd 1741459379 dd 958394255 dd 505874286 dd 260380768 dd 132163268 dd 66589974 dd 33424039 dd 16744533 dd 8380427 dd 4192257 dd 2096640 dd 1048448 dd 524256 dd 262136 dd 131070 dd 65536

By the use of exactly the same Magic Numbers for log and antilog, further savings can be made. We save 64 bytes by using the same table.

Square root

There are two aspects to creating an efficient square root. Firstly, there must be a good first guess. We studied the case of a number such as 1.21 having its virtual 1 removed, giving 0.21, having this halved to 0.105 and having the 1 restored to give 1.105.

This style of first guess was the author`s answer to Stavros.

It was pointed out that the exponent must be halved also, and that the two processes can be combined.

However, if the argument itself is not altered, we will be heading towards the same result for even and odd powers of two.

It is unavoidable that for even powers of two we must shift the argument one place to the right. This loses a single bit, equivalent to half a bit in the square root.

These routines are assumed to be 9.6329 places of decimals to be used with an 8-decimal display, and having about one and a half places of hidden precision. This is quite usual, and gives us the freedom to have an error of plus or minus 2 at the internal level.

Tests would normally be carried out to ensure that one can find the root, the root again X number of times and then square the result X number of times without corrupting the data.

It is far better to work this way than to fuss over an error of 1, to slow the system with correction routines, and still end up with an ill-conditioned math that corrupts the data when repeated.

Here is a fast square root that puts it all together:

exponent equ 4 mantissa equ 4160749568 start: mov eax,exponent mov ebx,mantissa mov ecx,ebx add ecx,ecx shr eax,1 jc odd shr ebx,1 clc odd: rcr ecx,1 push eax stc rcr ecx,1 mov edx,ebx div ecx add ecx,eax rcr ecx,1 mov edx,ebx div ecx add ecx,eax rcr ecx,1 mov edx,ebx div ecx add ecx,eax rcr ecx,1 call display pop eax call single int 32

display: push ecx mov dx,leader+256 mov ah,9 int 33 pop eax mov cx,32 thirtytwo: push cx xor dl,dl mov ebx,eax add eax,eax adc dl,dl add eax,eax adc dl,dl add eax,ebx adc dl,0 add eax,eax adc dl,dl push eax add dl,48 mov ah,2 int 33 pop eax pop cx loop thirtytwo ret single: push eax mov dx,crlf+256 mov ah,9 int 33 pop edx add dl,48 mov ah,2 int 33 ret leader: db "0.$" crlf: db 13,10,"$"

Tangent

The tangent is sine over cosine.