Showing posts with label article. Show all posts
Showing posts with label article. Show all posts

Thursday, April 16, 2026

Primitive Pythagorean triples with the help of trigonometry

 Not a puzzle, just a quick note.

Primitive pythagorean triples are of the form $m^2 - n^2, 2mn, m^2 + n^2$ with $m,n$ relatively prime and of opposite parities.

We will use trigonometry to help derive this.

Let $a,b,c$ with $a^2 + b^2 = c^2$ be a primitive pythagorean triplet. This means $\gcd(a,b,c) = 1$. Note that this also implies that $a,b,c$ are pairwise coprime. Thus exactly one of $a,b,c$ is even. It cannot be $c$ as that would imply $4$ divides $a^2 + b^2$ for $a,b$ odd (which is not possible). 

Thus assume $b$ is even. 

We can now prove that 

$\textbf{Proposition:}$ With $b$ even, the highest power of $2$ that divides $b$ is different from the highest power of $2$ that divides $a+c$. 

$\textbf{Proof:}$ Since $a,c$ are odd and relatively prime, we must have that $\gcd(a+c, c-a) = 2$. Thus $a+c = 2u$ and $c-a = 2v$ for some odd $u,v$. Now if $b = 2w$ for some odd $w$ (if we assume that the highest power of $2$ that divides $b$ is same as $a+c$), then $b^2 = (c+a)(c-a) \implies w^2 = uv$. Since $w^2 = 1 \mod 4$, this implies $u = v$ mod $4$, and thus $a = u - v = 0 \mod 4$, not possible.

Now consider the triangle $$B = (0,0), C = (a, 0), A = (a,b)$$

Draw the anglular bisector $BD$ of angle $B$ such that $D$ lies on $AC$. By the angular bisector theorem, $CD/DA = a/c$ and hence $CD = \frac{ab}{a+c}$. Thus we have that

$$\tan \left(\frac{B}{2}\right) = \frac{CD}{BC} = \frac{b}{a+c} = \frac{n}{m}$$

For some $n, m$ relatively prime.  Since the highest power of $2$ that divides $b$ is different from that of $a+c$, $n,m$ are of opposite parities.

Now using $\tan 2x = \frac{2 \tan x}{1 - \tan^2 x}$ we get

$$\frac{b}{a} = \tan B = \frac{2 n/m}{1 - (n/m)^2} = \frac{2nm}{m^2 - n^2}$$

We can show that $2mn$ and $m^2 -n^2$ are relatively prime if $m,n$ are relatively prime and of opposite parities: if a prime $p > 2$ divides the numerator, it has to divide one of $m,n$ and cannot divide the denominator also. The denominator is also odd, so their $\gcd$ cannot be $2$.

Thus $$b = 2mn, a = m^2 - n^2, c = m^2 + n^2$$

Thursday, February 16, 2023

A proof of Pythagoras Theorem

Here is a potentially new proof of the Pythagoras theorem I discovered recently. I could not find any prior proofs like these (including in the cut-the-knot pythagorean proof dump), so if you have seen it somewhere, please let me know.





ABC is the right triangle with sides a,b,c.


AD is the angular bisector of CAB and DE is the perpendicular from D to AB.


CD = x, DB = a-x.

It is easily seen that triangles ACD and ADE are congruent and so AE = AC = b and BE = c-b


Area of ACD + Area of ADB = Area of ABC gives

$bx + cx = ab $

And so $x = ab/(b+c)$

Now triangles DEB and ABC are similar and so

$$DE/AC = BE/BC$$

i.e.

$$x/b = (c-b)/a$$


$$a/(b+c) =(c-b)/a$$

Which gives

$$ a^2 = c^2 - b^2$$




Thursday, February 9, 2023

Base b and generating functions

 Not a puzzle, just a quick observation noting the power of generating functions.


Say $b \ge 2$ is a positive integer. We give a quick proof using generating functions that every positive integer can be written uniquely in base $b$ with digits $0, 1, \dots, b-1$.


Let $n \ge 1$. Consider 

$$ G(x) = \prod_{k=0}^{n} (1 + x^{b^k} + x^{2b^k} + \dots + x^{(b-1)b^k})$$

Observe that the coefficient of $x^N$ shows the numbers of ways of writing $N$ is base $b$.

Now if $u = x^{b^k}$, then

$$ (1 + x^{b^k} + x^{2b^k} + \dots + x^{(b-1)b^k}) = (1 + u + u^2 + \dots u^{b-1})  = \frac {1 - u^b}{1-u}$$


Thus

$$G(x) = \frac{1-x^b}{1-x}\frac{1-x^{b^2}}{1-x^b}\frac{1-x^{b^3}}{1-x^{b^2}}\dots\frac{1-x^{b^{n+1}}}{1-x^{b^n}}$$


This telescopes to


$$ G(x) = \frac{1 - x^{b^{n+1}}}{1-x}$$

$$ = 1 + x + x^2 + \dots + x^{b^{n+1} - 1}$$

Thus every integer $N$ in the range $1$ to $b^{n+1} - 1$, can be written in base $b$ using at-most $n+1$ digits, taken from $0,1,2, \dots, b-1$. Since the coefficient of $x^N$ is $1$, the representation is unique.


Saturday, November 27, 2021

4 people in a square with hats

I was given this puzzle by Arvind Hariharan and he has a write up with a different solution approach here: https://www.starvind.com/math/hats-to-leave-you-scratching-your-head/.

The puzzle (apparently from Riddler) is as follows:

Four people are trying to escape from a room. Guards have placed a hat on each person’s head, and each hat is one of three colors: red, yellow or blue. The four people are arranged at the vertices of a square, with an obstacle in the middle. Each person can see the hats on the heads of those on adjacent vertices of the square, but they cannot see the hat of the person diagonally across from them. They also do not know the color of the hat on their own head.

Each person must guess the color of the hat on their own head. If at least one person guesses correctly, they can all escape the room together. No communication is allowed once the hats are placed on their heads, but they can coordinate on a strategy beforehand. They also know how they will be arranged in the square.

How can they be guaranteed to escape the room?


Please stop reading if you wish to solve the puzzle by yourself.


The reason for writing this article is to discuss a strange approach to solving this problem: Giving a proof of existence of a strategy using Linear Algebra!


Since there are $3$ colours, we can work mod $3$, with $0, 1, 2$ being the colours and $x + y = x + y \mod 3$ (This is the field $\mathbb{F}_3$). Note that $2x = -x$ etc.

Let us assume the four people are $A,B,C,D$ and they are given the hat colours $a, b, c, d$, with $A,C$ and $B,D$ being diagonally opposite. Thus $A$ and $C$ only have access to $b,d$ and $B, D$ only have access to $a, c$.

Let us further assume that there is strategy where each person computes some linear combination of the colours they can see and uses it as their guess. 

For eg: $A$ guesses $x_a b + y_a d$ for some $x_a, y_b$ specific to $A$.

Now introduce the offset variable $\delta_a$ which is the amount $A$'s guess is off from their actual value. i.e.

$$ \delta_a = x_a b+ y_a d - a = x_a b + y_a d + 2a$$

Writing this equation for each of $A,B,C,D$ we see that the $\delta$'s can be captured as a matrix equation as follows

$$S \begin{bmatrix}a\\ b \\ c \\ d\end{bmatrix} =  \begin{bmatrix}\delta_a \\ \delta_b \\ \delta_c \\ \delta_d \end{bmatrix} $$

where $S$ (short for strategy) is the 4x4 matrix 

$$   \begin{bmatrix} 2 & x_a & 0 & y_a\\ x_b & 2 & y_b & 0\\ 0 & x_c & 2 & y_c\\x_d & 0 & y_d & 2 \end{bmatrix}$$

$S$ is a viable strategy if for every possible $a,b, c, d$, the atleast one of the $\delta$'s is $0$.

Now comes the key observation:

Lemma 1: For any $x, y$ at least one of $x, y, x + y, x+2y$ is $0$ (note $\mathbb{F}_3$ working mod 3). This is easily proved.

This fact can be stated in linear algebra terms

Lemma 2: Each vector in the space spanned by $u = \begin{bmatrix} 1 & 0 & 1 & 1 \end{bmatrix}$ and $v = \begin{bmatrix}0 & 1 & 1 & 2 \end{bmatrix}$ has at least one co-ordinate which is $0$, i.e. for any $x, y$, the vector $xu + yv$ has at least one co-ordinate which is $0$.

Lemma 3: Now notice that if you now have two arbitrary values $p$ and $q$ and you wanted a 1x4 vector with one of the coordinates as $p$ and other as $q$ (in some given positions), you can find (a unique) vector of the form $w = x u + y v$ which has co-ordinates $p$ and $q$ in the given positions (the other two coordinates are not "free"), as you are essentially solving $x + y = p, x + 2y = q$ or some such pair of equations from Lemma 1.

Now look at the columns of $S$

$$   \begin{bmatrix} 2 & x_a & 0 & y_a\\ x_b & 2 & y_b & 0\\ 0 & x_c & 2 & y_c\\x_d & 0 & y_d & 2 \end{bmatrix}$$

Each column has $0$ and $2$ in some positions, thus by Lemma 3, there is an $S$ where each column of $S$ can be written as a vector of the form $x u + y v$.

Thus we see that for such an $S$, the offset vector $\begin{bmatrix}\delta_a \\ \delta_b \\ \delta_c \\ \delta_d \end{bmatrix}$, which is essentially a linear combination of the columns of $S$ is thus also a linear combination of $u$ and $v$ and by Lemma 2, has to have at least one co-ordinate which is $0$. 

Hence such an $S$ is a strategy which solves the puzzle. This $S$ can actually be computed easily, which we will leave to the reader.

Note that this generalizes, for example you could have a strategy for 6 people with 5 colours where each person could see all but one.

Sunday, July 19, 2020

From fair to irrational and back

This article is kind of related to the previous post about random numbers.

The main questions we want to answer:

Can we simulate a biased coin with a fair coin? What about the reverse?

We are more interested in the case where the biased coin needs to have an irrational probability.

We will consider the classic and a not so classic approach.

(We will also explicitly look at $\frac{1}{\sqrt{2}}$ and $\frac{1}{e}$ in this post).

From fair to irrational

So suppose we have a random bit generator which generates bit 0 (and 1) with probability $\frac{1}{2}$. Also assume that each call takes constant time.

Can we use this to simulate a bit generator that generates bit 0 with probability $p$ where $p$ is some irrational in $(0, 1)$? 

One condition we will require is that the biased bit generator will run in expected constant time, without that it is not of much use.

For motivation for the method: Say we had a way to generate a real number $r$ in $(0, 1)$ and compare that against $p$, we could then have an easy algorithm to return 0 when $r \lt p$ and return 1 when $r \gt p$. This will return $0$ with probability $p$.

The trick is to have a representation of $p$ and that of $r$ which can be compared while we are in the process of generating $r$. The most obvious representation that comes to mind is the base-2 representation!

Base-2 generator


So say we had a base-2 representation of $p = 0.d_1 d_2 d_3 \dots $

To generate $r = 0.r_1 r_2 r_3 \dots $ we can randomly generate the bits $r_1, r_2, \dots$ and compare against $d_1, d_2, \dots$ while we are generating the $r_i$. With a very high likelihood there will be some finite $k$ where $r_k \neq d_k$ for the first time. We can stop and return bit 0 or bit 1 depending on whether $r_k \gt d_k$ (and so $r \gt p$) or $r_k \lt d_k$ (and so $r \lt p$).

This can be written in code as follows:


# Given bitstream (base 2 representation) generates random bit with
# pbt of 0 = number represented by bitstream.
def random_bit_base2(bias_bitstream, fair_random_bit):
  # We generate a random number and if < given number
  # return 0, else return 1.
  for bit in bias_bitstream():
    fair_bit = fair_random_bit()

    # the number represented by the fair_bit seq is different!
    if fair_bit != bit:

      # This basically is returning 0
      # if number formed by fair_bit seq < bitstream seq.
      # else returning 1.
      return fair_bit



The probability of generating a 0 is

$$\frac{d_1}{2} + \frac{1}{2}\times\frac{d_2}{2} + \dots = p$$

The expected number of random number calls is $\sum P(X \gt n) = \sum \frac{1}{2^n} = 2$. [See this math.se post]

The problem with this approach is generating the base-2 representation of the required bias $p$, but for practical purposes (depending on the application) we can choose to only work with the first "few" (eg 128) digits of $p$. Very likely $p$ can be approximated via infinite series etc and the digits placed in a lookup table.

For certain $p$ it is pretty easy to do though. For eg:


Problem 1: Can you write a simple algorithm to lazily generate the base-2 representation of $\frac{1}{\sqrt{2}}$?

Hint: If $x^2 \lt 2^{2n+1} \lt (x+1)^2$ then we only need to consider $2x$ or $2x+1$ for $2^{2n+3}$.


This base-2 version generator is in fact pretty well known (Googling will give similar results here).

There are other approaches like base-n instead of base-2. Choosing the base $n$ will depend on the bias $p$. For eg, with the existence of spigot algorithms for $\pi$, we could use base-16 for say $p = \frac{\pi}{16}$.

We will now discuss the not-so-classic approach.

Alternative to base-n: factorial base


If we assume that instead of a fair bit generator, we had a random integer generator which given input $K$ returned one of $\{0, 1, \dots, K-1\}$ with probability $\frac{1}{K}$ then for constants like $\frac{1}{e}$ we could do something other than base-n.

We can use a representation which uses factorials. For a given $r \in (0, 1)$ there is a representation of $r$ as

$$ r = \frac{a_2}{2!} + \frac{a_3}{3!} + \dots  = \sum_{n=2}^{\infty} \frac{a_n}{n!}$$

And $a_m \in \{0, 1, \dots , m - 1\}$.

This representation (say $(a_2, a_3, \dots)$) is finite only when $r$ is rational (we also need $a_k \neq k-1$ infinitely often, but we can ignore that for our purposes).

For eg.

$$ e - 2 = \sum_{n=2}^{\infty} \frac{1}{n!} = (1, 1, \dots )$$

It can also be shown that to compare $r = (a_2, a_3, \dots) $ against $p = (b_2, b_3, \dots)$ it is enough to find the first $k$ where $a_k \neq b_k$ and compare the "digits".

To generate $a_k$ we just generate a random integer in $\{0, 1, \dots, k-1\}$ using the random integer generator we have.

The code would look something like this:


# Given digits in factorial base
# i.e num = \sum_{n=2}^{\infty} d_n/n!
# Generate a random bit with probability num.
def random_bit_base_factorial(digit_stream, fair_random_int):
  # We generate a random number and if < given number
  # return 0, else return 1.
  base = 2

  # We generate a random digit (0, ..., base - 1)
  # and then increment the base.
  for digit in digit_stream():
    fair_digit = fair_random_int(base)
    base = base + 1
    if fair_digit < digit:
      return 0
    if fair_digit > digit:
      return 1

Say we need the probability $p = \frac{b_2}{2!} + \frac{b_3}{3!} + \dots $.

The probability of returning 0 is

$$ \frac{b_2}{2} + \frac{1}{2} \times \frac{b_3}{3} + \frac{1}{2} \times \frac{1}{3} \times \frac{b_4}{4} + \dots = \frac{b_2}{2!} + \frac{b_3}{3!} + \dots = p$$

The expected number of fair random number calls is $\sum \frac{1}{n!} = e - 1 \lt 2$.

[It might look like this approach is better than base-2 because it converges faster and  uses fewer random number calls but we need to consider the implementation of the random integer generator, which might be using the random bit generator under the covers.]

This representation has the same problem as generating the digits in base-n of the required probability. For certain constants it is easy to do, though.

For eg:

Problem 2: Show that $$ \frac{1}{e} = \frac{2}{3!} + \frac{4}{5!} + \dots = \sum_{n=1}^{\infty} \frac{2n}{(2n+1)!}$$


A note regarding the "actual" runtime (applies to both base-n and factorial)


We have conveniently ignored the time taken to generate both the digits of the bias $p$ and the random digit of the generated real number $r$ but that is not a huge problem.

In fact we can prove the following:

Problem 3: If the $n^{th}$ digit generation of the bias $p$ and random number $r$ takes polynomial time in $n$, show that the expected runtime of the biased random bit generator is $\mathcal{O}(1)$.

Hint: Use the generating functions and their derivatives.

So if we hardcode the first say 128 digits in a lookup table and then switch to generating the rest, we will get a very good practical algorithm.

From Irrational back to Fair


This is actually not too hard, since if we generate two bits, the probability of getting 0,1 is the same as the probability of getting 1,0. 

So the algorithm is simple, call the biased bit generator twice. return 0 if we get 0,1. return 1 if we get 1,0. Continue if we get 0,0 or 1,1. This has finite expected runtime too.



Wednesday, July 15, 2020

Using math as an interviewer.

One of my go to interview questions is

"Generate a random binary tree".

We are only interested in structure of the tree and number of nodes will be < 1000.

I say that the ideal case is all trees are generated with the same probability (uniformly random), but explain that any algorithm where the probability of every tree is non-zero will do. Two trees are distinct if they are structurally different.

[Generating a random binary tree where each tree is equally like is not an easy problem (definitely too hard for an interview, though I have had a couple of candidate who did this!) and if you want to try it here is an old post about it: https://ruffnsluff.blogspot.com/2015/03/generate-random-binary-tree.html]


I make sure to tell them I don't need uniform (all equal probability) because it is a hard problem.

 Inspite of the explicit mention that it is a hard problem and not what I am looking for, candidates come up with all sorts of algorithms and claim it is uniformly random.

As an interviewer you need to be able to quickly tell if the algorithm is truly generating uniformly random or not.

I use a "trick" which has worked with almost every candidates' solution so far.

All I need to do is look at the random number generator calls their algorithm is making.

Typically the probabilities they generate are of the form $\frac{1}{n}$ or $\frac{1}{2}$ or $\frac{1}{3}$ or some such limited set like that.

The number of trees on $n$ nodes is the Catalan number: $ \frac{\binom{2n}{n}}{n+1}$ which tends to have many weird prime number factors (see, Bertrand's postulate). For eg, for $n=3$ we expect the probability of each tree to be $\frac{1}{5}$. For $n=4$ the probability of each tree needs to be $\frac{1}{14}$ (7 is the weird prime here).

No matter how the algorithm multiplies and adds these probabilities from the limited set together, they will fail to get the factor of the weird prime in the the denominator and thus can never generate the probability we require! I don't have to walk through their complicated algorithm and compute a tree etc.

[Caveat: Here we are assuming that the algorithm always terminates in a finite number of steps, no matter how the random numbers are generated. Section about this added at the end.]

For eg, consider this algorithm to generate a random tree.

1. Add root to tree.
2. Pick an existing node in the tree randomly that has child space
  2.1 Randomly add a child node to it.
3. Continue till tree has n nodes.


Does this generate the trees uniformly?

The probabilities involved are $\frac{1}{k}$ where $k \le n$. No matter how you combine these you will not get the the weird prime factor (the one between $n$ and $2n$) in the denominator. So this cannot possibly result in a uniform distribution.

 This looking at the denominators of the random number calls has been pretty useful when trying to see if some algorithm/process generates the probabilities we need.


As another example, here is a typical flawed algorithm to shuffle an array randomly, so that each permutation is equally likely.


for (int i = 0; i < array.length; i++) {
  int r = random(array.length); // random position where ith card goes.
  swap(array, i, r); // Swap a[i] and a[r]
}

Why is this not uniform? The random number generator only generates probability $\frac{1}{n}$. The number of permutations is $n!$ and we can never achieve $\frac{1}{n!}$ no matter how we combine those $\frac{1}{n}$s (for $n \gt 2$).

Why do we require the algorithm always terminate?


Consider the following

do {
  int bit1 = random_bit();
  int bit2 = random_bit();
  int result = bit1 + 2*bit2;
  if (result < 3) return result; //0, 1 or 2.
} while (true);

This generates a number in $\{0, 1, 2\}$ with probability $\frac{1}{3}$, even though the random number calls only generates probabilities of $\frac{1}{2}$ (the random bit calls).

Why does the denominator thing not work? Because we now have an infinite series for the probability calculation and that could give any result (even irrational numbers! See newer blog post.).

The probability of getting a $0$ is

$$\frac{1}{4} + \frac{1}{4^2} + \frac{1}{4^3} + \dots  = \frac{1}{3}$$

This corresponds to getting a 0, or a 4 (so while loop continues) and then a 0 or a 4,4, 0 and so on.

Note that the algorithm not halting with a probability zero is different from the algorithm always halting!

Friday, October 24, 2014

Stairs, Fibonacci numbers and O(log n) algorithms.

Introduction

You probably have heard of Fibonacci numbers.

They are the numbers in the sequence $1,1,2,3,5,8,13,\dots$ where each number is the sum of previous two numbers. So the next number after $13$ is $8+13=21$, and then $13+21=34$ and so on.

A combinatorial way of looking at them is the following:

Suppose you wanted to go up stairs having n steps, going up 1 or 2 at time. How many ways can you do this?

For instance, if there were 3 steps, you could go up as follows

$1,1,1$
$1,2$
$2,1$

where $1,2$ means go up one step, then go up two steps.

That the total number is a Fibonacci number is easy to see. If the number for $n$ steps was $F_n$, then the first number of steps is either one or two, and thus $F_n = F_{n-1} + F_{n-2}$.

We get $F_0 = 1, F_1 = 1, F_2 = 2, F_3 = 3, F_4 = 5, \dots$.

This combinatorial approach is a useful way [among many other approaches, like the matrix based approach] of looking at Fibonacci numbers. For instance, try proving the following identity combinatorially:

$$ \binom{n}{0}F_0 + \binom{n}{1}F_1 + \dots + \binom{n}{n}F_n = F_{2n} $$

($\binom{n}{k}$ is the binomial coefficient: $n$ choose $k$)

We will use this combinatorial approach to demonstrate an $O(\log n)$ arithmetic operations algorithm to compute the $n^{th}$ Fibonacci number. [The matrix approach also gives an $O(\log n)$ algorithm]

Note, we are talking about $O(\log n)$ arithmetic operations, rather than $O(\log n)$ time. This is because the size of the Fibonacci numbers becomes substantial. Since $F_n \sim \varphi^n$, $\varphi$ being the golden ratio, $F_n$ will be $\Theta(n)$ bits in size.

Adding $F_n$ and $F_n$ will take $\Theta(n)$ time, and multiplying $F_n$ and $F_n$ would typically be slower than that, its run-time depending on the multiplication algorithm used.

[Please stop reading if you want to try coming up with an algorithm yourself]

An $O(\log n)$ algorithm

So, we want to compute $F_n$, the number of ways to climb a stair of $n$ steps, going up one or two at a time. Assume for the moment, that you don't know that $F_n$ is a Fibonacci number.

To get an $O(\log n)$ algorithm, an effective technique is to reduce the problem in half (e.g: binary search).

So, if we could compute $F_n$ in terms of $F_{n/2}$ (or similar) we would have an $O(\log n)$ algorithm.

So let us try to compute $F_{2n}$ in terms of $F_n$.

Of all the ways to go up $2n$ steps, each possibility will have us end up at step $n-1$ or $n$ at some intermediate stage.

So we count the total number by splitting into two: possibilities that go through step $n-1$ but not step $n$, and possibilities that go through step $n$.

We get the following identity:

$$F_{2n} = F_{n-1}^2 + F_{n}^2 $$

And similarly:

$$ F_{2n+1} = F_{n}F_{n-1} + F_{n+1}F_{n}$$

Thus we can compute $F_{2n}$ if we knew $F_{n-1}$ and $F_{n}$, and we could compute $F_{2n+1}$ if we knew $F_{n-1}, F_{n}$ and $F_{n+1}$.

A naive recursive implementation of these identities would give an algorithm which requires $\Omega(n)$ arithmetic operations.

A common way to speed up such recursive algorithms is to use memoization: cache the previously computed values and use those instead of recomputing them.

As is frequently the case, using memoization gives us an exponential speedup, leading to an $O(\log n)$ arithmetic operations algorithm!

To prove that, consider what values we need, to compute the four Fibonacci numbers $F_{2n+1}, F_{2n+2}, F_{2n+3}, F_{2n+4}$:

$$F_{2n+1}, F_{2n+2}, F_{2n+3}, F_{2n+4} \to F_{n-1}, F_{n}, F_{n+1}, F_{n+2}\quad (1)$$

And to compute $F_{2n}, F_{2n+1}, F_{2n+2}, F_{2n+3}$:

$$F_{2n}, F_{2n+1}, F_{2n+2}, F_{2n+3} \to F_{n-1}, F_{n}, F_{n+1}, F_{n+2}\quad(2)$$

Thus starting with $F_{2n}$ (or $F_{2n+1}$) we can see that ultimately, we need to only compute a total of $O(\log n)$ such $F_i$.

For example, starting with $F_{8n+4}$:

$$F_{8n+4} \to F_{4n+1}, F_{4n+2} \to F_{2n-1}, F_{2n}, F_{2n+1} \to F_{n-2}, F_{n-1}, F_{n}, F_{n+1}$$

To compute $F_{8n+4}$, we need to compute two Fibonacci numbers of half the size ($F_{4n+1}, F_{4n+2}$). To compute those two, we need to compute three more of half their size ($F_{2n-1}, F_{2n}, F_{2n+1}$), and then four more. Based on $(1)$ and $(2)$ above, the amount of new numbers added to compute then stays at four, and never increases. At each step we halve the size, and so the total number of numbers is $O(\log n)$.

The number of arithmetic operations associated with computing $F_i$ is $O(1)$, if we use a memoization table, and so the whole algorithm will be using $O(\log n)$ arithmetic operations.

We can get rid of the memoization table by trying to compute the tuple of 4 elements instead: $T_n = (F_{n}, F_{n+1}, F_{n+2}, F_{n+3})$ and noticing that $T_n$ can directly be computed using $T_{\lfloor n/2 \rfloor -1}$.

[If you look at the matrix approach, you can view it as computing a tuple of three elements]

By massaging the identities a bit more, and using $F_{n+2} = F_{n+1} + F_{n}$, we can in fact, reduce the size of the tuple to two.

Further optimizations can be done, to reduce the number of multiplications (as they are slower than additions for large sized integers). See the GNU gmp Fibonacci implementation which actually does that by using similar identities to what we had earlier.

Conclusion

We looked at a combinatorial way of looking at Fibonacci numbers, and used that to give an $O(\log n)$ arithmetic operations algorithm.

This algorithm naturally falls out by applying two basic algorithm techniques: divide and conquer (by halving the input size) and caching (using memoization on the recursive version) and does not require any 'major leaps' (like the matrix approach).

Monday, September 22, 2014

Iterators part three: A General Binary Tree Iterator

[This is part three in a series about iterators]

Introduction

If you have heard of binary trees, you probably know what in-order, pre-order and post-order traversals are. These traversals are usually defined in a recursive fashion.

For instance, in-order is defined: traverse left sub-tree, visit node, traverse right sub-tree.

These kind of definitions lead to easy recursive implementations.

A common interview question is to convert these recursive method to iterative methods, which can be solved in a general fashion by doing what the compilers do: simulate the recursion using a call stack.

A possible followup to that is to add the ability to do the traversal in a lazy fashion: provide a next interface, which will return the visited nodes one by one (with next being reasonably efficient). [One could use yield, if the language provides it, but we will assume that it is unavailable.]

If you search the web, you will probably see different iteration solutions, custom made for in-order, pre-order etc.

A General Lazy Iterator for Binary Trees

In this article we will discuss a general lazy iterator for binary trees, for traversals defined recursively. This includes not just in-order etc, but also crazy traversals like traverse left node, traverse right node, visit node, visit node, traverse left node.

These general traversals can be represented as a string consisting of the letters "L", "R" and "V".

For instance, in-order traversal is "LVR", L for traverse left sub-tree, V for visit node and R for traverse right sub-tree. Post-order is "LRV" and the crazy traversal defined above is "LRVVL".

So, given a traversal method in the form a traversal string (like "LVR" and which contains at least one V), and the tree, we would like to create a lazy iterator, which will allow us to visit the nodes of the given tree in the given traversal order.

If you were to write a recursive method for the traversal, the code might look something like this (python):

# traversal_method is a string like "LVLR".
def traverse(tree, traversal_method):
    if tree == None:
        return
    for method in traversal_method:
        if method == "L":
            traverse(tree.Left, traversal_method)
        if method == "R":
            traverse(tree.Right, traversal_method)
        if method == "V":
            Visit(tree.Value) 

This recursive method can be converted to an iterative method, by maintaining a stack of (tree,  traversals_left) pairs.

For example, if the traversal required is "LVR", we first push the (root, "LVR") pair. Once the left sub-tree is traversed because of the "L" (started by pushing the left child), we update that to (root, "VR"), we then visit the node, update that node to (root, "R"), then traverse the right sub-tree.

Once a node's traversals_left becomes empty, we pop it off the stack, pick up what is on the top, examine what needs to be done and continue.

This is what the iterative version might look like (python):


# traversal_method is a string like "LVLR".
def traverse(tree, traversal_method):
    if tree == None:
        return
    stack = [(tree, traversal_method)]
    while len(stack) > 0:
        node, traversal_left = stack[-1]
        del stack[-1] # pop
        if node == None:
            continue
        # in python L[-1] is the last element of the list.
        # and L[1:] is the list without the first element.
        if len(traversal_left[1:]) != 0
            # push.
            stack.append((node, traversal_left[1:]))
        if traversal_left[0] == "L":
            stack.append(node.Left, traversal_method)
        if traversal_left[0] == "R":
            stack.append(node.Right, traversal_method)
        if traversal_left[0] == "V":
            visit(node.Value)        

In languages that support yield, an easy way to make this lazy would be to replace visit(node.Value) above with yield node.Value.

Assuming we don't have such a language feature, we want to implement a next interface, which can give us the nodes one at a time, and in a reasonably efficient manner (both worst case and amortized per call).

[We could probably blindly simulate what the C# compiler/python interpreter does behind the scenes for implementing yield, but the result might turn out to be an unreadable mess and we will probably have to spend some effort making it readable, anyway.]

We implement next as two steps.

First step is move, which will do the iteration of pushing/popping nodes into/off the stack, till we get to a node which needs to be visited, resulting in the appropriate node being at the top of the stack.

The second step of next would be to update the top of the stack indicating that a visit was just performed (and returning the tree node to the caller).

We also maintain the invariant that only the move step will change the size of the stack.

A C++ implementation might look like: [Click here to expand/collapse]

And here is a sample usage of this iterator
// Determine if two nodes of a binary search tree
// sum to target.
bool has_sum(Tree *tree, int target) {
  if (tree == NULL) return false;

  Traverser ascending(tree, "LVR");
  Traverser descending(tree, "RVL");

  Tree *left = ascending.next();
  Tree *right = descending.next();

  while (left->Value <= right->Value;) {
    int sum = left->Value + right->Value;
    if (sum == target) {
      return true;
    }
    if (sum > target) {
      right = descending.next();
    } else {
      left = ascending.next();
    }
  }
  return false;
} 

Analysis of next

If the size of the traversal description string is T [O(1) for most useful cases], and the maximum depth of tree is H, then next would run in O(TH) time, and would use O(TH) space.

Even though each call to next is O(TH), it is O(T) amortized per node visited (averaging over the whole traversal), as each node and each edge of the tree is touched only O(T) times.

Conclusion

Even though any recursive implementation which returns a list of elements, can theoretically be modified into a lazy iterative version, explicitly doing it for a binary tree was an interesting exercise.

This approach can also be generalized to general trees, where instead of "L" and "R", we can use child numbers etc.

Sunday, September 14, 2014

Iterators part two: Using the Odometer

[This is part two in a series about iterators.]



Introduction

Like the odometer of a car, the odometer we consider is basically a counter, with one generalization: each position could cycle through a different number of "digits" (each position starts at digit 0).

[Digits in quotes, because we could have one position go from, say, 0 to 999, each number is considered a digit, in base 1000.]

The general odometer works just like the one in the car.

Every time you need to increment the odometer, you increment the digit in the rightmost position (going back to 0 if you are at the largest possible digit for that position). If the increment results in that position cycling back to 0, you increment the position to the left of it. If that results in cycling back of that position to 0, you move one position left and so on.

Basically, simulating hand addition with carries.

When all the digits reach the maximum value, we cannot increment anymore. An odometer can easily be implemented as an iterator, to provide a has_next, next interface with next incrementing the counter by one.

Here is how an implementation might look like (C++):
#include <vector>

// Represents a generalized odometer.
class Odometer {
public:
  // 0 is the right most position.
  Odometer(const std::vector<int>& max_digits) {
    _max_digits = max_digits;
    for (int i = 0; i < max_digits.size(); i++) {
      _odometer.push_back(0);
    }
    _overflow = false;
  }

  bool has_next() const {
    for (int i = 0; i < _odometer.size(); i++) {
      if (_odometer[i] != 0 ) {
        return true;
      }
    }
    // All zeroes. If we overflowed,
    // we are at the end.
    return false || !_overflow;
  }

  // returns the current, but increments before
  // returning.
  std::vector<int> next() {
    _overflow = true;
    std::vector<int> current = _odometer;
    for (int i = 0; i < _odometer.size(); i++) {
      _odometer[i] = (_odometer[i] + 1) % _max_digits[i];
      if (_odometer[i] != 0) {
        break;
      }
    }
    return current;
  }

private:
  bool _overflow;
  std::vector<int> _odometer;
  std::vector<int> _max_digits;
}; 
 
We will now discuss how we can do a iterator based generation of permutations using an odometer (and briefly, Cartesian products and power-set enumeration).

Generating Permutations

[If you want to try solving the problem of implementing a has_next, next type of iterator for permutations using an odometer, please stop reading]

While the method discussed below is not optimal (unlike Heap's method, [in fact, not even asymptotically optimal]), and generates permutations in a "weird" order (unlike Narayana Pandita's method of generating in lexicographic order), it is an interesting application of odometers.

Say we want to provide an iterator to generate all the permutations of $1,2, \dots, n$, with a has_next, next interface.

The idea is to use a suitable odometer as a black-box, and map each odometer state to a permutation uniquely.

Since there are $n!$ permutations, an odometer that suggests itself:

The odometer has $n-1$ digits, $d_1, d_2, \dots, d_{n-1}$ with $d_1$ being the left-most.

The maximum possible digit achievable by $d_1$ is $n-1$, and so total $n$ digits for $d_1$, which are $\{0,1,2, \dots, n-1\}$.

Maximum digit for $d_2$ is $n-2$ and so on (max for $d_i$ being $n-i$).

This odometer can count upto $n\times(n-1)\times\dots\times 1 = n!$, exactly the number of permutations.

To generate a permutation given the odometer $[d_1, d_2, \dots, d_{n-1}]$, we can use the following method.

Start with $n$ empty circles in a row.

Write $1$ in the $(d_1+1)^{th}$ circle from the left.

Write $2$ in the $(d_2+1)^{th}$ empty circle from the left.

Write $3$ in the $(d_3+1)^{th}$ empty circle from the left.

and so on.

Once you have written $n-1$, there will be one empty circle left, where you write $n$, resulting in a permutation of $1,2,\dots, n$.

Notice that in the resulting permutation, the number of numbers greater than $i$, which appear to the left of $i$ is precisely $d_i$. In other words $[d_1, d_2, \dots, d_{n-1}]$ is the inversion table of the permutation, which tells you how many inversions a given number has.

It is not difficult to see that each odometer state generates a unique permutation.

If the odometer can be implemented as an iterator with has_next, next, so can the permutation generator.

There is an easy $O(n^2)$ time algorithm to map the odometer state to the permutation. Even $O(n\log n)$ is achievable.

Generating Cartesian Products


The Cartesian product of sets $A_1, A_2, \dots, A_n$ is the set of tuples $(a_1, a_2, \dots, a_n)$ such that $a_i \in A_i$.

If the size of $A_i$ is $s_i$, then the number of such tuples is $s_1\times s_2\times \dots \times s_n$.

An odometer for this is $[d_1, d_2, \dots, d_n]$ where $d_i$ cycles through $s_i$ digits. Given an odometer state $[d_1, d_2, \dots, d_n]$, we pick the $(d_i+1)^{th}$ element from $A_i$ for the corresponding element of the resulting tuple.

The mapping can be done in $O(n)$ time, assuming the $k^{th}$ element of $A_i$ can be picked in $O(1)$ time.

Generating Subsets



If we want to generate all the subsets of $\{a_1, a_2, \dots, a_n\}$, we can represent each set as a bit-vector $[v_1, v_2, \dots, v_n]$, where $v_i = 1$ iff $a_i$ is in the set.

We can use the odometer $[d_1, d_2, \dots, d_n]$ with each $d_i$ cycling through exactly two digits.

We pick $v_i = d_i$.

This mapping can be done in $O(n)$ time.

Conclusion


The generalized odometer seems to be a useful tool for use in implementing has_next, next type of iterators for certain combinatorial objects.

Also, given a number $M$ we can directly construct (and vice-versa) the odometer corresponding to $M$ increments, without having to increment $M$ times.

This direct mapping allows us to compute the $M^{th}$ permutation, Cartesian product, or subset directly (without having to perform $M$ increments). This would be useful to parallelize the generation and processing of those objects.


Monday, September 8, 2014

Iterators part one: Lazy shuffling with the Knuth Shuffle

[This is part one in a series about iterators.]

tl;dr: Generate a random permutation of $1,2 \dots, n$ lazily. Skip to the section titled Lazy Shuffle if you already know how the Knuth Shuffle works.

Introduction

Suppose you are working on something where you needed to create a random permutation of $1,2,\dots, n$ and process the random permutation one element at a time, potentially stopping the processing long before you run out of the n elements. How many elements you need to process is dependent on the elements you see (i.e. you don't know it before hand).

For example, n could be a billion, and you process the first 1000 elements in some cases, 100 in some etc.

A well known way to create a random permutation/shuffle is the Knuth Shuffle. In the Knuth shuffle, you are given an array, and the algorithm creates a random, in-place permutation of the array, in $O(n)$ time.

The idea of the Knuth shuffle is actually quite simple.

An explanation of Knuth Shuffle


The array consists of two portions, the shuffled portion and the yet to be shuffled portion. Initially the whole array is yet to be shuffled.

The shuffled portion will grow from the right to the left (i.e. initially empty, then $(A[n-1])$, then $(A[n-1], A[n-2])$ etc). You now pick a random object from the yet to be shuffled portion and move it to the location which will grow the shuffled portion by one. This will result in the displacement of a yet to be shuffled object. That object you move to the empty space created by the newly selected object that was just moved to the shuffled portion.


Lazy Shuffle


We could apply the Knuth shuffle to our problem, by creating an array $A$ of size n, such that $A[i] = i$, performing the Knuth shuffle on it, and process the elements in the order $A[0], A[1], \dots$.

Since n could be a billion, the time and space costs of this could be wasteful.

A practical approach would be to repeatedly generate a random number from $1$ to $n$ and check if you have generated that number before. If you have seen it before, you regenerate. This is a big improvement, but still wasteful, with the potential to run for a long time.

We can actually make the Knuth Shuffle work for us in a lazy fashion, without having to pay the $O(n)$ space or time costs.

Assuming we need to process only $m$ elements ($m$ unknown to us), there is a way to make only $m$ random number calls, and use $O(m)$ space, and run in $O(m)$ expected time.

We can achieve this by simulating the array in the Knuth shuffle by a hashtable.

[We can also use a sparse array implementation of the array, which will use $\Theta(n)$ space, but the running time will still be $O(m)$]

If we are simulating the array $A$, then the key will be $i$, and the value will be $A[i]$. The hashtable will only contain entries for elements that have been touched (read or written to). The first time we try to read $A[i]$, the hashtable will create an entry for it, setting the value to be $i$.

This allows us to provide the hasNext, next style iterator for the shuffled elements.

Here is what the code might look like (C++).

#include <map>  // For unordered_map
#include <random>  // Assume exists with some API.

// Sample usage:
// LazyShuffle ls(1000000);
// while (ls.has_next()) {
//     int v = ls.next();
//     Status s = process(v);
//     if (s.Done()) break;
// }
class LazyShuffle {
public:
    // generates a random permutation of 0 to n-1.
    LazyShuffle(int n)
        : _to_shuffle(n) {
    };

    bool has_next() {
        return _to_shuffle > 0;
    }

    int next() {
        // Do the Knuth shuffle here.
        // [yet to be shuffled] [shuffled]

        // Generate a random number from
        // 0 to _to_shuffle-1;
        _to_shuffle--;
        int selected = _random.generate(0, _to_shuffle);

        // Now swap this with element which 
        // grows the shuffled part
        swap(selected, _to_shuffle);

        return get(_to_shuffle);
    }

private:
    void set(int idx, int value) {
        _array[idx] = value;
    }

    int get(int idx) {
        if (_array.find(idx) == _array.end()) {
            _array[idx] = idx;
        }
        return _array[idx];
    } 

    void swap(int idx1, int idx2) {
        int val1 = get(idx1);
        int val2 = get(idx2);
        set(idx1, val2);
        set(idx2, val1);
    }

    // We use unordered_map, but we can use a hashtable
    // with expected O(1) time operations.
    std::unordered_map<int, int> _array; 

    // The number of elements not shuffled.
    int _to_shuffle;

    // Assume exists
    random::Generator _random;
};


Wednesday, September 3, 2014

A surprising property of power strings

[This came about trying to come up with a linear time algorithm to determine if a string is a power string]

Definition of power string.

A string $x$ is called a power string, if and only if there is some string $z$, such that $$x = z^m\quad \text{and}\quad m \gt 1$$

In other words, $x$ is $z$ concatenated to itself $m$ times.

For example, $nananana$ is the power string $(na)^4$, while  $nanananananabatman$ is not a power string.

[If you want to try finding a linear time algorithm for the detection of power strings, please stop reading]

The surprising property of power string is given by the following:

Proposition: Power strings are the only strings which are a cyclic shift of themselves.

We ignore the trivial shifts, by the length of the string which rotates every string into itself.

That power strings are cyclic shifts of themselves is easy to see. We show two proofs of the converse here.

Proof 1)

[Even though this proof seems verbose, it is not really difficult]

Suppose \(x\) is a string of length \(n\), over the alphabet \(A\), which is a cyclic shift of itself. You can basically treat \(x\) as a function \(f:\{0,1,2,\dots, n-1\}\to A\).

If \(x\) is a cyclic shift of itself, then for some \(1 \le k \le n-1\), \(f\) satisfies $$f(x+k \mod n) = f(x), \forall x \in \{0,1,2,\dots, n-1\}$$
In order to show that \(x\) is a power string, we need to show that there is a \(p\) such that \(p\) divides \(n\) and \(f(x+p \mod n) = f(x), \forall x \in \{0,1,2,\dots, n-1\}\)
 
Now if \(n = qk - k_1\) where \(q \gt 1\) and \(0 \le k_1 \lt k \), the we must have that

$$ f(x) = f(x + k \mod n) = f(x+ qk \mod n)$$
$$ = f(x + n + k_1 \mod n) = f(x+ k_1 \mod n)$$
Thus we get a sequence: \(k = k_0 \gt k_1 \gt k_2 \gt \dots \) such that $$f(x+k_i) = f(x), \forall x \in \{0,1,2,\dots, n-1\}$$
 Since \(k_i\) is a strictly decreasing sequence, we must end up at \(0\) at some point, say \(k_s = 0\). Then we can choose \(p = k_{s-1}\).

Proof 2)

We use Lyndon's theorem which says that two strings \(x,y\) commute (\(xy = yx\)) iff they are power strings of the same string, i.e. \(x = z^r\) and \(y = z^s\).

If \(x\) is a cyclic shift of itself, then we have that \(x = uv\) (concatenation of strings \(u\) and \(v\)) such that \( uv = vu\).

Linear Time Algorithm

This property allows you to determine if a string is a power string in \(O(n)\) time, using any linear time string matching algorithm like Boyer-Moore. All you need to do is check if \(x\) is a sub-string of \(xx\) with the first and last letters removed from \(xx\).

Saturday, August 30, 2014

The \(n^{th}\) root of \(n\)

For the sake of a quick, "pure" math post.

It is a well known fact that \(\lim_{n \to \infty} \sqrt[n]{n} = 1\)

There are multiple proofs of that. One such proof is to use the inequality of Arithmetic and Geometric means

$$ \frac{x_1 + x_2 + \dots + x_n}{n} \ge \sqrt[n]{x_1\cdot x_2\cdots x_n} \quad (x_i \ge 0)$$

Choose \(x_1 = x_2 = \dots = x_{n-2} = 1\) and \(x_{n-1} = x_n = \sqrt{n}\)

Then we get,

$$\frac{n-2 + 2\sqrt{n}}{n} \ge \sqrt[n]{n} \ge 1 $$

And so,

$$ 1 - \frac{2}{n} + \frac{2}{\sqrt{n}} \ge \sqrt[n]{n} \ge 1$$

By the sandwich theorem, the limit is one.

Friday, August 29, 2014

Hamiltonian Paths in Graph Powers

[Based on a puzzle given to me by Anand Ganesh]

Determining whether a connected graph \(G\) has a Hamiltonian Path is NP-Complete.

A little known fact is that determining if the square of a connected graph \(G\) has a Hamiltonian Path is also NP-Complete. In the square of \(G\), two vertices are adjacent iff they are at a distance at most two in \(G\).

This post is about the cube (distance at most three) of a connected graph.

It can be shown that the cube of any connected graph always contains a Hamiltonian Path, and there is a simple proof of that, which demonstrates the power of using a stronger induction hypothesis.

It is enough to prove that the cube of a tree is Hamiltonian, as given a graph, we can just consider a spanning tree.

Given a tree, assume it is a rooted tree (pick an arbitrary vertex as the root), so we can talk about children and parents etc. Now we try to arrange the vertices of the tree as \(v_1, v_2, \dots, v_n\) such that the distance between \(v_k\) and \(v_{k+1}\) is at most three (i.e. \(v_1, v_2, \dots, v_n\) is a Hamiltonian Path in the cube of the graph).

We now use the following induction hypothesis:

Proposition: The \(n \ge 2\) vertices of a rooted tree \(T\) can be arranged as \(v_1, v_2, \dots, v_{n}\) such that \(v_k\) and \(v_{k+1}\) are at a distance at most three in \(T\), and (this where we make the hypothesis stronger) \(v_1\) is the root, and \(v_n\) is a child of the root.

Now if the root has \(s\) sub-trees, we get the arrangements for the sub-trees (using the induction hypothesis), concatenate them, put the root at the end and reverse the resulting arrangement.

The resulting arrangement has the root as the first vertex, and a child of the root as the last. It is also easy to verify the distance at most three requirement: the distance between any child and any grand-child (possibly a child of a different child) is at most three.

This inductive proof can also be converted to a linear time algorithm for finding a Hamiltonian path in the cube of a graph.

Note that while trying to prove the property about Hamiltonian paths, we basically ended up proving the property about Hamiltonian cycles.

Related reading: http://en.wikipedia.org/wiki/Graph_power