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$.

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. GC picks a destination within the city to move to and moves 1/k of the distance in the direction towards the destination.
  2. S can moves 1/k of the distance to a building located on the corners of the city, in the direction of that building.
The question is: can S ever catch up to GC?  The result in [1] shows that the answer is yes and gives an explicit algorithm of which building S should swing from at each epoch in order for S to catch up to GC.

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.

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.

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 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$.

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)) $$.

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.