## Friday, November 18, 2016

### Transistor radios

I was listening to Van Morrison's classic hit "Brown Eyed Girl" and the line about transistor radios always appeals to me as I am fascinated by how song lyrics capture science and technology of the era. Another example is the song "Kodachrome" by Paul Simon and I have written about fractals in the song "Frozen" in an earlier post. When I was younger, I built a crystal radio receiver and was amazed that I can listen to AM radio stations with a device that has no external power source and is powered solely by the energy in the radio waves! I also had a vintage transistor radio that is shaped like a pack of Marlboro cigarettes and powered by a 9V battery. The potentiometers for volume and channel have such a satisfying feel. The transistor radio was born as transistors have supplanted vacuum tubes in the demodulation and amplification circuitry of a radio receiver, resulting in portable receivers that are much smaller and use much less power. Today, portable radios are so much smaller and is packed with integrated circuits as more signal processing is done in the digital domain. The goal of software defined radio (SDR) is to replace all analog signal processing with digital signal processing done by computer algorithms by moving the analog to digital (A/D) conversion as early as possible. However, current computer processors and A/D converters do not have the speed and bandwidth to process digital samples at RF frequencies, so some analog processing is still done in SDR to demodulate the signals to intermediate frequencies or baseband frequencies. So the transistor is still necessary in radios today. But then again, even if the ideal SDR, where all processing is done by software, is possible, since computer processors today are still filled with transistors, technically the term "transistor radios" will be here to stay for a long time!

## Saturday, November 5, 2016

### "Bots" creating news digests about "bots"

I use a news digest app on my phone to read a selection of important news in areas I selected. This is nice since it only lists news items that I am most likely interested in and saves me time in reading the day's news. I think the news summary is created automatically via computer algorithms, since many times the highlighted quote does not make sense or is attributed to the wrong person. In addition, the "to explore further" section sometimes points to unrelated and/or inappropriate Wikipedia articles. Sometimes, this is because the person in the news has the same name as someone who is much more famous, so the Wikipedia article is about the wrong (but more well known) person.

Today, in the "science" section, there is a summary digest article about news organizations utilizing "bots", or automatic algorithms to aggregate data and generate news article. It is quite amusing, since again the highlighted quote is attributed to the wrong person and the "to explore further" section points to a specific news paper and is tangentially related to the news article.

Today, in the "science" section, there is a summary digest article about news organizations utilizing "bots", or automatic algorithms to aggregate data and generate news article. It is quite amusing, since again the highlighted quote is attributed to the wrong person and the "to explore further" section points to a specific news paper and is tangentially related to the news article.

## Wednesday, September 28, 2016

### Add oil

In Chinese slang, the phrase "加油", literally translated as "add oil", means "put more effort" and is typically used as encouragement to try harder in order to succeed at something. I believe the origin come from the fact that we need to add oil (gasoline) to cars in order for it to move. As we are all expected to drive electric or hydrogen cars in the future, this might become an archaic slang in the not too distant future.

After reading this post, my wife asked me what the corresponding Chinese slang should be for electric vehicles. She said that "充電" which is the translation of "charging electricity" is not a good choice since "充電" is slang for "refresh" or "renew" and is typically used when one is tired or drained of energy.

**Addendum: October 4, 2016**After reading this post, my wife asked me what the corresponding Chinese slang should be for electric vehicles. She said that "充電" which is the translation of "charging electricity" is not a good choice since "充電" is slang for "refresh" or "renew" and is typically used when one is tired or drained of energy.

## Tuesday, September 27, 2016

### Continued fraction expansion of the square root of n: part II

In an earlier blog post, $\delta(n)$ is defined as the smallest term in the periodic part of the continued fraction of $\sqrt{n}$ and I showed that if $r$ is even then $\delta((\frac{rm}{2})^2+m) = r$ for all $m\geq 1$. and if $r$ is odd, then $\delta((rm)^2+2m) = r$ for all $m\geq 1$. Note that $\delta(n)$ is only defined if $n$ is not a perfect square.

If you look at the first few numbers $n$ that satisfy $\delta(n) = r$, it would appear that they all follow the quadratric equations above. However, not all integers $n$ such that $\delta(n) = r$ are of the forms above. In particular, if $r > 0$ is even, then $\sqrt{\frac{r^4}{4} + r^3 + 2r^2 + 3r + 2} = \sqrt{\frac{(r^2-2)^2}{4}+(r+1)^3}$ has continued fraction expansion $\left[\frac{(r+1)^2+1}{2};\overline{r+1,r,r+1,(r+1)^2+1}\right]$ and thus $\delta\left(\frac{r^4}{4} + r^3 + 2r^2 + 3r + 2\right) = r$ and it is not of the forms above.

Similarly, if $r$ is odd, then $\sqrt{r^4 + r^3 + \frac{5(r+1)^2}{4}}$ has continued fraction expansion $\left[\frac{(r+1)(2r-1)+2}{2};\overline{r,2r-1,r,(r+1)(2r-1)+2}\right]$ and thus $\delta\left(r^4 + r^3 + \frac{5(r+1)^2}{4}\right) = r$ and it is not of the forms above either.

If you look at the first few numbers $n$ that satisfy $\delta(n) = r$, it would appear that they all follow the quadratric equations above. However, not all integers $n$ such that $\delta(n) = r$ are of the forms above. In particular, if $r > 0$ is even, then $\sqrt{\frac{r^4}{4} + r^3 + 2r^2 + 3r + 2} = \sqrt{\frac{(r^2-2)^2}{4}+(r+1)^3}$ has continued fraction expansion $\left[\frac{(r+1)^2+1}{2};\overline{r+1,r,r+1,(r+1)^2+1}\right]$ and thus $\delta\left(\frac{r^4}{4} + r^3 + 2r^2 + 3r + 2\right) = r$ and it is not of the forms above.

Similarly, if $r$ is odd, then $\sqrt{r^4 + r^3 + \frac{5(r+1)^2}{4}}$ has continued fraction expansion $\left[\frac{(r+1)(2r-1)+2}{2};\overline{r,2r-1,r,(r+1)(2r-1)+2}\right]$ and thus $\delta\left(r^4 + r^3 + \frac{5(r+1)^2}{4}\right) = r$ and it is not of the forms above either.

## Wednesday, September 21, 2016

### Paul Erdős and Kevin Bacon

In his famous 1967 paper, Stanley Milgram describes a study he conducted that shows that people are related to each other via a very small of number of acquaintances. This led to the phrase "six degrees of separation" being coined by John Guare in his play of the same name. Mathematicians study a similar concept called Erdős numbers. Two persons are linked if they have co-authored a mathematical paper together and a person's Erdős number is the minimum number of links between him/her and Paul Erdős. Thus Paul Erdős has Erdős number 0, A person other than Erdős who has written a paper with Erdős has Erdős number 1. A person who does not have Erdős number $\leq 1$ and has written a paper with a person with Erdős number 1 will have Erdős number 2, etc. If you have not written a paper with anyone with a (finite) Erdős number, then your Erdős number is $\infty$.

As a consequence of a drinking game, there is a similar notion among actors, called the Bacon number. Kevin Bacon has Bacon number 0 and the authors who are his co-stars in a movie have Bacon number 1, etc.

The fascinating aspect is that the Erdős number and the Bacon number of most people whose number is finite is relatively small, which is typically referred to as the "small world effect". There are various websites that lets you type in a name and it will attempt to find the Erdős or Bacon number of this person.

There is an additional notion of an Erdős-Bacon number which is the sum of a person's Erdős number and Bacon number. As of now the lowest Erdős-Bacon number appears to be 4.

There is something unsatisfactory about the definition of the Erdős-Bacon number. In particular, both Paul Erdős and Kevin Bacon have very large (and possibly infinite) Erdős-Bacon number. As of today, according to this link and this link, Paul Erdős has Bacon number $\infty$ and Kevin Bacon has Erdős number $\infty$.

There is one way to remedy this injustice. Kevin Bacon should publish a math paper with someone with Erdős number 1 (unfortunately Paul Erdős died in 1996) and starred in a movie with people who appeared in the documentary "N Is a Number: A Portrait of Paul Erdős". Since several mathematicians in the above documentary have Erdős number 1, both these activities can be combined if Kevin Bacon makes a documentary of how he collaborated with one (or more) of these mathematicians on a math paper. This will ensure that both Paul Erdős and Kevin Bacon have Erdős-Bacon number 2 (the lowest possible) and the universe will be in order again.

So Mr. Bacon, if you are reading this, please make that your next project!

As a consequence of a drinking game, there is a similar notion among actors, called the Bacon number. Kevin Bacon has Bacon number 0 and the authors who are his co-stars in a movie have Bacon number 1, etc.

The fascinating aspect is that the Erdős number and the Bacon number of most people whose number is finite is relatively small, which is typically referred to as the "small world effect". There are various websites that lets you type in a name and it will attempt to find the Erdős or Bacon number of this person.

There is an additional notion of an Erdős-Bacon number which is the sum of a person's Erdős number and Bacon number. As of now the lowest Erdős-Bacon number appears to be 4.

There is something unsatisfactory about the definition of the Erdős-Bacon number. In particular, both Paul Erdős and Kevin Bacon have very large (and possibly infinite) Erdős-Bacon number. As of today, according to this link and this link, Paul Erdős has Bacon number $\infty$ and Kevin Bacon has Erdős number $\infty$.

There is one way to remedy this injustice. Kevin Bacon should publish a math paper with someone with Erdős number 1 (unfortunately Paul Erdős died in 1996) and starred in a movie with people who appeared in the documentary "N Is a Number: A Portrait of Paul Erdős". Since several mathematicians in the above documentary have Erdős number 1, both these activities can be combined if Kevin Bacon makes a documentary of how he collaborated with one (or more) of these mathematicians on a math paper. This will ensure that both Paul Erdős and Kevin Bacon have Erdős-Bacon number 2 (the lowest possible) and the universe will be in order again.

So Mr. Bacon, if you are reading this, please make that your next project!

## Thursday, August 25, 2016

### The 2016 Rio Olympics

It has been a thrilling Summer Olympics in the last 2 weeks and we enjoyed watching many of the events on TV. While watching the swimming competition, we noticed that Nathan Adrian looks surprisingly similar to Chow Yun Fat (周潤發), especially when he smiles. Chow Yun Fat is one of my favorite Hong Kong actors and I grew up watching him in several TVB TV series, with 北斗雙雄 being my favorite series (having watched it several times over the years). Perhaps they can cast Mr. Adrian when they need to make a biopic of Mr. Chow.

Another thing during the Olympics that I like is that I can listen to a special CD. Many years ago, (around the 1996 Atlanta Olympics I believe), I took some pictures on film (yes, there used to be such a thing as photographic film) and went to the local drugstore to get them developed. There was a special promotion from Kodak that included a CD titled "The Sound and the Spirit" with music from various Olympic games. Since then I would play this CD every 4 (sometimes 2) years.

Another thing during the Olympics that I like is that I can listen to a special CD. Many years ago, (around the 1996 Atlanta Olympics I believe), I took some pictures on film (yes, there used to be such a thing as photographic film) and went to the local drugstore to get them developed. There was a special promotion from Kodak that included a CD titled "The Sound and the Spirit" with music from various Olympic games. Since then I would play this CD every 4 (sometimes 2) years.

### The Hartman-Grobman Linearization Theorem

*Theorem*: In the neighborhood of a hyperbolic fixed point, a smooth vector field or a diffeomorphism is topologically conjugate to its linear part.

This result was proved by Grobman and Hartman independently around 1959-1960 and basically states that the dynamics near a hyperbolic fixed point is essentially the same as the dynamics of its linearization which we can characterize completely from the eigenvalues pattern. This is true for both continuous-time dynamics (vector field) or discrete-time dynamics (diffeomorphism).

Here is a sketch of the standard proof for the case of a diffeomorphism. First, we need the following simple fact for linear maps in Banach spaces: if $F$ is an invertible contraction, then $I+F^{-1}$ is also invertible. This can be seen as follows. $I+F^{-1} = F^{-1}(I+F)$. If $I+F$ is not invertible, then there exists $x\neq y$ such that $x+F(x) = y+F(y)$. This implies that $x-y = F(y)-F(x)$, i.e. $\|x-y\| = \|F(y)-F(x)\|$, contradicting the fact that $F$ is a contraction. Therefore $I+F$ is invertible, and thus $I+F^{-1}$ is invertible since it is the product of two invertible maps.

Consider a diffeomorphism $f$ with a hyperbolic fixed point at $0$. Let $A$ be the linear part of $f$ at $0$. We want to find a homeomorphism $h = I+\delta$ such that $fh = hA$. As we are interested only at $f$ near a neighborhood of $0$, we can assume that $f$ can be written as $f = A+\phi_1$ such that $\phi_1$ is bounded and have a small Lipschitz constant. Furthermore, $\phi_1$ can be chosen small enough such that $A+\phi_1$ is a homeomorphism. Consider the equation $(A+\phi_1)h

= h(A+\phi_2)$. After using the fact that $h=I+\delta$ and some manipulation, we get the following Eq. (1):

\[\delta - A^{-1}\delta(A+\phi_2) = A^{-1}(\phi_2-\phi_1(I+\delta))

\]

Next we argue that the linear operator $H: \delta \rightarrow \delta - A^{-1} \delta(A+\phi_2)$ is invertible.

By hyperbolicity of $A$, we can decompose the phase space into the stable subspace $W^s$ and the unstable subspace $W^u$. Since $W^s$ and $W^u$ are invariant under $A^{-1}$, if $\delta$ is a bounded function into $W^s$ and $W^u$ then $H(\delta)$ is also a bounded function into $W^s$ and $W^u$ respectively. Split $\delta = \delta^s + \delta^u$ into two functions $\delta^s$ and $\delta^u$ which maps into $W^s$ and $W^u$ respectively.

The map $\delta^s \rightarrow A^{-1}\delta^s(A+\phi_1)$ is invertible with inverse $\delta^s \rightarrow A \delta^s (A+\phi_1)^{-1}$ since $A \delta^s (A+\phi_1)^{-1} = A^s\delta^s(A+\phi_1)^{-1}$ the map

$\delta^s \rightarrow A \delta^s (A+\phi_1)^{-1}$ is a contraction and therefore

the map $\delta^s \rightarrow \delta^s - A^{-1}\delta^s(A+\phi_1)$ is invertible based on the fact discussed before. The same thing can be done with the $W^u$ and this implies that $H$ is invertible.

Coming back to Eq. (1) above, we get

\[ \delta = H^{-1}A^{-1}(\phi_2-\phi_1(I+\delta)) = \psi(\delta)\]

For small $\phi_1$ and $\phi_2$, $\psi$ is a contraction and thus for given $\phi_1$ and $\phi_2$ there exists a unique $\delta$ and hence a unique $h$. It can be shown that $h$ is a homeomorphism and by choosing $\phi_2 = 0$ we get the desired result.

**References**

D. M. Grobman, "Homeomorphisms of systems of differential equations," Doklady Akademii Nauk SSSR, vol. 128, pp. 880–881, 1959.

P. Hartman, "A lemma in the theory of structural stability of differential equations," Proc. AMS, vol. 11, no. 4, pp. 610–620, 1960.

## Wednesday, August 17, 2016

### Rounding the k-th root of n

Consider the problem of finding the $k$-th root of a number $n\geq 0$ and rounding it to the nearest integer, i.e. find $[\sqrt[k]{n}]$, where $[x]$ is $x$ rounded to the nearest integer. This can be easily computed in many computer languages using floating point arithmetic, but care must be taken for large $n$ to ensure enough significant digits are available. On the other hand, languages such as Python has built-in support for integers of arbitrary sizes and will automatically allocate more space to fit the number under consideration. This can be used to compute $[\sqrt[k]{n}]$ using only integer arithmetic without worrying whether there are enough precision in the floating point representation.

Let $i$ be the largest integer such that $i \leq \sqrt[k]{n}$. The number $i$ can be computed using integer arithmetic with an iterative Newton's method.

Since $n \geq 0$, $[\sqrt[k]{n}] = i+1$ if $\sqrt[k]{n}-i \geq \frac{1}{2}$ and $[\sqrt[k]{n}] = i$ otherwise. The condition $\sqrt[k]{n}-i \geq \frac{1}{2}$ is equivalent to $2^k n \geq (2i+1)^k$ which can be computed using integer arithmetic.

A simple python function using the gmpy2 module to implement this is the following:

from gmpy2 import iroot

def round_root(n,k): # round(k-th root of n), n >= 0

i = iroot(n,k)[0]

return int(i) + int(2**k*n >= (2*i+1)**k)

The gmpy2 module also includes the functions isqrt_rem and iroot_rem. The function isqrt_rem(n)returns a pair of numbers $i,j$ such that $i$ is the largest integer $\leq \sqrt{n}$ and $j = n-i^2$.

Similarly, iroot_rem(n,k)returns a pair of numbers $i,j$ such that $i$ is the largest integer $\leq \sqrt[k]{n}$ and $j = n-i^k$.

Since

\begin{eqnarray*}(2i+1)^k &=& (2i)^k + (2i+1)^{k-1} + \\

&&(2i+1)^{k-2}2i + \cdots + (2i+1)(2i)^{k-2} + (2i)^{k-1}\end{eqnarray*}

the condition can be rewritten as:

\begin{eqnarray*}2^k j &\geq &(2i+1)^{k-1} + (2i+1)^{k-2}2i + \cdots + (2i+1)(2i)^{k-2} + (2i)^{k-1}\\ & \geq & \sum_{m=0}^{k-1} (2i+1)^{k-1-m}(2i)^m \end{eqnarray*}

For $k=2$, this is reduced to: $4j \geq 4i + 1$. A python function implementing $[\sqrt{n}]$ is:

from gmpy2 import isqrt_rem

def round_sqrt(n): # round(square root of n), n >= 0

i, j = isqrt_rem(n)

return int(i) + int(4*(j-i) >= 1)

Similarly, for $k=3$, the condition is reduced to $8j \geq 6i(2i+1)+1$.

Let $i$ be the largest integer such that $i \leq \sqrt[k]{n}$. The number $i$ can be computed using integer arithmetic with an iterative Newton's method.

Since $n \geq 0$, $[\sqrt[k]{n}] = i+1$ if $\sqrt[k]{n}-i \geq \frac{1}{2}$ and $[\sqrt[k]{n}] = i$ otherwise. The condition $\sqrt[k]{n}-i \geq \frac{1}{2}$ is equivalent to $2^k n \geq (2i+1)^k$ which can be computed using integer arithmetic.

A simple python function using the gmpy2 module to implement this is the following:

from gmpy2 import iroot

def round_root(n,k): # round(k-th root of n), n >= 0

i = iroot(n,k)[0]

return int(i) + int(2**k*n >= (2*i+1)**k)

The gmpy2 module also includes the functions isqrt_rem and iroot_rem. The function isqrt_rem(n)returns a pair of numbers $i,j$ such that $i$ is the largest integer $\leq \sqrt{n}$ and $j = n-i^2$.

Similarly, iroot_rem(n,k)returns a pair of numbers $i,j$ such that $i$ is the largest integer $\leq \sqrt[k]{n}$ and $j = n-i^k$.

Since

\begin{eqnarray*}(2i+1)^k &=& (2i)^k + (2i+1)^{k-1} + \\

&&(2i+1)^{k-2}2i + \cdots + (2i+1)(2i)^{k-2} + (2i)^{k-1}\end{eqnarray*}

the condition can be rewritten as:

\begin{eqnarray*}2^k j &\geq &(2i+1)^{k-1} + (2i+1)^{k-2}2i + \cdots + (2i+1)(2i)^{k-2} + (2i)^{k-1}\\ & \geq & \sum_{m=0}^{k-1} (2i+1)^{k-1-m}(2i)^m \end{eqnarray*}

For $k=2$, this is reduced to: $4j \geq 4i + 1$. A python function implementing $[\sqrt{n}]$ is:

from gmpy2 import isqrt_rem

def round_sqrt(n): # round(square root of n), n >= 0

i, j = isqrt_rem(n)

return int(i) + int(4*(j-i) >= 1)

Similarly, for $k=3$, the condition is reduced to $8j \geq 6i(2i+1)+1$.

## Wednesday, July 27, 2016

### A way for Spiderman to catch the Green Goblin

Some years ago, we worked on a problem in digital halftoning for which we could prove a result [1] which states that there is a variant of the vector error diffusion algorithm such that the sum of the errors is bounded. This implies that the average error decreases to 0 as the number of errors to be averaged increases. This is paraphrased as saying that averaging a large area of the halftoned image is similar to the average of the original image over the same area, which is what we expect a halftoning algorithm to behave. In the course of this research, I came up with the following scenario to which this result also provide a solution to, and I have used it to describe this problem to a lay audience. In celebration of Spiderman's well received introduction to the Marvel Cinematic Universe, I thought I'll describe it here:

Consider a city whose shape is a convex polygon, with a building at each corner of this polygon. The main villain Green Goblin (GC) is loose in the city and it is up to Spiderman (S), our friendly neighborhood superhero, to catch him. At the start, GC and S are located at different places in the city. Because of fatique or lack of fuel, at each time epoch, both GC and S are moving less and less. However, whereas GC can move arbitrarily within the city, S (being a webslinger) can only move toward a building at the corner of the city, along the line connecting S and this building.

More precisely, time is divided into epochs numbered 1, 2, 3, ...

At the k-th epoch:

[1] R. Adler, B. Kitchens, M. Martens, A. Nogueira, C. Tresser, C. W. Wu, "Error bounds for error diffusion and related digital halftoning algorithms," Proceedings of IEEE International Symposium on Circuits and Systems, 2001, pp. II-513-516.

Consider a city whose shape is a convex polygon, with a building at each corner of this polygon. The main villain Green Goblin (GC) is loose in the city and it is up to Spiderman (S), our friendly neighborhood superhero, to catch him. At the start, GC and S are located at different places in the city. Because of fatique or lack of fuel, at each time epoch, both GC and S are moving less and less. However, whereas GC can move arbitrarily within the city, S (being a webslinger) can only move toward a building at the corner of the city, along the line connecting S and this building.

More precisely, time is divided into epochs numbered 1, 2, 3, ...

At the k-th epoch:

- GC picks a destination within the city to move to and moves 1/k of the distance in the direction towards the destination.
- S can moves 1/k of the distance to a building located on the corners of the city, in the direction of that building.

**References**[1] R. Adler, B. Kitchens, M. Martens, A. Nogueira, C. Tresser, C. W. Wu, "Error bounds for error diffusion and related digital halftoning algorithms," Proceedings of IEEE International Symposium on Circuits and Systems, 2001, pp. II-513-516.

## Saturday, July 23, 2016

### Lava, donuts and printing

We have just returned from vacation touring the British isles. In particular, we visited the Giant's causeway in Belfast, Northern Ireland. The area is covered by more than 40,000 basalt columns with many of them having a shape close to a hexagonal column (similar to the shape of a pencil). They were formed when lava were cooling and shrinking. Why the hexagonal shapes? This might be due to the fact that a hexagonal lattice arrangement of circles on a plane is the densest arrangement of circles of the same size [1]. The hexagonal lattice arrangement is also the covering of circular disk with the least overlap. Many arrangements of objects in nature such as honeycombs have this hexagonal arrangement. The Voronoi regions of this arrangment are hexagons. Since this arrangement is periodic in two (non perpenticular) axes, one can view this as a periodic packing on the torus, provided the density matches the dimension of the torus. When the density does not match, we don't have a hexagonal packing on the torus, and it is not clear what the densest packing is. This has applications in digital halftoning [2]. The only difference in the halftoning application is that the circle centers are points on a discrete grid on the torus. We studied 2 algorithms to generate such packing on the torus. The first algorithm is the Direct Binary Search (DBS) algorithm [3] and generate patterns like this:

The second algorithm is based on the Riesz energy minimization theory of Hardin and Saff [4] and we were able to obtain patterns that are more uniform than DBS:

In both cases, they look like the patterns found on Giant's causeway:

In 3 dimensions, the densest packing is the face-centered cubic (FCC) lattice packing (also known as the "canonball" packing). It was conjectured to be the densest packing by Kepler in 1611, and it was only proved relatively recently by Hales with a proof first announced in 1998 and the correctness of the lengthy proof checked by computer in 2014.

[1] J. Conway and N. J. A. Sloane, "Sphere Packings, Lattices and Groups", Springer, 3rd Edition, 1998.

[2] C. W. Wu, B. Trager, K. Chandu and M. Stanich, "A Riesz energy based approach to generating dispersed dot patterns for halftoning applications," Proceedings of SPIE-IS&T Electronic Imaging, SPIE, vol. 9015, pp. 90150Q, 2014,

[3] J. P. Allebach, “DBS: retrospective and future directions,” in Proceedings of SPIE, 4300, pp. 358–376, 2001.

[4] D. P. Hardin and E. B. Saff, “Discretizing manifolds via minimum energy points,” Notices of the AMS 51(10), pp. 1186–1194, 2004.

The second algorithm is based on the Riesz energy minimization theory of Hardin and Saff [4] and we were able to obtain patterns that are more uniform than DBS:

In both cases, they look like the patterns found on Giant's causeway:

In 3 dimensions, the densest packing is the face-centered cubic (FCC) lattice packing (also known as the "canonball" packing). It was conjectured to be the densest packing by Kepler in 1611, and it was only proved relatively recently by Hales with a proof first announced in 1998 and the correctness of the lengthy proof checked by computer in 2014.

**References:**[1] J. Conway and N. J. A. Sloane, "Sphere Packings, Lattices and Groups", Springer, 3rd Edition, 1998.

[2] C. W. Wu, B. Trager, K. Chandu and M. Stanich, "A Riesz energy based approach to generating dispersed dot patterns for halftoning applications," Proceedings of SPIE-IS&T Electronic Imaging, SPIE, vol. 9015, pp. 90150Q, 2014,

[3] J. P. Allebach, “DBS: retrospective and future directions,” in Proceedings of SPIE, 4300, pp. 358–376, 2001.

[4] D. P. Hardin and E. B. Saff, “Discretizing manifolds via minimum energy points,” Notices of the AMS 51(10), pp. 1186–1194, 2004.

## Wednesday, July 20, 2016

### Musical tastes across the generations

One day my son and I were eating at a restaurant and we heard an instrumental piece playing over the speakers. He said to me: "doesn't this sound familiar?". I listened to it for a while and indeed it sounded familiar. It brought back memories from my childhood, as the music is from a song from Hong Kong that I remember when I was much younger. But as my son grew up in the USA, how come it is familiar to him as well? Some quick internet searching reveals that the song is a 1963 instrumental piece called "Washington Square" by the Village Stompers, a band from New York City. It turns out the music was used in 1978 by Sam Hui (許冠傑), one of my favorite Hong Kong singers from my childhood, in a song called 學生哥, roughly translated as "student". It is a catchy song encouraging students to study and become productive members of society. One of the reason Sam Hui's songs are popular is because they use Cantonese slang a lot, as contrasted to more formal Cantonese used by others. Sometimes the slang words are purely verbal, and there are no written Chinese logograms to describe these words, so he would put English words to phonetically describe these words and many of his song lyrics are full of such words.

Anyway, the reason my son heard of this song that was composed before I was born is because it has been covered in a more modern version called "Faidherbe square" by ProleteR and was used in a Youtube video my son was watching. Of course he prefers the modern version while I prefer the Hong Kong version, but we both liked the original 1963 composition. It is remarkable this music has delighted multiple generations and I am glad we both discovered this music together.

Anyway, the reason my son heard of this song that was composed before I was born is because it has been covered in a more modern version called "Faidherbe square" by ProleteR and was used in a Youtube video my son was watching. Of course he prefers the modern version while I prefer the Hong Kong version, but we both liked the original 1963 composition. It is remarkable this music has delighted multiple generations and I am glad we both discovered this music together.

### My teenage son wants to go grocery shopping with us

The other day I was going to the grocery store to do some shopping and my son said to me: "I want to go with you to the grocery store." This has never happened before. He is also spending more time outside. I believe it is all due to this life changing program that my son has downloaded on his phone. It is a proven system (proven within the last week or so) that has gotten many people exercising more and spending more time outside. The basic premise is simple; the goal is to catch small monsters (so called pocket monsters) while wandering around the neighborhood. All I can say is Bravo to the creators of this app!

## Monday, June 20, 2016

### Meta recurrences

In my freshman year at college, my roommate Joe introduced me to the book "Gödel, Escher, Bach: an Eternal Golden Braid" by Douglas Hofstadter that I really enjoyed reading. In this book, Hofstadter considers sequences whose recurrence relation is such that the arguments themselves depend on the terms of the sequence. For instance, the classical Fibonacci sequence is defined by

Several variants of meta Fibonacci sequences have been analyzed in Refs. [1-4].

For some simpler meta recurrences, much more can be said about them. For instance Golomb showed that sequences satisfying the recursion

In particular, he showed that the recursion

Next we show that this last recursion can be generalized as follows. For an integer

The cases of

For

[1] S.M Tanny "A well-behaved cousin of the Hofstadter sequence," Discrete Math., 105 (1992), pp. 227–239.

[2] Callaghan, Joseph, John J. Chew III, and Stephen M. Tanny. "On the behavior of a family of meta-Fibonacci sequences." SIAM Journal on Discrete Mathematics 18.4 (2005): 794-824.

[3] B. Balamohan, A. Kuznetsov and S. Tanny, On the behavior of a variant of Hofstadter's Q-sequence, J. Integer Sequences, Vol. 10 (2007), #07.7.1.

[4] C. Deugau and F. Ruskey, Complete k-ary Trees and Generalized Meta-Fibonacci Sequences, J. Integer Seq., Vol. 12.

*f(n) = f(n-1) + f(n-2)*, with*f(1) = f(2) = 1*. Hofstadter considered a meta Fibonacci sequence defined by*q(1) = q(2) = 1*, and*q(n) = q(n-q(n-1)) + q(n-q(n-2))*. Not much is known about this sequence, except that it behaves wildly. It is not even known whether the sequence is well defined!Several variants of meta Fibonacci sequences have been analyzed in Refs. [1-4].

For some simpler meta recurrences, much more can be said about them. For instance Golomb showed that sequences satisfying the recursion

*a(n) = a(n - a(n-1))*is eventually periodic when it is well-defined. He also showed several meta recurrences with explicit solutions, such as*a(n) = a(a(n))*.In particular, he showed that the recursion

*a(n) = a(n-a(n-1)) + 1*with*a(1) = 1*has a simple solution described by the sequence 1, 2, 2, 3, 3, 3, 4, 4, 4, 4, ....Next we show that this last recursion can be generalized as follows. For an integer

*k*, define the meta recurrence*a(n) = a(n-a(n-k))+k*with initial conditions*a(i) = i*for*i = 1,...,k.*This recurrence has a simple solution and it can be constructed starting with*km*copies of*k(m+1)*, followed by*km+1*,*km+2*, ...,*k(m+1)*, as*m*varies from 0, 1, 2, ... It is straightforward to verify that this construction satisfies the recurrence relation. Thus for each*k*the sequence corresponding to the meta recurrence is well defined, every positive integer is in the sequence, and every integer not a proper multiple of*k*appears once. If*t*is a multiple of*k*, then*t*appears*t-k+1*times.The cases of

*k=2*and*k=3*are shown in OEIS sequences A193358 and A274213 respectively.For

*k = 4*, the sequence is: 1, 2, 3, 4, 8, 8, 8, 8, 5, 6, 7, 8, 12, 12, 12, 12, 12, 12, 12, 12, 9, 10, 11, 12, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 16, 13, 14, 15, 16, ....**References**[1] S.M Tanny "A well-behaved cousin of the Hofstadter sequence," Discrete Math., 105 (1992), pp. 227–239.

[2] Callaghan, Joseph, John J. Chew III, and Stephen M. Tanny. "On the behavior of a family of meta-Fibonacci sequences." SIAM Journal on Discrete Mathematics 18.4 (2005): 794-824.

[3] B. Balamohan, A. Kuznetsov and S. Tanny, On the behavior of a variant of Hofstadter's Q-sequence, J. Integer Sequences, Vol. 10 (2007), #07.7.1.

[4] C. Deugau and F. Ruskey, Complete k-ary Trees and Generalized Meta-Fibonacci Sequences, J. Integer Seq., Vol. 12.

### Continued fraction expansion of the square root of n

Results from Galois imply that the continued fraction expansion of $\sqrt{n}$ is of the form $[a_0; \overline{a_1,a_2,\cdots, a_2, a_1, 2a_0}]$ (see also Ref.[1], page 469).

In OEIS sequences A031710, A031712, A031713, A031749 and several other sequences, the least number among the periodic part of the continued fraction expansion is considered, i.e what is $\delta(n) \triangleq \min (a_1,a_2,\cdots, a_2, a_1, 2a_0)$?

If $t$ is an integer that divides $2k$, then it is straightforward to verify that the continued fraction expansion of $\sqrt{(km)^2+tm}$ for $m \geq 1$ is equal to $[km; \overline{\frac{2k}{t}, 2km}]$. This implies the following result:

If $r$ is even, then $\delta((\frac{rm}{2})^2+m) = r$ for all $m\geq 1$.

If $r$ is odd, then $\delta((rm)^2+2m) = r$ for all $m\geq 1$.

[1] Charles Smith, A Treatise on Algebra, Fifth Edition, Macmillan and Co., 1896.

In OEIS sequences A031710, A031712, A031713, A031749 and several other sequences, the least number among the periodic part of the continued fraction expansion is considered, i.e what is $\delta(n) \triangleq \min (a_1,a_2,\cdots, a_2, a_1, 2a_0)$?

If $t$ is an integer that divides $2k$, then it is straightforward to verify that the continued fraction expansion of $\sqrt{(km)^2+tm}$ for $m \geq 1$ is equal to $[km; \overline{\frac{2k}{t}, 2km}]$. This implies the following result:

If $r$ is even, then $\delta((\frac{rm}{2})^2+m) = r$ for all $m\geq 1$.

If $r$ is odd, then $\delta((rm)^2+2m) = r$ for all $m\geq 1$.

**References:**[1] Charles Smith, A Treatise on Algebra, Fifth Edition, Macmillan and Co., 1896.

## Saturday, May 28, 2016

### Generating functions of sequences

The generating function of a sequence

$\{a(n)\}_{n=0}^\infty = \{a_0, a_1, a_2, \cdots\}$ for integers $n\geq 0$ is defined as $f(x) = a_0 + a_1x + a_1x^2 + \cdots$.

For a sequence that satisfies the linear recurrence relation

$a(n) = c_0a(n-1) + c_1a(n-2) + \cdots + c_{m}a(n-m-1)$

with initial terms $a(0) = a_0, a(1) = a_1, \cdots , a(m) = a_m$,

the general form of the generating function is given by:

$$ f(x) = \frac{a_mx^m + \sum_{i=0}^{i< m}\left(a_ix^i- c_ix\sum_{j=0}^{j < m-i}a_jx^{i+j}\right)}{

1-x\sum_{i=0}^{i \leq m}c_ix^{i}}

$$

The sequence $a(n) = bc^n+d$ satisfies the recurrence relation $a(n) = (c+1)a(n-1)-ca(n-2)$ with

generating function:

$$ (b(1-x) + d(1- cx ))/((1-x )(1-cx)) $$.

$\{a(n)\}_{n=0}^\infty = \{a_0, a_1, a_2, \cdots\}$ for integers $n\geq 0$ is defined as $f(x) = a_0 + a_1x + a_1x^2 + \cdots$.

For a sequence that satisfies the linear recurrence relation

$a(n) = c_0a(n-1) + c_1a(n-2) + \cdots + c_{m}a(n-m-1)$

with initial terms $a(0) = a_0, a(1) = a_1, \cdots , a(m) = a_m$,

the general form of the generating function is given by:

$$ f(x) = \frac{a_mx^m + \sum_{i=0}^{i< m}\left(a_ix^i- c_ix\sum_{j=0}^{j < m-i}a_jx^{i+j}\right)}{

1-x\sum_{i=0}^{i \leq m}c_ix^{i}}

$$

The sequence $a(n) = bc^n+d$ satisfies the recurrence relation $a(n) = (c+1)a(n-1)-ca(n-2)$ with

generating function:

$$ (b(1-x) + d(1- cx ))/((1-x )(1-cx)) $$.

### How to lose weight without working out

Today, my son said: "Since we are spinning along with the earth, shouldn't the centrifugal force make us feel lighter?" The answer is yes, since the centrifugal force cancels out a little bit of the gravitational pull of the Earth. In fact, that is exactly what a satellite does. It spins so fast around the earth that the centrifugal force cancels out the gravitational pull so it is "floating" in space. The centrifugal force is a virtual or fictitious force denoting the opposite of the force (the centripetal force) needed to achieve the circular motion and this centripetal force is supplied by the force of gravity. There is another interpretation where the satellite is forever falling towards earth, but the motion also has a velocity component that is perpendicular to the force of gravity and its direction changes (due to the gravitational acceleration), but the magnitude doesn't.

How much lighter do you weight on the equator, where you are spinning the fastest, versus on the poles where you are not spinning at all. There are 2 reasons why you weight is different on the poles versus the equator. The first reason is the centrifugal force described above.

This force is given by Newton's second law: $F = ma$ along with $a = \omega^2r$. $\omega$ is the angular velocity of the earth rotation, i.e. $\frac{2\pi}{T}$, where $T$ is about 24 hours = $86400s$. The radius at the equator is about $r$ = 6378137 m. This results in a=0.03373 which is about 0.34% of the gravitational acceleration of $g =9.8 \frac{m}{s^2}$

The second reason is that the earth is not a perfect sphere but is squashed on the poles since the centrifugal force of the spinning earth pulls the matter of the earth from the center of mass and this happens more at the equator than at the poles. Thus you are further away from the center of mass of the Earth on the surface of the equator than on the surface of the pole. This affects the value of $g = \frac{MG}{r^2}$ where $M$ is the mass of the Earth and $G$ is the universal gravitational constant. This adds another 0.67% to the difference resulting in you (and any mass) weighting about 1% less on the equator than at the poles.

How much lighter do you weight on the equator, where you are spinning the fastest, versus on the poles where you are not spinning at all. There are 2 reasons why you weight is different on the poles versus the equator. The first reason is the centrifugal force described above.

This force is given by Newton's second law: $F = ma$ along with $a = \omega^2r$. $\omega$ is the angular velocity of the earth rotation, i.e. $\frac{2\pi}{T}$, where $T$ is about 24 hours = $86400s$. The radius at the equator is about $r$ = 6378137 m. This results in a=0.03373 which is about 0.34% of the gravitational acceleration of $g =9.8 \frac{m}{s^2}$

The second reason is that the earth is not a perfect sphere but is squashed on the poles since the centrifugal force of the spinning earth pulls the matter of the earth from the center of mass and this happens more at the equator than at the poles. Thus you are further away from the center of mass of the Earth on the surface of the equator than on the surface of the pole. This affects the value of $g = \frac{MG}{r^2}$ where $M$ is the mass of the Earth and $G$ is the universal gravitational constant. This adds another 0.67% to the difference resulting in you (and any mass) weighting about 1% less on the equator than at the poles.

## Sunday, May 15, 2016

### Interesting number: 55

The number 55 = 5 x 11 = (10 x 11)/2 is a semiprime, a Fibonacci number, a triangular number, a square pyramidal number, a squarefree number, a palindrome in base 10, and a Stirling number of the second kind.

## Tuesday, May 3, 2016

### Ironic song while driving this morning

While driving to work this morning, the sky was gray and overcast and there is a light drizzle. Sometime during the commute, the song "I can see clearly now" by Johnny Nash came on the music streaming service. I got a chuckle out of it as the lyrics are in direct contradiction to what I was seeing around me.

## Thursday, April 14, 2016

### Interesting number: 348588115

The number 348588115 = 5 × 23 × 3031201 has some interesting properties. It is a sphenic number (i.e. it is the product of 3 distinct prime numbers). Reversing the decimal digits also results in a sphenic number: 511885843 = 7 × 953 × 76733. In addition, the decimal digits also are a concatenation of the decimal digits of 2 other sphenic numbers:

'348588115' = '34858' + '8115' and 34858 = 2 × 29 × 601, 8115 = 3 × 5 × 541.

'348588115' = '34858' + '8115' and 34858 = 2 × 29 × 601, 8115 = 3 × 5 × 541.

## Tuesday, March 15, 2016

### Grocery shopping and Braess' paradox

Some time ago, our local grocery store introduced self-checkout, where we can scan our grocery items ourselves. My wife and I like to use this option since there are several stations open and not many people are using it, so we can get through the checkout quicker. However, recently we noticed that the lines at the self-checkout are getting longer, while the lines at the checkout clerks are sometimes shorter than the self-checkout lines. I believe this is an example of Braess' paradox. Braess noticed that when a new road is added to a network of existing roads, the congestion can actually increase. His reasoning is that each driver follows a selfish routing algorithm and chooses the shortest route to the destination. If the new road is a shortcut, then everyone who can use this road will choose it, and this can actually increase the congestion through the network. This seems to be happening with the self-checkout lines as well. Everyone assumes it is quicker to use them and this leads to longer lines. One way to achieve optimal throughput in a road network is if there is a centralized control to route each driver, but this would require information about the speed and travel plans of all drivers and it might favor some drivers over others for the overall good. To implement this would require cars to talk to each other using technologies such as connected vehicles, vehicle-to-vehicle communications and dedicated short range communication (DSRC).

Back to our grocery shopping trips. Yesterday the self-checkout lines was quite long so we went to the checkout clerk who is much faster than us in scanning the groceries (plus there was an additional clerk to bag our groceries), and we got done much quicker than if we had used the self-checkout line. Sometimes one or more of the self-checkout stations are out of order which is similar to lane closures on a multi-lane highway. If more people realize that regular checkout is faster, we could have another Braess' paradox in the making and lead to longer lines at the regular checkout lines. In the long run, this can lead to cycles of short and long lines forming alternately at the self-checkout lanes and at the regular checkout stations! Fortunately this is unlikely to happen since the 2 types of checkout lines are in close proximity so the shopper has a good sense of the congestion of the entire system. But imagine a mega-store where for some odd reason the regular checkout lines are far away from the self-checkout lines! The shopper needs to decide which line is more optimal based on incomplete information and what other shoppers would do. On the other hand, the speed of checkout may not be the only objective. Some people prefer self checkout because they feel like they don't make mistakes and make better bagging decisions (like not having eggs at the bottom of the bag) if they do it themselves.

Back to our grocery shopping trips. Yesterday the self-checkout lines was quite long so we went to the checkout clerk who is much faster than us in scanning the groceries (plus there was an additional clerk to bag our groceries), and we got done much quicker than if we had used the self-checkout line. Sometimes one or more of the self-checkout stations are out of order which is similar to lane closures on a multi-lane highway. If more people realize that regular checkout is faster, we could have another Braess' paradox in the making and lead to longer lines at the regular checkout lines. In the long run, this can lead to cycles of short and long lines forming alternately at the self-checkout lanes and at the regular checkout stations! Fortunately this is unlikely to happen since the 2 types of checkout lines are in close proximity so the shopper has a good sense of the congestion of the entire system. But imagine a mega-store where for some odd reason the regular checkout lines are far away from the self-checkout lines! The shopper needs to decide which line is more optimal based on incomplete information and what other shoppers would do. On the other hand, the speed of checkout may not be the only objective. Some people prefer self checkout because they feel like they don't make mistakes and make better bagging decisions (like not having eggs at the bottom of the bag) if they do it themselves.

## Monday, March 14, 2016

### Happy Pi Day 2016!

Today, March 14, is Pi Day. As I wrote in an earlier post this year's Pi Day is specially nice since 3.14.16 is Pi rounded to 4 places after the decimal point.

## Saturday, February 27, 2016

### Van der Pol Oscillator

The original neon-bulb oscillator circuit studied by Van der Pol and Van der Mark in their 1927 paper consists of a neon-bulb connected to a capacitor, resistor, a bias voltage source and a sinuosoidal signal voltage source. The operation of this circuit without the sinusoidal signal can be described as follows: the neon bulb is not conducting and the voltage source charges up the the capacitor until a certain point at which the neon bulb starts conducting (and glowing) and discharges the capacitor until the voltage is low enough to turn off the neon bulb and the cycle starts anew. Looking at the $v$-$i$ characteristic of the neon bulb we see that it is a function with two local extrema just like the cubic polynomial. The role of the bias voltage source is to add power to the system, since the neon bulb by itself does not supply any power. In particular, by looking at the equivalent $v$-$i$ characteristic of the neon bulb in parallel with the bias voltage, it effectively shift the characteristic so it looks more like the cubic and thus can oscillate.

Except for the external driving source, this circuit is similar to the circuits considered in a previous post. To write state equations, we model the neon-bulb as a nonlinear resistor in parallel with an small linear inductor and a small linear capacitor.

When the sinusoidal source is turned off, the oscillator oscillates freely with a frequency $\omega_1$. The addition of the sinusoidal voltage source add another frequency $\omega_2$ to the system. One can think of the Poincaré map (generated by sampling the flow at frequency $\omega_2$) as a map on the circle, the circle being the limit cycle. In fact, for certain parameters, we see behavior similar to those obtained from the circle map.

The complete state equations of the circuit is given by:

\[ \begin{array}{lcl} \frac{dv}{dt} & = &

-\frac{1}{RC}v - \frac{1}{C}i + \frac{E}{RC} \\

\frac{di}{dt} & = & \frac{1}{L}v - \frac{f(i)}{L}

- \frac{E_0 \sin(\omega t)}{L}

\end{array}

\]

By adjusting the capacitor $C$, we change the free running frequency of the oscillator, i.e. changing the ratio between the two frequency. We observe frequency locking where the oscillator oscillates at a submultiple of the driving frequency for a range of capacitor values. If we vary the driving frequency, we observe the Devil's staircase, similar to the one observed in the circle map.

Van der Pol, B. and Van der Mark, J., “Frequency demultiplication”, Nature, 120, 363-364, 1927.

Kennedy, M. P. & Chua, L. O., “Van der Pol and chaos,” IEEE Trans. Circuits Syst. CAS-33, 974–980, 1986.

Kennedy, M. P., Krieg, K. R. & Chua, L. O., “The Devil’s staircase: The electrical engineer’s fractal,” IEEE Trans. Circuits Syst. CAS-36, 1133–1139, 1989.

Except for the external driving source, this circuit is similar to the circuits considered in a previous post. To write state equations, we model the neon-bulb as a nonlinear resistor in parallel with an small linear inductor and a small linear capacitor.

When the sinusoidal source is turned off, the oscillator oscillates freely with a frequency $\omega_1$. The addition of the sinusoidal voltage source add another frequency $\omega_2$ to the system. One can think of the Poincaré map (generated by sampling the flow at frequency $\omega_2$) as a map on the circle, the circle being the limit cycle. In fact, for certain parameters, we see behavior similar to those obtained from the circle map.

The complete state equations of the circuit is given by:

\[ \begin{array}{lcl} \frac{dv}{dt} & = &

-\frac{1}{RC}v - \frac{1}{C}i + \frac{E}{RC} \\

\frac{di}{dt} & = & \frac{1}{L}v - \frac{f(i)}{L}

- \frac{E_0 \sin(\omega t)}{L}

\end{array}

\]

By adjusting the capacitor $C$, we change the free running frequency of the oscillator, i.e. changing the ratio between the two frequency. We observe frequency locking where the oscillator oscillates at a submultiple of the driving frequency for a range of capacitor values. If we vary the driving frequency, we observe the Devil's staircase, similar to the one observed in the circle map.

**References:**Van der Pol, B. and Van der Mark, J., “Frequency demultiplication”, Nature, 120, 363-364, 1927.

Kennedy, M. P. & Chua, L. O., “Van der Pol and chaos,” IEEE Trans. Circuits Syst. CAS-33, 974–980, 1986.

Kennedy, M. P., Krieg, K. R. & Chua, L. O., “The Devil’s staircase: The electrical engineer’s fractal,” IEEE Trans. Circuits Syst. CAS-36, 1133–1139, 1989.

### Triangular numbers

A triangular number $T(n)$ is defined as $n(n+1)/2$. It is equal to the binomial coefficient $\left(\begin{array}{c}n+1\\2\end{array}\right)$ and is also equal to $0+1+2+\cdots+n$. When is an integer $m\geq 0$ equal to a triangular number $T(n)$ for some integer $n\geq 0$? We have $n(n+1) = 2m$. Using the quadratic formula to solve the quadratic equation $n^2+n-2m = 0$ for $n$, we obtain one solution for $n$:

$$ n = \frac{\sqrt{8m+1}-1}{2} $$

The other solution for $n$ is negative (or complex) so we will not use that.

If $8m+1$ is not the square of an integer, then $n$ is not an integer. If $8m+1$ is the square of an integer, then $8m+1$ is odd and thus $\sqrt{8m+1}$ is odd as well. This means that $n = \frac{\sqrt{8m+1}-1}{2}$ is an integer. Thus $m$ is a triangular number if and only if $8m+1$ is the square of an integer.

$$ n = \frac{\sqrt{8m+1}-1}{2} $$

The other solution for $n$ is negative (or complex) so we will not use that.

If $8m+1$ is not the square of an integer, then $n$ is not an integer. If $8m+1$ is the square of an integer, then $8m+1$ is odd and thus $\sqrt{8m+1}$ is odd as well. This means that $n = \frac{\sqrt{8m+1}-1}{2}$ is an integer. Thus $m$ is a triangular number if and only if $8m+1$ is the square of an integer.

## Sunday, February 14, 2016

### What not to listen to while driving

We were driving on the road today and listening to "In the shadow of the mountain king" from Peer Gynt by Edvard Grieg when my son remarked: "This is not a good piece of music to listen to when driving, it makes you drive faster and faster." My wife replied: "The same is true for Bolero". I tend to agree. Luckily, my driving is not affected by the tempo of the music and we reached our destination safely.

## Friday, February 12, 2016

### The curse and blessing of exponential growth

Exponential growth is a common phenomena occurring all around us. Sometimes it is good, other times not. Simply speaking, we talk about exponential growth if some quantity is increasing as an exponential function of another quantity, i.e. y = ab

^{x}.**The Good:**- Typically when you put your money
*n*years. But because of compound interest, i.e. your interest will also earn interest, you will receive 5% of \$105 = \$5.25 in interest the second year and the money in the bank will be will be \$100*1.05^{n}after*n*years, i.e. an exponential growth. If you can appreciate how fast exponential growth is, this is the best reason for anyone to save money early. - In cryptography, making the number of possibilities that an encryption key can take on large is an important idea to make cryptographic techniques secure. This is akin to choosing a long password. The number of combinations grows exponentially with the number of digits. For instance, a n-digit decimal number has $10^n$ possible combinations.

**The not so Good:**- In the analysis of computer algorithms, many known algorithms for solving hard NP-complete problems (such as 3SAT) have time complexity exponential in the problem size. This means that the size of the problem that can be solved will be limited.
- Many models of growth assumes an exponential growth rate, but this cannot be sustained indefinitely in a finite universe. Some examples are Moore's Law or models of population growth. Models of chaotic systems such as the Lorenz system and Chua's circuit also experience exponential divergence of trajectories on a small scale, but on a larger scale the trajectories remain bounded.

## Monday, February 8, 2016

### Parity of binomial coefficients

Consider the binomial coefficient C(n,m) = n*(n-1)*...*(n-m+1)/(m*(m-1)*...*2*1). When is this number odd and when is it even? In this post, we use Lucas' theorem to provides a simple characterization of this. Lucas' theorem states that for integers $m$, $n$ and prime $p$, the following relationship holds:

\[ C(n,m) \equiv \prod_{i=0}^{k} C(n_i, m_i) \mod p \]

where $n_i$ and $m_i$ are the digits of $n$ and $m$ in base $p$ respectively, with $n_0$ and $m_0$ being the least significant digits.

When $p=2$, $n_i$ and $m_i$ are the bits in the binary expansion of $n$ and $m$ and $C(n_i, m_i)$ is $0$ if and only if $n_i < m_i$. This implies that $C(n,m) \equiv 0 \mod 2$ if and only if $n_i < m_i$ for some $i$.

The truth table for $n_i < m_i$ is:

\[\begin{array}{ c | c || c }

n_i & m_i & n_i < m_i \\ \hline

0 & 0 & 0 \\

0 & 1 & 1 \\

1 & 0 & 0 \\

1 & 1 & 0 \\

\end{array}\]

i.e., it is logically equivalent to $m_i AND ( NOT n_i)$. If we think of $AND$ and $NOT$ as bitwise operations on the binary representation of numbers, then we have shown that:

$C(n,m) \equiv 0 \mod 2$ if and only if $m AND ( NOT n)$ is not zero.

A simple Python function that implements this statement is

def bincoefmod2(n,m):

# computes binomial(n,m) mod 2

return 0 if ~n & m else 1

or equivalently:

def bincoefmod2(n,m):

# computes binomial(n,m) mod 2

return int(not ~n & m)

Incidentally, for bits $n_i$ and $m_i$, $n_i < m_i$ is logically equivalent to $NOT (m_i\Rightarrow n_i)$.

Next we consider C(n,m)C(a,b) mod 2. Clearly this is equivalent (mod 2) to (C(n,m) mod 2)*(C(a,b) mod 2). Thus $C(n,m)C(a,b) \equiv 0 mod 2$ if and only if $m AND ( NOT n)$ is not zero or $b AND ( NOT a)$ is not zero. This in turns implies the following:

\[\begin{array}{l}

C(n,m)C(a,b) \equiv 0 mod 2 \Leftrightarrow \\

(m AND ( NOT n)) OR (b AND ( NOT a)) \ne 0

\end{array}\]

where OR is a bitwise operation. The corresponding Python function is:

def bincoefprodmod2(n,m,a,b):

# computes binomial(n,m)*binomial(a,b) mod 2

return 0 if (~n & m)|(~a & b) else 1

or equivalently:

def bincoefprodmod2(n,m,a,b):

# computes binomial(n,m)*binomial(a,b) mod 2

return int(not (~n & m)|(~a & b))

\[ C(n,m) \equiv \prod_{i=0}^{k} C(n_i, m_i) \mod p \]

where $n_i$ and $m_i$ are the digits of $n$ and $m$ in base $p$ respectively, with $n_0$ and $m_0$ being the least significant digits.

When $p=2$, $n_i$ and $m_i$ are the bits in the binary expansion of $n$ and $m$ and $C(n_i, m_i)$ is $0$ if and only if $n_i < m_i$. This implies that $C(n,m) \equiv 0 \mod 2$ if and only if $n_i < m_i$ for some $i$.

The truth table for $n_i < m_i$ is:

\[\begin{array}{ c | c || c }

n_i & m_i & n_i < m_i \\ \hline

0 & 0 & 0 \\

0 & 1 & 1 \\

1 & 0 & 0 \\

1 & 1 & 0 \\

\end{array}\]

i.e., it is logically equivalent to $m_i AND ( NOT n_i)$. If we think of $AND$ and $NOT$ as bitwise operations on the binary representation of numbers, then we have shown that:

$C(n,m) \equiv 0 \mod 2$ if and only if $m AND ( NOT n)$ is not zero.

A simple Python function that implements this statement is

def bincoefmod2(n,m):

# computes binomial(n,m) mod 2

return 0 if ~n & m else 1

or equivalently:

def bincoefmod2(n,m):

# computes binomial(n,m) mod 2

return int(not ~n & m)

Incidentally, for bits $n_i$ and $m_i$, $n_i < m_i$ is logically equivalent to $NOT (m_i\Rightarrow n_i)$.

Next we consider C(n,m)C(a,b) mod 2. Clearly this is equivalent (mod 2) to (C(n,m) mod 2)*(C(a,b) mod 2). Thus $C(n,m)C(a,b) \equiv 0 mod 2$ if and only if $m AND ( NOT n)$ is not zero or $b AND ( NOT a)$ is not zero. This in turns implies the following:

\[\begin{array}{l}

C(n,m)C(a,b) \equiv 0 mod 2 \Leftrightarrow \\

(m AND ( NOT n)) OR (b AND ( NOT a)) \ne 0

\end{array}\]

where OR is a bitwise operation. The corresponding Python function is:

def bincoefprodmod2(n,m,a,b):

# computes binomial(n,m)*binomial(a,b) mod 2

return 0 if (~n & m)|(~a & b) else 1

or equivalently:

def bincoefprodmod2(n,m,a,b):

# computes binomial(n,m)*binomial(a,b) mod 2

return int(not (~n & m)|(~a & b))

### Chinese New Year 2016

Today, February 8, 2016, is Chinese New Year. Just for fun, I looked through OEIS to see what interesting facts I can find about the leap year 2016. Here are some facts about the number 2016:

- 2016 is the sum of the divisors of both the 42th and the 47th generalized pentagonal number.
- 2016 = 3^3+...+9^3 is the sum of 7 positive cubes.
- The sum of divisors of 2016 squared is a multiple of 2016.
- The 2016th Fibonacci number is a multiple of the sum of divisors of 2016 (= 6552)
- Both 2016 and 2016^2 has 0 as smallest digit and 6 as largest digit.
- The product of the proper divisors of 2016 is a multiple of the sum of the proper divisors of 2016.
- Both the number of divisors of 2016 and the number of positive integers relatively prime to 2016 that are less than or equal to 2016 are perfect squares.
- 2016 = 63*64/2 is a triangular number, i.e, it is the sum 1+2+...+63.
- 2016 = 6!+36^2.
- 2016 (in base 10) is divisible by the sum of their digits in every base from 2 through to 16.

Subscribe to:
Posts (Atom)