# Multiplication 128 x 64 bits

I was fooling around with the Collatz sequence a bit and found out that long is only sufficient for starting values below 8528817511L. As I wanted to optimize my code from this answer for computations up to $10^{10}$, I needed either BigInteger or to run my own. As the range is limited and just few operations are needed, I tried my own. My special-case multiplication is 5x faster, so I guess it wasn't a useless exercise.

I'm interested in optimization tips and general comments (apart from formatting as it's fine, at least according to my conventions).

Below I'm giving only the multiplication code and supporting methods. The whole class can be found on github.

/**
* A class representing a 127 bit positive number
* (negative numbers are unsupported and the highest bit is currently unusable).
* Only those operations which were needed for this toy project are implemented.
*
* <p>Speed is the goal, general usability is not given.
*/
@Getter @EqualsAndHashCode final class MutableLong128 implements Cloneable {
/**
* Multiply {@code this} by the argument which must be non-negative.
* Overflow may or may not be detected.
*/
public void multiply(long x) {
checkArgument(x>=0);

long y0 = low;
long y1 = high;

// The variables p0 to p3 represent the product where each of them is to be treated as an unsigned quantity.
// The weight of pi is 2**32 rather then 2**64.

long p0 = unsignedLow(x) * unsignedLow(y0);
long p1 = unsignedLow(x) * unsignedHigh(y0);
long p2 = unsignedHigh(x) * unsignedHigh(y0);

// Before another product can be added to p1, its higher part must be transferred to p2.
p2 += unsignedHigh(p1);
p1 = unsignedLow(p1) + unsignedHigh(x) * unsignedLow(y0);

// p3 represents the highest part and a possible overflow gets ignored.
long p3 = unsignedLow(x) * unsignedHigh(y1) + unsignedHigh(x) * unsignedLow(y1);

// Similarly preparing p2 to absorb another product.
p3 += unsignedHigh(p2);
p2 = unsignedLow(p2) + unsignedLow(x) * unsignedLow(y1);

assert assertCorrectProduct1(new long[] {p0, p1, p2, p3}, x, y0, y1);

// Final cleaning of upper halves.
p1 += unsignedHigh(p0);
p0 = unsignedLow(p0);
p2 += unsignedHigh(p1);
p1 = unsignedLow(p1);
p3 += unsignedHigh(p2);
p2 = unsignedLow(p2);

assert assertCorrectProduct1(new long[] {p0, p1, p2, p3}, x, y0, y1);

low = composeLowHigh(p0, p1);
high = composeLowHigh(p2, p3);
}

public BigInteger toBigInteger() {
BigInteger result = BigInteger.valueOf(low);
if (high==0) return result;
assert high>0;
}

/** If an overflow happened, {@code this} has already become unusable, but who cares. */
private void checkOverflow() {
if (high<0) throw new ArithmeticException("Overflow or negative number");
}

/**
* Used to verify the intermediate steps of the multiplication.
*
* The array represents the product starting with the least significant part.
* The parts are threated as unsigned 64 bit quantities.
* The weight of each subsequent part is 2**32 times bigger (not 2**64 times).
*
* The other arguments are just like in {@link #setToProduct1(long, long, long)}.
*/
private boolean assertCorrectProduct1(long[] parts, long x, long y0, long y1) {
final BigInteger expected = new MutableLong128(y0, y1).toBigInteger().multiply(BigInteger.valueOf(x));
final BigInteger actual = composeBigInteger(parts);
assert actual.equals(expected) : actual.toString(16) + " " + expected.toString(16);
return true;
}

/**
* The array represents the result starting with the least significant part.
* The parts are threated as unsigned 64 bit quantities.
* The weight of each subsequent part is 2**32 times bigger (not 2**64 times).
*/
private BigInteger composeBigInteger(long[] parts) {
BigInteger result = BigInteger.ZERO;
for (int i=parts.length; i-->0; ) {
final long p = parts[i];
}
return result;
}

/** Extract the 32 least significant bits treated as unsigned number. */
private static long unsignedLow(long x) {
return x & 0xFFFFFFFFL;
}

/** Extract the 32 most significant bits treated as unsigned number. */
private static long unsignedHigh(long x) {
return x >>> 32;
}

/** Compose a long from two unsigned parts. The upper 32 bits of low must be zero. */
private static long composeLowHigh(long low, long high) {
return (high << 32) + low;
}

private static final BigInteger TWO_64 = BigInteger.ONE.shiftLeft(64);

/** The lower 64 bits treated as unsigned. */
private long low;
/** The upper 64 bits (currently only 63 are usable). */
private long high;
}
• @JS1 I guess I could figure this out... comparison of unsigned longs can be done with two "xors" and a single compare. Just undelete your answer. I can't see how you can throw away carries from p2, could you elaborate? – maaartinus Dec 12 '14 at 1:04
• I undeleted my answer now, with a much better answer that is faster than even my previous one. And it's even portable to Java now! – JS1 Dec 12 '14 at 21:07
• Btw if the only thing you are using this for is the Collatz sequence, you could probably write a very efficient "multiply by 3" function. – JS1 Dec 13 '14 at 12:50
• @JS1 I'll have a look at it soon -- I'm looking forward to it. Sure, I could multiply by 3 much simpler, but I don't need it! I need to multiply by powers of three, namely $\{1, 3, 9, \ldots, 3^{17}\}$. That's a very substantial optimization allowing me to make 17-34 steps at once. I'll submit the code for a CR later, now you can see it at github only. – maaartinus Dec 13 '14 at 14:43
• @JS1 This brought me to the idea to replace the multiplication by some hard-coded shifts. The exponent is known and many needed powers of 3 are rather sparse, e.g., 3, 9, 81=0x51, $3^7$=0x88B, $3^8$=0x19A1. – maaartinus Dec 13 '14 at 15:01

After many iterations, where I tried stuff in C that didn't port to Java cleanly, I came up with the following, that should work equally well in Java as in C. When tested in C versus the OP's function (ported to C), this function was about 60% faster. I don't know if that speed will translate to Java or not. Here is the Java version:

public void multiply(long x)
{
long y0  = unsignedLow (low);
long y1  = unsignedHigh(low);
long x0  = unsignedLow (x);
long x1  = unsignedHigh(x);
long y12 = composeLowHigh(y1, high);

// The low 64-bits of the output can be computed easily.
low = x * low;

// The upper 64 bits of the output is a combination of several
// factors.  These are the first two:
high = (high * x0) + (y12 * x1);

// Now handle the factors coming from the rest:
long p01 = x0 * y1;
long p10 = x1 * y0;
long p00 = x0 * y0;

// Add the high parts directly in.
high += unsignedHigh(p01);
high += unsignedHigh(p10);

// Account for the possible carry from the low parts.
long p2 = unsignedHigh(p00) + unsignedLow(p01) + unsignedLow(p10);
high += unsignedHigh(p2);
}

### An explanation of the algorithm

To explain the algorithm I use, first it would be good to look at how to do the multiplication the "normal way", which is the way the OP's function works. The normal way involves using 32-bit multiplies to create 64-bit partial products, and adding the partial products together.

Multiplication of 64-bit x against 128-bit y using 32-bit multiplies:

U(n) means upper 32 bits of n.
L(n) means lower 32 bits of n.

x0 = L(x)
x1 = U(x)
y0 = L(y.low)
y1 = U(y.low)
y2 = L(y.high)
y3 = U(y.high)

p00 = x0 * y0
p01 = x0 * y1
p02 = x0 * y2
p03 = x0 * y3
p10 = x1 * y0
p11 = x1 * y1
p12 = x1 * y2
p13 = x1 * y3

(              Bits in the result                 )
0.......32......64......96......128.....160.....192
[------p00------]
[------p01------]
[------p10------]
[------p02------]
[------p11------]
[------p03------]
[------p12------]
[------p13------]

First notice that p13 is not needed because we are not looking for any bits beyond bit 127 of the result. That leaves 7 multiplies (notice the OP's function has exactly 7 multiplies), and a bunch of adds/shifts. The formula for the resulting low and high could be written as:

low  = p00 + L(p01)<<32 + L(p10)<<32
high = U(p01) + U(p10) + p02 + p11 + L(p03)<<32 + L(p12)<<32 + Carry
Carry = U(U(p00) + L(p01) + L(p10))

This algorithm uses this idea to speed up the computation:

We can use a 64 bit x 64 bit multiplication when we only need the bottom 64 bits of the product.

What does this mean? There are several places where we can save multiplies and adds by multiplying more than 32-bits at a time. What happens when we multiply 64 bits against 64 bits? Suppose we define these 64 bit values:

x01 = x
y01 = y.low
y12 = U(y.low) + L(y.high)<<32
y23 = y.high

What does multiplying x01 * y01 produce? The full 128-bit result in terms of the partial products listed above would be:

x01 * y01 = p00 + p01<<32 + p10<<32 + p11<<64

But since the result is only 64 bits, the actual 64-bit result would be:

x01 * y01 = p00 + L(p01)<<32 + L(p10)<<32

But if you notice, this is exactly what we needed to compute for low. So the first computation of my algorithm is exactly that:

low = x * low;

This one line solves half the problem with a single multiplication. The rest of the function is dedicated to finding the upper 64 bits. Let's look at how the equation for the upper 64 bits can be simplified:

high = U(p01) + U(p10) + p02 + p11 + L(p03)<<32 + L(p12)<<32 + Carry

We can use the same idea to combine some of the multiplies together. Except instead of multiplying two 64-bit numbers together, we can multiply a 32-bit number by a 64-bit number. The value p02 + L(p03)<<32 is simply x0 * y23. Similarly, p11 + L(p12)<<32 is x1 * y12. So this has been reduced to:

high = U(p01) + U(p10) + x0*y23 + x1*y12 + Carry
Carry = U(U(p00) + L(p01) + L(p10))

This is exactly what is computed in my function. It would be nice if there was a better way of taking care of the "Carry", but as of right now I'm not sure if there is anything better than the straightforward way. I did try a different way but it ended up being slower so I reverted that change.

In the end, there are 6 multiplies which is only one less than the original function, but there are only 6 adds which is significantly less the original. There are also 6 truncate operations where we take the upper or lower 32-bits, and that is also a lot less than in the original function.