A fast routine? You're probably better off doing floating point than 64 bit integer math on a processor with 16 bit multiplication.

Assuming d0 and d2 are the high words:

128 bit result: d0*d2*2^64+(d0*d3+d1*d2)*2^32+d1*d3. This uses 32 bit multiplication, a problem I assume you've already solved.

I don't have a division routine, but it sounds positively painful.

Continuing from my last post, here's what I've gotten so far. I decided to start with a plain 64x64->128 multiplication routine but with a few optimization tricks.

First of all, this routine will have to execute 16 multiply instructions, since there are four 16-bit words per 64-bit value (4*4=16), as I noted in my previous post. To make the rest of this post make more sense, I've numbered all 16 of these 32-bit products:
00 11 22 33     ae af ag ah
44 55 66 77  =  be bf bg bh
88 99 aa bb     ce cf cg ch
cc dd ee ff     de df dg dh

The first little optimization I thought to do is to push the first 4 32-bit products onto the stack (storing the 128-bit result on the stack). The first four products are 00 (88|55|22) (dd|aa|77) ff, where (...) provides a few different options in that position. I wanted to minimize the number of swaps needed to multiply these products (swaps are needed to put the appropriate 16-bit value in the lower 16 bits of the register that contains that value), so I chose the products 00 88 dd ff. This sequence requires only two swaps. (I actually calculate them in reverse order--ff dd 88 00--so that 00 ends up at the top of the stack (lowest address) because it's the most significant word, and we're storing values in big endian). From here I was just thinking about calculating the remaining 12 products and add them to the current 128-bit result, with carry, using inline loops (or unrolled loops). Either that, or push each of them onto the stack as individual partial products and add them all with carry.

The second optimization I thought of is to combine or "pack" the products together into partial products like the first one. This gave me the following seven partial products:
0088ddff
.44ccee.
.1199bb.
..55aa..
..2277..
...66...
...33...Then I would add up only seven partial products rather than 13 (or worse, 16). Note that the ordering of the products within each partial product is somewhat arbitrary. I could have used the partial products ..5577.. and ..22aa.. instead of ..55aa.. and ..22aa.., for example. They add up the same either way.

The last optimization tick comes from a recent realization of mine: the maximum product of two unsigned 16-bit values is 0xFFFE0001. What this means is that staggered long (32-bit) values can be added together without carry. Here's an example:
FFFE0001
+   FFFE0001
+       FFFE0001
----------------
FFFEFFFFFFFF0001This addition can be extended to any number of staggered longs without carry.

This optimization resulted in the following layout:
first partial product:
0088ddff +
.44ccee.

second partial product:
.1199bb. +
..55aa..

third partial product:
..2277.. +
...66...

fourth partial product:
...33...

The first partial product is calculated by adding the first set of products (00 88 dd ff) and the second set of products (44 cc ee) without carries, since they are staggered from each other. The second partial product is calculated similarly and is then added to the first partial product (using addition with carries). The third and fourth partial products are also calculated and added to the result in turn. Afterwards, the 128-bit result on the stack is popped off into the registers %d0-%d3 (with a single movem.l (%sp)+,%d0-%d3 instruction) and the routine returns.

Using this staggered-add trick, I just have to add up 4 partial products with carry instead of 16 or even 13. Hopefully I can find a few other tricks like this to help speed up the routine even more.

The stack usage of this routine will be 256 bits, or 32 bytes. I might be able to reduce that by storing part of the 128-bit result in registers and part of it on the stack, which will surely complicate it more but might also make it a bit faster as well.



UPDATE: I finally implemented a working multiply routine, almost exactly as I described it here. It takes a bit more than 32 bytes of stack space, though, since it has to save a few registers. My first test was to multiply 0x000000E8D4A51000 (1e12 decimal) by 0x000000000280DE80 (42e6 decimal) to get the product 0x000000000000000246DDF97976680000 (42e18 decimal), which the routine managed to compute correctly after a few bug fixes. I'll need to test some larger values to ensure that all bits are touched with non-zero values.

Right now the routine takes 2326 cycles, which isn't too bad but not that great either. The 16 multiplies account for just over half of that (1184 cycles), so I may be able to reduce the total by a few hundred cycles with some additional optimizing.

As I've been posting here at Cemetech, I have a working routine now. It's not necessarily the most optimized but it's pretty quick (it runs in 1824 cycles).

Here's a screenshot showing a little test program multiplying a few inputs:

If you're so inclined, you can double-check the results yourself.


On a somewhat but not too related note, the Windows calculator sucks. It can't do fractions in "Programmer" mode, and that's the only mode that supports different bases (eg, hexadecimal), which means it can't calculate hexadecimal fractions like 1.6A09E667F. Microsoft must think that programmers don't (or can't) use fractions. I'm a programmer, and I do in fact use fractions at work! It also uses weird keyboard shortcuts like Y for exponentiation or F9 for negate. The only sensible keyboard shortcuts are the keys 0 through 9 and +-/*. Lame.

Fortunately, I use Linux at home, so I don't have these problems at home. gcalctool FTW!
</rant>

Um, fractions aren't defined for INTEGER data types, which are all the Windows calculator deals with in programmer mode.

That was basically the point of my gripe: that it supports only integers in Programmer mode. It would be nice if it supported fractions in hexadecimal, octal, and binary. (Hexadecimal/octal/binary are not integer types; they are bases).

For example, if I want to see what 2.375 is in binary, I use the standard Gnome calculator (gcalctool) and switch to binary. 2.375 decimal is 10.011 in binary, 2.6 in hexadecimal, or 2.3 in octal. It's really handy for calculating fixed-point values, for example (which is a common task in some areas of PROGRAMMING), but the Windows calculator can't do this basic task.

Quote from: Qwerty.55 on August 07, 2011, 12:14:15 am

Um, fractions aren't defined for INTEGER data types, which are all the Windows calculator deals with in programmer mode.

That was basically the point of my gripe: that it supports only integers in Programmer mode. It would be nice if it supported fractions in hexadecimal, octal, and binary. (Hexadecimal/octal/binary are not integer types; they are bases).

For example, if I want to see what 2.375 is in binary, I use the standard Gnome calculator (gcalctool) and switch to binary. 2.375 decimal is 10.011 in binary, 2.6 in hexadecimal, or 2.3 in octal. It's really handy for calculating fixed-point values, for example (which is a common task in some areas of PROGRAMMING), but the Windows calculator can't do this basic task.

I know that hexadecimal and the rest are bases. I've worked in some pretty crazy bases before (base -5i, anyone?)

All I was saying was that the Windows calculator only handles one data type in that particular mode and the operations you'd like aren't defined for that data type. Calculator is like Notepad. It works, but it doesn't have a lot of nice features. I'm frankly surprised they even put Integer data types in there. Personally, I never need more than integer types, but a cursory glance at Google shows some calculators you can use.

How do you work in base -5i

so 17 = -(4i * (5i)¹ + 3 * (5i)⁰) = -[4i, 3]_5i?

I couldn't come up with an expression which used natural numbers lower then 5 only, but I may be missing something...