# February 16, 2015

## Swiftly Secure, Part 3 – Dividing Huge Numbers

Posted in iPhone development, mac development tagged , , , , , at 9:05 pm by tetontech

Let’s face reality. Computer division is ugly. It’s not the case, but it might seem to some people that those who conceived the computers we use today dropped the ball. Addition and Multiplication algorithms at both code and hardware levels are fairly straight forward, even if there is significant computational cost involved. Division, on the other hand, can seem like it was ‘tacked on’ afterwards. To think this thought is to do a disservice to the computational founders. The problem didn’t arise with them. It came from how we humans have done division for millennia. The age-old long hand grade school division algorithm taught to young people around the world is the problem. It is a guess-and-check process and computers are notoriously bad at guessing. It is no wonder the founders struggled. What else did they have to base division on?

When I set out to find a division algorithm to implement for Swiftly Secure and this posting, the other common division method, inverting the divisor, jumped right out of old and new computer science patents and research. The inversion idea is that a / b is the same as a * 1/b and it may be easier to calculate 1/b than to do the original division. Modern research into division seems to be fixated on optimizing variations on this approach. Each variation seems to share the same weakness as the grade school algorithm–they end up being guess-and-check algorithms trying to outdo Isaac Newton’s division algorithm or its more modern variants and replacements. This doesn’t resolve the problem, it shifts it to an area where guessing is hopefully easier.

I’m going to make a bold statement here. Computers are bad at guessing. I wanted to come up with an algorithm where no guessing was going on to use in Swiftly Secure. I think I’ve found one. It works for integer division for single word integers and huge numbers that require multiple words to represent them. I haven’t tested the algorithm for applicability to floating point division yet.

The idea behind the solution is this–since computers work in binary maybe there is something in the binary representations of the dividend, a, the divisor, b, and the quotient, q (a / b = q) that allows calculation of the quotient without guessing. There is. Finding it was a significant process and like all good answers to interesting questions, after I found it it the solution seemed obvious. Here it is simplified.

The difference between the index of the highest on-bit of a (the dividend) and b (the divisor) is often the index of the highest on-bit of q (the quotient).

A) if a / b = q

B) then  i[a] – i[b] = i[q] (where i means the index of the highest on-bit)

Where statement B isn’t true, the difference calculated is always one too large. When the difference is off is completely predictable, calculable, and derivable from  the value of the dividend (I’ll show how to do this later in this posting).

Interesting, but h0w does this help us not guess and check? Instead of guessing and checking by attacking the dividend and/or the divisor and then seeing if the quotient is ‘close enough’, as is done by Newton’s method and its decedents, this realization allows us to attack the quotient, q, directly. As an example lets look at a simple one-word division.

127 / 17 = 7 Remainder 8

decimal                 binary (bigendian)              highest on-bit index (zero based)

127                        01111111                                               6

17                         00010001                                            4

7                          00000111                                              2

i = 6,

i = 4, and

i = 2

6-4 = 2

The answer’s highest on-bit index is calculable. This implies an algorithm that reduces the dividend while increasing the quotient.

1) calculate the quotient’ highest on-bit index.

2) calculate the numeric value of the quotient’s highest on-bit index

3) add the numeric value calculated to a running quotient accumulator

4) subtract from the dividend the numeric value * the divisor

5) loop back to 1) until the updated dividend is less than the divisor

6) the final value of the dividend is the division’s remainder

No guessing is involved nor is there a ‘stop when close enough’ guess as in Newton’s and many other divisor inversion methods. The code for this algorithm, even for huge multi-word unsigned integers, lays out easily in Swift.

```func divide(var dividend:[UInt64], divisor:[UInt64])
->([UInt64],[UInt64],NSError?){
var quotient = [UInt64]()
/*
* sanity checking for division by zero, empty arrays, and
* where the divisor > dividend goes here
*/
//5) loop until the divisor > dividend
while compare(dividend, divisor) >= 0{
//1) calculate the highest on-bit index
let highestQuotientBit = calculateHighBitIndex(dividend, divisor)
//2) calculate the numeric value of the highest on-bit
let dividendPartialValue = shift_1_LeftBy(highestQuotientBit)
//3) update the value of the quotient collector
//4 subtract the calculated numeric value from the dividend
let tooSmallBy = multiply(dividendPartialValue, divisor)
dividend = subtract(dividend, tooSmallBy)
}
//the amount left in the dividend is the remainder
}
```

In this algorithm I take advantage of the multiplication algorithm discussed and implemented in a previous post. There are several additional helper functions I’ve created to make the divide function’s implementation easier to understand.

The most interesting of these helper functions is calculateHighBitIndex. It calculates the index of the highest on-bit of the quotient for any divisor-dividend pair. To help you understand it, I’ll show the equations for calculating when the difference between the dividend and divisor’s highest on-bit index, described above, is too high by 1.

when a / b = q

if  a <= (b *  (1 << (i[a] – i[b]))) – 1

then i[a] – i[b] is too large by 1

This is a pattern I found by examining every integer division possible using unsigned 8 bit integers. I’ve been running an exhaustive check of all possible UInt64 divisions. It hasn’t failed yet but is obviously taking a long time. I’m also working on a mathematical proof but won’t bore you with it here.

Describing the equations using words, if 1 is left-shifted by the highest on-bit index difference, the result is then multiplied by the dividend, and reduced by 1, we get a value for the largest dividend where the highest on-bit index difference is too large. The source code is as follows.

```//calculate the highest on-bit of a multi-word quotient
//given multi-word dividends and divisors
func calculateHighBitIndex(var dividend:[UInt64], var divisor:[UInt64])
->UInt64{
var calculatedIndex:UInt64 = 0
if dividend != divisor{
//cleanup the dividend and the divisor before processing
//by removing any leading zero valued array elements
//dividend's highest bit location ( log2 dividend)
//for 64 bit array elements
let dividendLog2 = UInt64(floorLogBase2(dividend)
+ (dividend.count - 1 ) << 6)
//divisor's highest bit location (log2 divisor)
//for 64 bit array elements
let divisorLog2 = UInt64(floorLogBase2(divisor)
+ (divisor.count - 1) << 6)
//calculate the inter-highest on-bit index
calculatedIndex = dividendLog2 - divisorLog2
//calculate the maximum dividend value yielding the error
var maxLocation:[UInt64] = subtract(multiply(divisor,
shift_1_LeftBy(calculatedIndex)),
[UInt64(1)])
//check if the calculated index needs to be reduced by 1
if compare (dividend, maxLocation) <= 0{
calculatedIndex--
}
}
return calculatedIndex
}
```

I created the floorLogBase2 function because the built-in log function works on doubles, is imprecise, and there is an easy, efficient algorithm to calculate log base 2 for unsigned ints. The code for this function is

```func floorLogBase2(var x:UInt64) ->Int{
//calculate the numeric value of the highest bit (1,2,4,8,etc.)
x |= (x >> UInt64(1))
x |= (x >> UInt64(2))
x |= (x >> UInt64(4))
x |= (x >> UInt64(8))
x |= (x >> UInt64(16))
x |= (x >> UInt64(32))
let numericValueOfHighestBit = (x - (x >> 1))

//lookup log base 2 of the value
//it is faster to look up the log than calculate it
return bitLocations[numericValueOfHighestBit]!
}
```

The last line of the floorLogBase2 function uses a Dictionary to look up the logarithm. The Dictionary’s key is 2^n and the value is n.

The other helper functions (compare, add, subtract, shift_1_leftBy, and removeLeadingZeros) aren’t that interesting. The source code for them will be part of the open source Swiftly Secure library I’ll be putting on GitHub when its complete.

When I started looking for a way to divide huge multi-word numbers, I wanted to find a way to do computerized division without the guesses required by other modern methods. I found one and a way to implement it in Swift. As I wrote the Swift code, I kept optimal computational practice in mind, but I make no claims regarding how ‘optimal’ this algorithm is compared to others. I plan on doing theoretical (big O) and actual comparisons a little later. First I have to find a way to efficiently raise huge multi-word numbers to huge multi-word powers without repeatedly multiplying the base number by itself. This is the last missing piece to completing Swiftly Secure’s implementation of the socialist millionaire protocol.

[UInt64]^[UInt64] is my next task and posting.

# February 10, 2015

## Swiftly Secure, Part 2 – Multiplying Huge Numbers

Posted in iPhone development, mac development tagged , , , , , , , at 12:11 am by tetontech

In my previous post I mentioned how two 64 bit numbers could be multiplied together without ever overflowing. The approach taken was a variation of the grade school multiplication algorithm. In this posting I continue my path to an implementation of the Socialist Millionaire protocol used to securely communicate between 2 or more peers. In order use socialist millionaire I will have to multiply numbers that are much too large to fit into a 64 bit integer. This means I’ll need to multiply arrays of numbers by arrays of numbers.

There are algorithms out there for doing multi-word multiplication. Warren, at hackers delight, lays out a C version of Knuth’s grade school algorithm. It works, but has a design point I was unwilling to follow. It requires the arrays being multiplied to consist of half-word numbers rather than full word numbers. I wanted to explore code using arrays of unsigned 64 bit numbers to represent the very long numbers. It seemed a waste to only use 32 bits per number.

With that in mind I sat down to figure out how to apply the multiplication algorithm used in the last post in this new situation. I wasn’t surprised when I found was easily adaptable from multiplying two Swift UInt64’s to multiplying two Swift UInt64 arrays (Please see the graphics and discussion in the previous post to understand the algorithm basics).

Part of the new algorithm implementation includes adding two UInt64’s in situations where they could readily overflow. That meant I needed to create an add method to handle this correctly. It isn’t that complicated but does require the use of bitwise operators and comparisons.

```func add(leftAddend:UInt64, rightAddend:UInt64)
->(UInt64,UInt64){
var carryAmount:UInt64 = 0
let sum:UInt64 = 0
//check for overflow
carryAmount = 1
}
else{
}
return (carryAmount,sum)
}
```

This implementation of add checks to see if overflow will happen before doing the addition. If there is no danger of overflow the two addends are combined. If overflow will happen the amount of overflow is stored and the sum is set to the maximum value of a UInt64.

In either case the add function returns the sum and a number representing how much too large the real sum is if it was to fit into a UInt64. We’ll use this carry amount in the multi-word multiplication algorithm.

Another change I made from Warren’s implementation was to have the arrays of numbers be in big endian format. This didn’t add significantly to the complicatedness of the algorithm but is how people ‘normally’ think of numbers. This can make it easier for some to see how the input relates to the algorithm’s output.

If the algorithm is implemented correctly some simple multiplications and their results should match these below. All arrays of numbers are base 2^64 where each array element represents a 64 bit ‘digit’ in big endian bit order.

• multiply([UInt64.max/2 + 1],)     => [1,0]
•  multiply([1,0,0],[1,1])                          => [0,1,1,0,0]
• multiply([2,2],[UInt64.max/2 + 1])  => [1,1,0]
• multiply([UInt64.max, UInt64.max, UInt64.max], [UInt64.max, UInt64.max])          => [18446744073709551615, 18446744073709551615, 18446744073709551615, 18446744073709551615, 1]

In this implementation of the grade school algorithm discussed in the previous post the ‘bit’ size of the multiplicand and the multiplier are unknown. This means we can’t use the nice trick with the << operator and bit masks to create portions of the product we can add together. Instead we’ll collect the product portions directly into an array of UInt64’s.

If we create an array for the resultant product sized to be the sum of the sizes of the multiplicand and the multiplier it will always be large enough to hold the product. Sometimes it will be 1 digit too large but that will always be in the most significant digit location when it does happen. You can see an example of this if you look at the second bullet point of the list above. This is OK since a zero value in that place doesn’t affect the value of the number anyway.

The multi-word multiplication implementation of the algorithm looks like this.

```func multiply(multiplicand:[UInt64], multiplier:[UInt64])
->[UInt64]{
let multiplicandSize  = multiplicand.count
let multiplierSize = multiplier.count
let productSize = multiplicandSize + multiplierSize
//create an array that is large enough
//to hold the product
var product = [UInt64](count: productSize,
repeatedValue: 0)
//to do multiplication
for index in stride(from: multiplierSize - 1,
to: -1, by: -1){
let multiplierElement = multiplier[index]
for subIndex in stride(from: multiplicandSize - 1,
to: -1, by: -1){
let multiplicandElement = multiplicand[subIndex]
//calculate the product and the carry amount
let partialProduct:[UInt64] =
multiplyWithOverflow(multiplicandElement, multiplierElement)
//calculate the location of the product
let productIndex = (subIndex + index + 1)
//calculate the index of the carry amount
let carryIndex = productIndex - 1

//calculate the amount to be inserted
//into the product and carry locations.
//Add the current calculated product to any amount that
//was previously carried forward.
partialProduct)
//the carrySum will never be able overflow.