The Knuth shuffle

The Knuth shuffle algorithm is an algorithm for randomly shuffling the elements of an array arr of size n given a function that generates a uniformly distributed random number in \mathcal{O}(1) time. The algorithm produces an unbiased permutation. That is, every permutation is equally likely.

Proof that the algorithm is unbiased
For any element e, the probability that it will be shuffled into position n-1 (last position at zero based indexing) is: the probability that it gets selected when i=n-1, which is \frac{1}{n}.

For any element e, the probability that it will be shuffled into position n-2 (second-to-last position) when i=n-2 is: the probability that it gets selected when i=n-2, which is \frac{1}{n-1} multiplied by the probability that it is not selected for the last position, which is 1 - \frac{1}{n} = \frac{n-1}{n}, giving  \frac{1}{n-1}\times\frac{n-1}{n} = \frac{1}{n}

For any element e, the probability that it will be shuffled into position n-3 (third-to-last position) when i=n-3 is: the probability that it gets selected when i=n-3, which is \frac{1}{n-2} multiplied by the probability that it is not selected for the last nor second-to-last position, which is 1 - \frac{2}{n} = \frac{n-2}{n}, giving \frac{1}{n-2}\times\frac{n-2}{n} = \frac{1}{n}.

This can easily be generalised for any other position.

So we can see that for any element e, the probability that it will be shuffled into the last, second-to-last, third-to-last, etc position is \frac{1}{n}.

 

Implement a ceil() function

In this post I will explain a less intuitive approach of implementing a ceil() function in Java, that is faster than known comparison approaches.

A few known facts
Fact 1
Dividing two integers a and b in Java will always floor the result:

well almost. Actually Math.floor(a/b) rounds towards negative infinity, where as a/b rounds towards zero:

Fact 2
\lfloor x \rfloor = m if m \le x \textless m+1
\lceil x \rceil = n if n - 1 \textless x \le n

Fact 3
\lfloor x \rfloor \le \lceil x \rceil

What we need is \lfloor (x + s) \rfloor == \lceil x \rceil
Since dividing two integers a and b in Java will always floor the result, we need to find some s such that \lfloor (x + s) \rfloor == \lceil x \rceil.

How do we find s?
Think about it this way: \lceil \frac{a}{b} \rceil \textgreater \lfloor \frac{a}{b} \rfloor, except when \frac{a}{b} is a whole number. So, we want to bump \frac{a}{b} to (or past) the next whole number, unless \frac{a}{b} is a whole number, by adding s. This way \lfloor \frac{a}{b} + s \rfloor == \lceil \frac{a}{b} \rceil. We want to find s so that for whole numbers, it won’t quite push them up to the next whole number, but for non-whole numbers it will.

For instance, assume a=6, b=3 and s=1 then \lfloor \frac{a}{b} + s \rfloor = \lfloor \frac{6}{3} + 1 \rfloor = 3 \ne \lceil \frac{a}{b} \rceil = \lceil \frac{6}{3} \rceil = 2, which is not the desired result. As you can see, by using s=1 we bumped \frac{a}{b} to far. So s must definitely be less than 1.

So how much s is just enough?

Assume a=7, b=3. Then \frac{a}{b} = \frac{7}{3}, which can be re-written as:
\frac{7}{3} = \frac{1}{b} + \frac{1}{b} + \frac{1}{b} + \frac{1}{b} + \frac{1}{b} + \frac{1}{b} + \frac{1}{b} = (\frac{1}{3} + \frac{1}{3} + \frac{1}{3}) + (\frac{1}{3} + \frac{1}{3} + \frac{1}{3}) + \frac{1}{3} = 2  + \frac{1}{3}
In other words, we get an integer part (2) and a left over (\frac{1}{3}). The left over part is always in the form of k\frac{1}{b}, with 0 \le k \textless b. Thus, the smallest non-whole number we can get from \frac{a}{b} is qb + \frac{1}{b} with qb being the integer part with q \in \mathbb{N}_0.

So assuming \frac{a}{b} is not a whole number, the smallest left over we can have is \frac{1}{b}. So in-order to bump \frac{a}{b} to the next whole number, it suffices to add s = 1 - \frac{1}{b}. This is \textless 1, which will not bump \frac{a}{b} to the next whole number in case \frac{a}{b} is a whole number, but still enough to bump \frac{a}{b} to the next whole number in case \frac{a}{b} is not a whole number.

Summing it up
So we know adding s = 1 - \frac{1}{b} to \frac{a}{b} will fulfill the equality: \lfloor (\frac{a}{b} + s) \rfloor = \lceil \frac{a}{b} \rceil. Thus we get:

 \lceil \frac{a}{b} \rceil = \lfloor (\frac{a}{b} + s) \rfloor = \lfloor (\frac{a}{b} + 1 - \frac{1}{b}) \rfloor = \lfloor \frac{a + b - 1}{b} \rfloor

 

Exponentiation by squaring

This is a rather short post. in which I explain how to we can implement a fast pow function to raise a number x to the power of n.

A rather naive approach would be to do the following:

This is ok, but the algorithm has a running time of \mathcal{O}(n). We can do better than that.

Let us consider an example.

Assume that we want to calculate x^8. We can rewrite it as follows:

x^8 = x \times x \times x \times x \times x \times x \times x \times x.

which requires 7 multiplications. What we can do to reduce the number of multiplications is to first compute:

x \times x = x^2

next we can compute

x^2 \times x^2 = x^4

and last

x^4 \times x^4 = x^8

which gives as a total of 3 multiplications.

Lets assume that we want to calculate x^{13}. We can rewrite it as follows:

x^{13} = x^8 \times x^4 \times x.

with 8, 4 and 1 being all powers of 2. The great thing about 8, 4 and 1 being a power of 2 is that we can easily get them by squaring (see previous example).

How do we find the numbers 8, 4 and 1? Or in other words, how can we find out what powers of 2 make up 13?
Well that’s what the binary representation of 13 tells us:

13_d = 1101_b

Now that we have our powers of 2, all we need to do is to keep squaring (that is compute x \times x, then (x \times x) \times (x \times x) , etc) and every time we reach a power of two (‘1’ in binary representation) we multiply the intermediate results. So in order to compute x^{13} we iterate over the binary representation of 13 from the least significant bit (right-most bit) to left-most and perform our squaring and multiplication.

Algorithm

  1. Set result = 1
  2. 1101  \rightarrow we have a power of 2 (2^0), so we multiply result with x: result = result \times x
  3. Square x: x = x \times x  \rightarrow x = x^2
  4. 1101  \rightarrow we don’t have a power of 2 so skip.
  5. Square x: x = x \times x  \rightarrow x = x^4
  6. 1101  \rightarrow we have a power of 2 (2^2), so we multiply result with x: result = result \times x
  7. 1101  \rightarrow we have a power of 2 (2^3), so we multiply result with x: result = result \times x

Code

The running time thus becomes \mathcal{O}(\log n)