Tautochrone curve on a Spherical surface

The Tautochrone problem is one of the classical problems in math. Both the Wikipedia and the Mathworld page present some great math on the problem and makes for an enjoyable read.

I was going through these pages one more time recently and I thought "All of this is fine. But is there a tautochrone curve on a sphere?". Turns out the answer is yes and this post is all about finding such a curve.

The simplest way to attack this problem is to use Able's solution idea. It says,

$\displaystyle \frac{ds}{dy}=\pm\sqrt{2g}\frac{T_0}{\pi}\frac{1}{\sqrt{y}}$

where $T_0$ is the constant time of travel, $g$ is the gravitational constant, $s$ is the arc length of the curve and $y$ is the measure of the gravitational potential energy.

Let's consider a unit sphere centered at the origin. Also, Let the plane tangent to the south pole of the sphere be the 'base' for measuring the gravitational potential.

As we are working on a spherical surface, we'll naturally use spherical geometry for this problem with $\theta$ as the polar angle and $\phi$ as the azimuthal angle. We have,

$y=1+\cos\theta$ which implies $dy=-\sin\theta d\theta$.

Therefore, Abel's equation becomes,

$\displaystyle ds=\sqrt{2g}\frac{T_0}{\pi}\frac{\sin\theta}{\sqrt{1+\cos\theta}}d\theta=\sqrt{2g}\frac{T_0}{\pi}\sqrt{1-\cos\theta}\text{ }d\theta$

The differential arc length in a spherical coordinate system is given by,

$\displaystyle \frac{ds}{d\theta}=\sqrt{1+\sin^2\theta\left(\frac{d\phi}{d\theta}\right)^2}$

Using this and simplifying, we get

$\displaystyle \sin^2\theta\left(\frac{d\phi}{d\theta}\right)^2=2g\left(\frac{T_0}{\pi}\right)^2(1-\cos\theta)-1$

Now, let $2g(T_0/\pi)^2=1/(1-\cos\lambda)$ where $\lambda$ is a constant. Then the above equation with further simplication becomes,

$\displaystyle \frac{d\phi}{d\theta}=\frac{1}{\sin\theta}\sqrt{\frac{\cos\lambda-\cos\theta}{1-\cos\lambda}}$

This is the differential equation that must be satisfied by the tautochrone curve.

We can solve this with the same substitution used in the Mathworld page. We then have,

$\displaystyle \phi=\cot(\lambda/2)\tanh^{-1}\left(\sqrt{\frac{\cos\lambda-\cos\theta}{1+\cos\lambda}}\right)-\tan^{-1}\left(\sqrt{\frac{\cos\lambda-\cos\theta}{1-\cos\lambda}}\right)$

The above expression should also have an arbitrary constant of integration, but let's say we only care about the shape of the curve and hence leave it out of this analysis.

Also, from the expression above, we see that the curve stats at $\theta=\lambda$ and spirals its way down to the south pole. We thus know, once we specify $T_0$, we indirectly specify where the curve should start from and vice versa.

Tautochrone curve with $\lambda=\pi/8$

We know the time of descent of the tautochrone curve is $T_0$. The arc length $L$ of this tautochrone curve is given by,

$\displaystyle L(\theta_0)=\int ds=\int \frac{ds}{d\theta}d\theta=\sqrt{2g}\frac{T_0}{\pi}\int\limits_{\theta_0}^\pi \sqrt{1-\cos\theta}d\theta=\sqrt{2g}\frac{T_0}{\pi}\cdot2\sqrt{2}\cos(\theta_0/2)$

or in terms of $\lambda$, we have,

$\displaystyle L(\theta_0)=2\frac{\cos\theta_0/2}{\sin\lambda/2}$

All of this absolutely amazing. But the story isn't complete yet. There is something more to this tautochrone curve that truly blows our mind which I'll share in the next post. Hope you enjoyed the post.

Mathematica code for the above curve:

Clear["Global`*"];

\[Alpha] = 45 Degree/2;

Phif[\[Theta]_] := 

  Cot[\[Alpha]/2] ArcTanh[Sqrt[(Cos[\[Alpha]] - Cos[\[Theta]])/(

     1 + Cos[\[Alpha]])]] - 

   ArcTan[Sqrt[(Cos[\[Alpha]] - Cos[\[Theta]])/(1 - Cos[\[Alpha]])]];

Cartf[\[Theta]_] := {Cos[Phif[\[Theta]]] Sin[\[Theta]], 

   Sin[\[Theta]] Sin[Phif[\[Theta]]], Cos[\[Theta]]};

Show[Graphics3D[{Opacity[1/5], White, Sphere[]}],

 ParametricPlot3D[Cartf[\[Theta]], {\[Theta], 0 Degree, 180 Degree}, 

  PlotStyle -> Thick, ColorFunction -> (Hue[0.8 #4] &)], 

 ParametricPlot3D[{Cos[t], Sin[t], 0}, {t, 0, 2 Pi}], 

 ParametricPlot3D[{Cos[t], 0, Sin[t]}, {t, 0, 2 Pi}], 

 ParametricPlot3D[{0, Cos[t], Sin[t]}, {t, 0, 2 Pi}], 

 ImageSize -> Large, Boxed -> False]

See here to paste this code into Mathematica with proper formatting.

Until then

Yours Aye

Me

Gergonnes's Solution to Apollonius' Problem

I was recently revisiting my earlier post on A nice result in Geometry and got 'distracted' into the Classical Apollonius problem. The fact that this problem has attracted the interests some of the most famous mathematicians throughout history tells us how beautiful this problem really is.

As always, the wiki page was a great source of information for me and it gives a very nice account on Gergonne's solution to the problem which, needless to say, is ingenious. Meanwhile, I was also consulting Cut-The-Knot's version of Gergonne's solution.

Both the solutions seemed similar but there was one subtle difference. To make it more clear (and readable), I (re)created the diagrams on Wikipedia using Geogebra. Before we dwell further, let's list out the key labels used in the construction.

In what follows , $i \in \{1,2,3\}$ and $X \in \{A,B\}$

$C_1,C_2,C_3$ - The given circles. In general, there are eight circles that are mutually tangent to all these circles
$C_O$ - Circle orthogonal to all the three given circles
$G$ - Radical center of the given circles. Also, center of the Orthogonal circle $C_O$
$C_A,C_B$ - A pair (out of the possible eight) of solution circles
$X_i$ - Points where the solution circle $C_X$ touches Circle $C_i$
$l_i$ - Line passing through $A_i$ and $B_i$. As $A_i$ and $B_i$ are inverses to each other w.r.t to $C_O$, this line also contains the Radical center $G$; also the polar of point $T_i$ defined below
$p$ - (One out of the possible four) Homothetic Axis. This line also happens to be the radical axis of the solution circles
$P_i$ - Pole of the Homothetic axis $p$ w.r.t to circle $C_i$
$J_i, K_i$ - Points of intersection of the Orthogonal circle $C_O$ and circle $C_i$
$T_i$ - Point of intersection of the line joining $J_i$ and $K_i$, and the Homothetic axis $p$; Also, the pole of line $l_i$ defined above
$O_1,O_2$ - Points of intersection between the Orthogonal circle $C_O$ and Homothetic axis $p$
$C_I$ - Circle with center at $O_1$ and passing through $O_2$

Gergonne's Solution - Wikipedia
(i) Construct $p$, the Homothetic axes (or axes of similitude) (ii) Construct the Radical center $G$ of the three given circles (iii) Construct the poles $P_i$ of the Homothetic axes w.r.t the given circles.

Points $P_i$ are the objects of interest in Wikipedia version of Gergonne's Solution

Now the line joining $G$ and $P_i$ cuts the circle $C_i$ in $A_i$ and $B_i$ giving six points in total. These six points determine the solution circles.

Gergonne's Solution - Cut-The-Knot
(i) Construct $p$, the Homothetic axes (or axes of similitude) (ii) Construct the circle orthogonal to the given circles $C_O$ (iii) Construct the Radical axis of the Orthogonal circle and each of the given circles (that is, the line joining $J_i$ and $K_i$); this intersects the Homothetic axes in $T_i$.

Points $T_i$ are the objects of interest in Cut-The-Knot version of Gergonne's Solution

Now construct tangent from $T_i$ to each of the given circles.These points then determine the solution circles.


The third step in these solutions were the source of my confusion. Especially, while Wikipeida's version needs an inversion (constructing the poles $P_i$), Cut-the-Knot's version needs no inversion at all. This baffled me and it took me some time (two days to be exact) to figure out how this was possible.

The explanation given in Cut-The-Knot did not seem to clarify this point. It quotes Four Touching Circles IV, which does not add anything in relation to the third step in the construction.

Gergonne's solution relies crucially on the inversive relation between poles and polars of conics. In particular, it exploits the fact that if the pole $P$ of a line $p$ w.r.t a conic lies on a line $q$, then the pole $Q$ of line $q$ w.r.t the same conic lies on line $p$. Wikipedia version (and most descriptions of Gergonne's solution) makes use of this relation.

Now, Cut-The-Knot's solution uses the mirror image of the same relation. That is, if the polar $p$ of a point $P$ contains a point $Q$, then the polar $q$ of point $Q$ w.r.t the same conic, contains the point $P$. And hence, it cleverly avoids any explicit inversion construction. Here's how it works.

From Wikipedia, we know that point $T_i$ is the pole of the (polar) line $l_i$ w.r.t circle $C_i$. Now the polar of point $T_i$ contains the point $G$. Hence, the polar of $G$ must contain $T_i$.

Here's comes the nice part. Given any pair of orthogonal circles, the polar of center of either circle w.r.t to the other circle is their radical axis. This is very obvious from the construction of polars and orthogonal circles.

So in essence, we are actually constructing a polar (rather than a pole in Wikipedia's version) of point. However, because of the clever choice, it just happens to be the point of intersection of circles already constructed avoiding the need for any inversion.

That was indeed neat. Hope you agree.

This whole process of revisiting one of the classical problems in Geometry was really nice experience for me. I greatly enjoyed creating the Geogebra drawing of the solution. In fact, I came up with something on my own.

Whenever, the Homothetic axis $p$ meets the Orthogonal circle $C_O$, something very interesting happens. As the solution circles are inverses of each other under $C_O$, the points of intersection between $p$ and $C_O$ also happens to be the points of intersection between $C_A$ and $C_B$.

Now draw a circle $C_I$ with center at one of the points of intersection (Point $O_1$) and passing through the other (Point $O_2$). Use this as a circle of inversion.

As the Orthogonal circle passes through the center of inversion, it becomes a straight line $C_O'$.
The three given circles, being orthogonal to the $C_O$, transforms into circles with their center on the line above.
The solution circles too passes through the inversion center and they too transform into straight lines ($C_A'$ and $C_B'$).
The (other) point of intersection $O_2$ of the solution circles passes through the inversion circle $C_I$ and hence this point is where the solution lines $C_A'$ and $C_B'$ meet. That is, point $O_2$ becomes the Homothetic center of the (inverted) given circles.

The given circles (not shown here) inverted under $C_I$ form a set of Homothetic circles with center $O_2$.
The Orthogonal circles and both the solution circles transforms to straight lines.

In summary, if $p$ meets $C_O$, we can construct a circle $C_I$ with center at $O_1$ and passing through $O_2$. Invert any of the given circle and draw tangent lines from $O_2$ to the inverted circle. Inverting back these tangent lines gives the required solution circles.

While constructing inversion is harder (needs more steps) while doing manually, this method is relatively simpler if we have an 'macro' for inversion (as in Computer graphics systems like Geogebra).

But the main drawback about this method is that it only works when the Homothetic axes intersects the Orthogonal circle which, sadly, is not guaranteed when the Homothetic axis contains two external Homothetic centers. Even in this case, we can actually develop a solution using the concept of mid-circles. This is not complicated but is more laborious. Gergonne's solution can be used in these cases.

That's all for today. Hope you enjoyed it.


Until then
Yours aye
Me

Error terms of the Exponential Function's limit form

One of the classical characterizations of the Exponential function is to define it by the following limit.

$e^x=\displaystyle\lim_{n \to \infty}\left(1+\frac{x}{n}\right)^n$

Characterizations of the exponential function gives the following as the first few error terms of this limi-expression,

$\displaystyle\left(1+\frac{x}{n}\right)^n=e^x\left(1-\frac{x^2}{2n}+\frac{x^3(8+3x)}{24n^2}-\cdots\right)$

stating that the numerators of the $n^k$ terms are polynomials in $x$ of degree $2k$. However, no mention is made about the computation of the coefficients or their relation with each other (if there is any).

That exactly is the point of this post. We try to derive an expression which allows us to calculate the coefficients as we desire. Considering that this definition of the Exponential function is all over the internet, I was disappointed that I could not find any reference regarding these coefficients.

Based on Wikipedia's expression, we suspect that,

$\displaystyle e^{-x}\left(1+\frac{x}{n}\right)^n=\sum_{i\geq j \geq 0}(-1)^{i}p_{i,j}\frac{x^{i+j}}{n^i}$

where we seek to compute the coefficients $p_{i,j}$. Now using the substitutions $x \to -nx$ and $n \to -1/n$ (in that order), we can see that,

$e^{-x/n}(1-x)^{-1/n}=\displaystyle\sum_{i\geq j \geq 0}p_{i,j}\frac{x^{i+j}}{n^j}$

But $e^{-x/n}(1-x)^{-1/n}=e^{-(x+\ln{(1-x)})/n}=\displaystyle\sum_{j=0}^\infty(-x+\ln{(1-x)^{-1}})^j\frac{1}{n^jj!}$

Comparing the coefficients of powers of $n$, we have,

$x^jp_j(x)=\displaystyle x^j\sum_{i=0}^\infty p_{i,j}x^i=\frac{1}{j!}(-x+\ln{(1-x)^{-1}})^j$

This readily gives a combinatorial representation. That is, $(i+j)!p_{i,j}$ is the number of derangements of $[1,2,3,\cdots,i+j]$ with exactly $j$ cycles. Therefore,

$\displaystyle (i+j)!p_{i,j}=\sum_{k=0}^j (-1)^{j-k}\binom{i+j}{i+k}\begin{bmatrix}i+k\\k\end{bmatrix}=\sum_{k=0}^j (-1)^{i+j-k}\binom{i+j}{i+k}s(i+k,k)$

where $\begin{bmatrix}n\\k\end{bmatrix}$(and $s(n,k)$) are the unsigned (and signed) Stirling numbers of the first kind.

Also, using the same derangement interpretation, we can easily see that,

$\displaystyle \sum_{j=0}^i i!p_{i-j,j}=i!\sum_{j=0}^i\frac{(-1)^j}{j!}\implies \sum_{j=0}^i p_{i-j,j}=\sum_{j=0}^i\frac{(-1)^j}{j!}$

On the other hand, we can derive a recurrence for $p_{i,j}$'s by manipulating $x^jp_j(x)$.

$jxp_j(x)=p_{j-1}(x)(-x+\ln(1-x)^{-1})=\displaystyle p_{j-1}(x)\left(\frac{x^2}{2}+\frac{x^3}{3}+\frac{x^4}{4}+\cdots\right)$

Comparing coefficients of powers of $x$, we have,

$p_{i,j}=\displaystyle\frac{1}{j}\sum_{k=2}^{i+1}\frac{1}{k}p_{i-k+1,j-1}=\frac{1}{j}\sum_{k=0}^{i-1}\frac{p_{k,j-1}}{i-k+1}$

Together with $p_{0,0}=1$ and $p_{i,j}=0$ when $i<j$ (or) $j=0$ enables us to calculate all the coefficients.

Further, we can differentiate $x^jp_j(x)$ to get a different recurrence.

$\displaystyle p_{i,j}=\frac{1}{i+j}p_{i-1,j-1}+\left(1-\frac{1}{i+j}\right)p_{i-1,j}$

This is really magical as we now have a probabilistic interpretation for the $p_{i,j}$'s.

Consider a bot that is initially at point $(i,j)$ with $i\geq j$ in the $x$-$y$ plane. The bot is programmed to move towards the origin such that at each step it either moves down-left to $(i-1,j-1)$ with some probability or moves left to point $(i-1,j)$ with the complementary probability.

If the probability of moving down-right is the inverse of the Manhattan distance between the bot's current position and the origin, then $p_{i,j}$ is the probability that bot reaches the origin without ever touching the $x$-axis or crossing the $x=y$ line.

This is simply wonderful as there does not seem to be a connection between the above probability and the error terms of the exponential function. Again the recurrence with the same boundary conditions as before gives us all the $p_{i,j}$'s.

We conclude the post with more terms of the error terms.

$\displaystyle\left(1+\frac{x}{n}\right)^n=e^x\left(1-\frac{x}{2}\frac{x}{n}+\left(\frac{x}{3}+\frac{x^2}{8}\right)\frac{x^2}{n^2}-\left(\frac{x}{4}+\frac{x^2}{6}+\frac{x^3}{48}\right)\frac{x^3}{n^3}+\left(\frac{x}{5}+\frac{13x^2}{72}+\frac{x^3}{24}+\frac{x^4}{384}\right)\frac{x^4}{n^4}-\cdots\right)$

Noting the linear term in $n$, I had the idea of using the generalized Binomial series to get the following alternate limit form. This is a better form of the limit expression as the resulting asymptotic series does not have linear term either in $n$ or in $x$.

$e^x=\displaystyle\lim_{n \to \infty}\left(1+\frac{x}{n}\right)^{n+x/2}$

That's all for today. Hope you enjoyed it.


Until then
Yours Aye
Me

A Binary Grid counting Problem - Part II

In this post, we consider the same problem as in the previous post but for grids of odd size. As much of the setup was discussed in the previous post, let's dive in right away.

Identity - Let $t_n$ denote total number of such colorings on an $n \times n$ grid. As discussed in the previous post, we have A001499 that gives us a recurrence. With $t_0=1$, we have

$t_n=n(n-1)t_{n-1}+\frac{1}{2}n(n-1)^2t_{n-2}$

Diagonal (Anti-Diagonal) symmetry - Let $d_n$ denote the required colorings that are symmetric along the main diagonal. Again A000986 gives us the following recurrence. Starting with $d_0=1$,

$d_n=2(n-1)d_{n-1}-(n-1)(n-3)d_{n-2}-\frac{1}{2}(n-1)(n-2)(n-3)d_{n-4}$

Horizontal (Vertical) reflection - This is really easy. Consider a $(2n+1) \times (2n+1)$ grid. The '$n+1$'th row will have a black square. As each column must have two black squares, there must be another black square in this column. If the second black square is on the top half of the grid, the horizontal reflection would result in an additional black square in the bottom half of the grid, making three squares in that column.Hence, there cant be any colorings that are fixed by horizontal/vertical reflection in an odd sized grid.

180 deg. rotation symmetry - Let $f_n$ denote the number of colorings of a $(2n+1) \times (2n+1)$ grid that are fixed by rotating the grid by 180 degrees. I found this to be the most interesting case. Proceeding like the previous post, we have,

$\begin{align}
f_n&=[x_1^2x_2^2\cdots x_n^2x_{n+1}t_1^2t_2^2\cdots t_n^2t_{n+1}]\displaystyle\prod_{i=1}^n(1+x_{n+1}t_i)(1+x_it_{n+1})\prod_{j=1}^n(1+x_it_j)^2\\
\end{align}$

Unfortunately, I couldn't simplify it any further. I computed the first few terms using Mathematica and searched for a generating function. And I found one.

$\displaystyle f(z)=\sum_{n=0}^\infty f_n \frac{z^n}{n!^2}=\frac{2ze^{-z}}{(1-4z)^{3/2}}$

Now let's derive a (holonomic) recurrence out of this. Let $f(z)=g(z)h(z)$ where $g(z)=ze^{-z}$ and $h(z)=(1-4z)^{-3/2}$. It is see that $(1-4z)h'=6h$ and $zg'=(1-z)g$. From this, we get,

$z(1-4z)f'(z)=6zf(z)+(1-z)(1-4z)f(z)$

Equating coefficients on both sides, we have,

$(n-1)f_n=n^2(4n-3)f_{n-1}+4n^2(n-1)^2f_{n-2}$ with $f_0=0$ and $f_1=2$.

90 deg. Clockwise (Counter-Clockwise) symmetry - Let $c_n$ denote the number of colorings of a $(2n+1) \times (2n+1)$ grid that are fixed by rotating the grid by 90 degrees (270 degrees). We proceed like in the previous post.

$c_n=\displaystyle[x_1^2x_2^2\cdots x_{2n+1}^2 y_1^2y_2^2\cdots y_{2n+1}^2]\prod_{i=1}^{n+1}\prod_{j=1}^n(1+x_iy_j\cdot x_jy_{2n+1-i}\cdot x_{2n+1-i}y_{2n+1-j}\cdot x_{2n+1-j}y_i)$

Now let $t_i=x_iy_ix_{2n+2-i}y_{2n+2-i}$. Then,

$\begin{align}
c_n&=\displaystyle[t_1^2t_2^2\cdots t_n^2t_{n+1}]\prod_{i=1}^{n+1}\prod_{j=1}^n(1+t_it_j)\\
&=[t_1^2t_2^2\cdots t_n^2t_{n+1}](1+(t_1+t_2+\cdots+t_n)t_{n+1}+O(t_{n+1}^2))\prod_{i=1}^n\prod_{j=1}^n(1+t_it_j)\\
&=n\cdot[t_1^2t_2^2\cdots t_{n-1}^2t_n]\prod_{i=1}^n\prod_{j=1}^n(1+t_it_j)\\
\end{align}$

Let $x_n=[t_1^2t_2^2\cdots t_{n-1}^2t_n]\prod_{i=1}^n\prod_{j=1}^n(1+t_it_j)$. We manipulate this expression hoping to get a recurrence.

$\begin{align}
x_n&=\displaystyle[t_1^2t_2^2\cdots t_{n-1}^2t_n]\prod_{i=1}^n\prod_{j=1}^n(1+t_it_j)\\
&=[t_1^2t_2^2\cdots t_{n-1}^2t_n](1+t_1t_n)^2(1+t_2t_n)^2\cdots(1+t_{n-1}t_n)^2(1+t_n^2)\prod_{i=1}^{n-1}\prod_{j=1}^{n-1}(1+t_it_j)\\
&=[t_1^2t_2^2\cdots t_{n-1}^2t_n](1+t_1t_n)^2(1+t_2t_n)^2\cdots(1+t_{n-1}t_n)^2\prod_{i=1}^{n-1}\prod_{j=1}^{n-1}(1+t_it_j)\\
&=[t_1^2t_2^2\cdots t_{n-1}^2t_n](1+(t_1+t_2+\cdots+t_{n-1})t_n+O(t_n^2))^2\prod_{i=1}^{n-1}\prod_{j=1}^{n-1}(1+t_it_j)\\
&=[t_1^2t_2^2\cdots t_{n-1}^2t_n](1+2(t_1+t_2+\cdots+t_{n-1})t_n+O(t_n^2))\prod_{i=1}^{n-1}\prod_{j=1}^{n-1}(1+t_it_j)\\
&=2(n-1)[t_1^2t_2^2\cdots t_{n-2}^2t_{n-1}]\prod_{i=1}^{n-1}\prod_{j=1}^{n-1}(1+t_it_j)\\
&=2(n-1)x_{n-1}\\
\end{align}$

But it's easy to see that $x_1=\displaystyle[t_1]\prod_{i=1}^1\prod_{j=1}^1(1+t_it_j)=[t_1](1+t_1^2)=0$

Therefore, like in the case of reflections, there are no colorings that are fixed by clockwise (anti-clockwise) rotations.

Now that we have all of the above, we can use Burnside's Lemma to solve the problem for grids of odd size. 

Hope you enjoyed this. See ya later with another interesting topic.


Until then
Yours Aye
Me

A Quest for Pi

For $m<n$, it is easy to show that,

$\displaystyle\frac{x^m}{(x-a_0)(x-a_1)\cdots(x-a_{n-1})}=\sum_{k=0}^{n-1}\frac{a_k^m}{x-a_k}\prod_{i=0\\i \neq k}^{n-1}(a_k-a_i)^{-1}$

Let's choose, $a_k=\eta^k$ where $\eta$ is the $n^{\text{th}}$ root of unity. Then we have,

$\displaystyle \frac{x^m}{1-x^n}=\sum_{k=0}^{n-1}\frac{1}{n\eta^{k(n-1)}}\frac{\eta^{km}}{\eta^k-x}=\frac{1}{n}\sum_{k=0}^{n-1}\frac{\eta^{km}}{1-x/\eta^k}$

Integrating $x$ in the above expression between $0$ and $t$ and comparing the real part, we get,

$\begin{align}
\displaystyle \sum_{j \geq 0}\frac{t^{jn+m}}{jn+m}&=\frac{1}{n}\sum_{k=0}^{n-1}\cos(2\pi km/n)\log(1-2\cos(2\pi k/n)t+t^2)^{-1/2}\\ &\quad +\frac{1}{n}\sum_{k=0}^{n-1}\sin(2\pi km/n)\tan^{-1}\left(\frac{\sin(2\pi k/n)}{t^{-1}-\cos(2\pi k/n)}\right)
\end{align}$

Comparing the Imaginary parts, we also get,

$\displaystyle \sum_{k=0}^{n-1}\cos(2\pi km/n)\tan^{-1}\left(\frac{\sin(2\pi k/n)}{t^{-1}-\cos(2\pi k/n)}\right)=\sum_{k=0}^{n-1}\sin(2\pi km/n)\log(1-2\cos(2\pi k/n)t+t^2)^{-1/2}$

Using the properties of Conjugates, we also have,

$\begin{align}
\displaystyle \sum_{j \geq 0}\frac{t^{jn+m}}{jn+m}&=\frac{1}{n}\log(1-t)^{-1}+\mathbf{1}_{n \text{ even}}\frac{(-1)^m}{n}\log(1+t)^{-1}\\ & \quad +\frac{2}{n}\sum_{k=1}^{\lceil n/2 \rceil-1}\sin(2\pi km/n)\tan^{-1}\left(\frac{\sin(2\pi k/n)t}{1-\cos(2\pi k/n)t}\right)\\ & \quad +\frac{2}{n}\sum_{k=1}^{\lceil n/2 \rceil-1}\cos(2\pi km/n)\log(1-2\cos(2\pi k/n)t+t^2)^{-1/2}
\end{align}$

On the other hand, using Vandermonde matrix or otherwise,

$\displaystyle \frac{1}{1-x/\eta^m}=\sum_{k=0}^{n-1}\eta^{-km}\frac{x^k}{1-x^n}$

Integrating the above,

$\displaystyle \log(1-\eta^{-m}t)^{-1}=\sum_{k=1}^n\eta^{-mk}P_{n,k}(t)$ where $\displaystyle P_{n,k}(t)=\sum_{j \geq 0}\frac{t^{nj+k}}{nj+k}$

Comparing real and imaginary parts, we have,

$\displaystyle \log(1-2\cos(2\pi m/n)t+t^2)^{-1/2}=\sum_{k=1}^n\cos(2\pi km/n)P_{n,k}(t)$

$\displaystyle \tan^{-1}\left(\frac{\sin(2\pi m/n)}{t^{-1}-\cos(2\pi m/n)}\right)=\sum_{k=1}^n\sin(2\pi km/n)P_{n,k}(t)$

The second one here is more interesting. Using $t=\cos(2\pi m/n)$, we get infinitely many representations for $\pi$.

$\displaystyle (n-4m)\frac{\pi}{2n}=\sum_{k=1}^n\sin(2\pi km/n)P_{n,k}(\cos(2\pi m/n))$

For example, with $m=1$ and $n=3$, we get a Bailey–Borwein–Plouffe type formula (BBP formula) for $\pi$ in Base-64 which can be verified with Wolfram Alpha.

$\displaystyle \frac{\pi}{(\sqrt{3}/2)^3}=\sum_{k=0}^\infty\left[\frac{1}{64^k}\left(\frac{4}{6k+1}+\frac{2}{6k+2}-\frac{2^{-1}}{6k+4}-\frac{4^{-1}}{6k+5}\right)\right]$

(or)

$\displaystyle \frac{32\pi}{3\sqrt{3}}=\sum_{k=0}^\infty\left[\frac{1}{64^k}\left(\frac{16}{6k+1}+\frac{8}{6k+2}-\frac{2}{6k+4}-\frac{1}{6k+5}\right)\right]$


Until then
Yours Aye
Me

Pólya's Enumeration Theorem with Applications - Part III

Ever since I first encountered Polya's Enumeration theorem (PET) to solve a problem from Project Euler, it has continued to fascinate on every occasion. I already wrote about them sometime back. Annoying Precision has a beautiful post about PET which I highly recommend. In this post, we are going to see three problems and the associated Cycle Indices. The first two are quite common whereas the third not so much.

The common setup for all the three problems is that we are given an infinite number of balls of $m$ different colors (all identical otherwise). Let the colors be $c_1,c_2,c_3,\cdots,c_m$. Note that we are going for the 'inventory' or the generating function $F(c_1,c_2,c_3,\cdots,c_m)$ rather than just the count.

Case 1 Well..in the first case, given the setup above, we are interested in choosing $n$ balls. This is easy as PET tells us that,

$F_1(c_1,c_2,\cdots,c_m)=Z_{S_n}(a_1,a_2,a_2,\cdots,a_n)$

where $S_n$ denotes the Symmetric group of order $n$ and $a_k=c_1^k+c_2^k+\cdots+c_m^k$.

For example, with $n=3$ and four colors black, white, green and red i.e. $c_1=W$, $c_2=B$, $c_3=R$ and $c_4=G$,

$\begin{align}
F_1(B,G,R,W)&=Z_{S_3}(B+G+R+W,B^2+G^2+R^2+W^2,B^3+G^3+R^3+W^3)\\
&=B^3 + B^2 G + B G^2 + G^3 + B^2 R + B G R + G^2 R + B R^2 +
 G R^2 + \\
&\quad R^3 + B^2 W + B G W + G^2 W + B R W + G R W + R^2 W + B W^2 + \\
&\quad G W^2 + R W^2 + W^3\\
\end{align}$

The use of Symmetric Group $S_n$ as we are 'choosing' the balls rather than 'ordering' them. Their Cycle indices are given by

$\displaystyle Z(S_n)=\sum_{i_1+2i_2+3i_3+\cdots+ni_n=n}\prod_{k=1}^n\frac{a_k^{i_k}}{k^{i_k}i_k!}$

The generating function of the same is given by,

$\displaystyle \sum_{n=0}^\infty Z(S_n)t^n=\text{exp}\left(a_1\frac{t}{1}+a_2\frac{t^2}{2}+a_3\frac{t^3}{3}+a_4\frac{t^4}{4}+\cdots\right)$


Case 2 In the second case, we are choosing $n$ balls such that each color can be used at most once. Now, with PET, we have,

$F_2(c_1,c_2,\cdots,c_m)=Z_{A_n}(a_1,a_2,a_2,\cdots,a_n)$

where $A_n$ denotes the Alternating group of order $n$.

For example, with $n=3$ and four colors black, white, green and red i.e. $c_1=W$, $c_2=B$, $c_3=R$ and $c_4=G$,

$\begin{align}
F_2(B,G,R,W)&=Z_{A_3}(B+G+R+W,B^2+G^2+R^2+W^2,B^3+G^3+R^3+W^3)\\
&=B G R + B G W + B R W + G R W\\
\end{align}$

$\displaystyle Z(A_n)=\sum_{i_1+2i_2+3i_3+\cdots+ni_n=n}(-1)^{i_2+i_4+\cdots}\prod_{k=1}^n\frac{a_k^{i_k}}{k^{i_k}i_k!}$

The corresponding generating function is given by,

$\displaystyle \sum_{n=0}^\infty Z(A_n)t^n=\text{exp}\left(a_1\frac{t}{1}-a_2\frac{t^2}{2}+a_3\frac{t^3}{3}-a_4\frac{t^4}{4}+\cdots\right)$

Case 3 In the last case, we will be choosing $n$ balls such that histogram of colors follows $\mathbf{n}=\{j_1,j_2,j_3,\cdots\}=\{1^{n_1},2^{n_2},3^{n_3},\cdots\}$ with $1\leq j_1 \leq j_2 \leq \cdots$ and $n=j_1+j_2+\cdots=n_1+2n_2+3n_3+\cdots$.

Computing the cycle index is tricky in itself. In fact, it is not even clear to me what 'group' are we dealing with here.

Any way to compute the cycle index, consider the polynomials $P_{\{j_1,j_2,j_3,\cdots,j_m\}}(a_1,a_2,a_3,\cdots)$ defined by the recursion (without the arguments $a_i$'s for clarity)

$\displaystyle P_{\{j_1,j_2,\cdots,j_m\}}=a_{j_m}P_{\{j_1,j_2,\cdots,j_{m-1}\}}-\sum_{i=1}^{m-1}P_{\{j_1,j_2,\cdots,j_{i-1},j_i+j_m,j_{i+1},\cdots,j_{m-1}\}}$

with the boundary condition $P_{\{\}}=1$.  Now the Cycle index is given by,

$\displaystyle Z(G_{\mathbf{n}})=\frac{1}{n_1!n_2!n_3!\cdots}P_{\{j_1,j_2,j_3,\cdots\}}(a_1,a_2,a_3,a_4,\cdots)$

The Generating function of Cycle indices of all 'histograms' can be shown to be,

$\displaystyle \sum_{\mathbf{n}} Z(G_{\mathbf{n}})\mathbf{t}^{\mathbf{n}}=\text{exp}\left(a_1C_1-a_2C_2+a_3C_3-a_4C_4+\cdots\right)$

where the sum on the LHS runs through all multisets $\mathbf{n}=\{1^{n_1},2^{n_2},3^{n_3}\cdots\}$ with $\mathbf{t}^\mathbf{n}=t_1^{n_1}t_2^{n_2}t_3^{n_3}\cdots$ and $C_n(t_1,t_2,t_3,\cdots)$'s are such that

$\begin{align}
\displaystyle C_n(t_1,t_2,t_3,\cdots)&=[z^n]\log (1-t_1z+t_2z^2-t_3z^3+\cdots)^{-1}\\
&=\sum_{k=1}^n\frac{1}{k}\hat{B}_{n,k}(t_1,-t_2,\cdots,(-1)^{k-1}t_{n-k+1})\\
&=\sum_{n_1+2n_2+3n_3+\cdots=n}(-1)^{n_2+n_4+\cdots}\binom{n_1+n_2+\cdots}{n_1,n_2,\cdots}\frac{\mathbf{t}^\mathbf{n}}{n_1+n_2+\cdots}\\
\end{align}$

where $\hat{B}_{n,k}(x_1,x_2,\cdots,x_{n-k+1})$ are the ordinary Bell Polynomials and the sum in the last equality runs through all integer partitions of $n$.

Note that (i) at $t_1=t$ and $t_i=0$ for $i>1$, that above reduces to the Case 2 and (ii) at $t_i=t$ for $i>0$ the above reduced to Case 1 as it should.

For example, we again have the four colors as in the previous cases but we choose $\mathbf{n}=\{1,1,2\}$. That is, we choose four balls of which two balls are of different colors and the remaining two balls are of the same color but that which different from the first two.

$\begin{align}
F_3(B,G,R,W)&=Z_{G_{\{1,1,2\}}}(B+G+R+W,B^2+G^2+R^2+W^2,B^3+G^3+R^3+W^3,B^4+G^4+R^4+W^4)\\
&=B^2 G R + B G^2 R + B G R^2 + B^2 G W + B G^2 W + B^2 R W + \\
&\quad G^2 R W +
 B R^2 W + G R^2 W + B G W^2 + B R W^2 + G R W^2\\
\end{align}$

where I've used $G_{\mathbf{n}}$ to denote the Group we are dealing with.

Now that we come to the end of the post, we make note of something beautiful in case it went unnoticed so far. In fact, it is also the reason why we looked at these three cases.

If we just set aside the notion of PET for a moment, what we have done amounts to expressing the homogenous symmetric functions, elementary symmetric functions and the monomial symmetric functions (respectively) in terms of the powersum symmetric functions in variables $B,G,R,$ and $W$. This reveals a nice connection between Cycle indices (along with their resp. generating functions) and ring of Symmetric functions.

In addition, All of the above have interesting applications in subset/multiset sums modulo a parameter, number of solutions to (multipartite) linear Diophantine equations, multipartite partitions and by extension to un-ordered factorizations which I invite the readers to explore. I found these extremely interesting and instructive to work all these out. Hope you enjoyed too.

(untidy) Mathematica code:

Clear["Global`*"];

MonomSymPoly[l_] := 

  MonomSymPoly[l] = 

   Which[Length[l] <= 0, 1, True, 

    Subscript[p, Last[l]] MonomSymPoly[Most[l]] - 

     Sum[MonomSymPoly[

       Sort[Most[l] + 

         Last[l] Table[Boole[i == j], {j, Length[l] - 1}]]], {i, 

       Length[l] - 1}]];

MultiSetCycIndex[l_] := 

  Factor[Expand[

    MonomSymPoly[l]/Times @@ Factorial[Length /@ Split[l]]]];

(*Sum[BellY[n,k,Table[Power[-1,j-1]Factorial[ j-1]Subscript[x, \

j],{j,n-k+1}]],{k,n}]/Factorial[n]*)

Li = {1, 2, 2};

Li = Sort[Li];

MultiSetCycIndex[Li]

Expand[% /. Subscript[p, j_] :> W^j + B^j + R^j + G^j]

Print["Hi"];

n = 4;

m = Total[Li];

Coeffs = Table[{Subscript[t, j], 0, 0}, {j, m}];

Do[

  Coeffs[[j, 3]] += 1;

  , {j, Li}];

(*Exp[Sum[Factor[Sum[Factorial[k-1]BellY[n,k,Table[Power[-1,j-1]\

Factorial[j]Subscript[t, j],{j,n-k+1}]],{k,n}]Subscript[p, \

n]/Factorial[n]],{n,m}]]//TraditionalForm*)

(*Sum[Factorial[k-1]BellY[n,k,Table[Power[-1,j-1]Factorial[j]\

Subscript[t, j],{j,n-k+1}]],{k,n}]/Factorial[n]*)

SeriesCoefficient[

 Exp[Sum[Sum[

     Factorial[k - 1] BellY[n, k, 

       Table[Power[-1, j - 1] Factorial[j] Subscript[t, j], {j, 

         n - k + 1}]], {k, n}] Power[-1, n - 1] Subscript[p, n]/

     Factorial[n], {n, m}]], Evaluate @@ Sequence[Coeffs]]

(*SeriesCoefficient[Exp[Sum[Factor[SeriesCoefficient[-Log[1+Sum[\

Subscript[t, i]Power[-x,i],{i,m}]],{x,0,k}]Subscript[p, \

k]],{k,m}]],Evaluate@@Sequence[Coeffs]]*)

SeriesCoefficient[

 Exp[Sum[Factor[

    SeriesCoefficient[

      Log[Power[

        1 + Sum[Subscript[t, i] Power[-x, i], {i, m}], -1]], {x, 0, 

       k}] Power[-1, k - 1] Subscript[p, k]], {k, m}]], 

 Evaluate @@ Sequence[Coeffs]]

Expand[% /. Subscript[p, j_] :> W^j + B^j + R^j + G^j]

Sum[Factor[

  SeriesCoefficient[

    Log[Power[1 + Sum[Subscript[t, i] Power[-x, i], {i, m}], -1]], {x,

      0, k}] Power[-1, k - 1] Subscript[p, k]], {k, m}]

Series[Log[

  Power[1 + Sum[Subscript[t, i] Power[-x, i], {i, m}], -1]], {x, 0, m}]

Series[Exp[Sum[Subscript[a, k] Power[t, k]/k, {k, 4}]], {t, 0, 4}]

Until then
Yours Aye
Me

A note on Partial fractions

I recently found a nice identity about Partial fractions on Interesting identities about Partial fractions. The same could also be found, albeit on a different form, on A note on partial fractions. So this post is actually a note on 'A note on Partial fractions'.

It says, if $u+v=$, then for integers $m,n \geq 1$,

$\displaystyle \frac{1}{u^mv^n}=\sum_{k=0}^{m-1}\frac{\binom{n-1+k}{k}}{u^{m-k}}+\sum_{l=0}^{n-1}\frac{\binom{m-1+l}{l}}{v^{n-l}}$

The wordpress page proves this with L'Hospital's rule whereas the webpage proves this with induction. But for someone interested in Probability, two-things-adding-upto-1, the binomial coefficient and then it immediately hit me.

It's actually the 'Problem of Points', one of the classical problems of probability where two friends $A$ and $B$ play a series of games where $A$ needs $m$ wins and $B$ needs $n$ wins to win the match. Assume $A$ has a probability $p$ of winning a each game whereas $B$'s probability is $q=1-p$.

Obviously,

$\mathbb{P}(\text{A wins the match})+\mathbb{P}(\text{B wins the match})=1$

Let's compute the probability of $A$ winning the match. As Pascal reasoned, $A$ needs a minimum of $m$ games and a maximum of $m+n-1$ games to win the match. By the law of total probability,

$\begin{align}
\displaystyle\mathbb{P}(\text{A wins the match})&=\sum_{k=m}^{m+n-1}\mathbb{P}(\text{A wins the match on the }k^{th}\text{ game})\\
&=\sum_{k=m}^{m+n-1}\mathbb{P}(\text{A won }m-1\text{ games out of }k-1\text{ games})\cdot\mathbb{P}(\text{A wins the }k^{th}\text{ game})\\
&=\sum_{k=m}^{m+n-1}\binom{k-1}{m-1}p^{m-1}q^{k-m}\cdot p\\
&=\sum_{k=0}^{n-1}\binom{m+k-1}{k}p^mq^k
\end{align}$

This is 'textbook' Negative Binomial.

Interchanging $m$ with $n$ and $p$ with $q$ gives the probability of $B$ winning the match. Therefore, we finally have,

$\displaystyle \sum_{k=0}^{n-1}\binom{m+k-1}{k}p^mq^k+\sum_{k=0}^{m-1}\binom{n+k-1}{k}q^np^k=1$

Dividing both sides by $p^mq^n$, gives the identity quoted at the start of this post. Tada!!

Note that Fermat's reasoning for the Problem of points produces the following equation. Though correct, it is not particularly useful as a partial fraction decomposition.

$\displaystyle \sum_{k=0}^{n-1}\binom{m+n-1}{m+k}p^{m+k}q^{n-1-k}+\sum_{k=0}^{m-1}\binom{n+m-1}{n+k}q^{n+k}p^{m-1-k}=1$

Hope you enjoyed this.


Until then
Yours Aye
Me

A note on Berlekamp Massey Algorithm

Let's say we are given a sequence of $2n$ terms and were told that these are generated by a linear recurrence of order $n$ with the first $n$ terms of the sequence being the initial terms. More specifically, given $n$ terms, we want a recurrence that generates the last $\lfloor \frac{n}{2}\rfloor$ terms with the first $\lceil \frac{n}{2} \rceil$ terms.

The problem of finding the recurrence given the terms can be solved by Berlekamp Massey (BM) Algorithm. Though there are not many sources that I could find for solving recurrences with BM, I found the CodeForces entry that explains much of the details. Yet another source that I found very interesting is Berlekamp Massey algorithm.

To make the algorithm more clear than the quoted page, let's take an example. Let the sequence of terms be $[1,2,4,10,24,50,124]$.

Looking at (only) the first three terms, we make a first guess that

$\text{Guess 1:}\quad a_n=2a_{n-1}$ (or) $a_n-2a_{n-1}=0$

But when we include the fourth term, we run into a problem. We have a discrepancy.

$\Delta_4=\text{Actual}-a_4=10-2a_3=10-8=2$

We stare into the (only) first four terms of the sequence and come with another guess.

$\text{Guess 2:}\quad a_n=2a_{n-1}+a_{n-2}$ (or) $a_n-2a_{n-1}-2a_{n-2}=0$

And this correctly generates the fourth term.

(two tiny details at this point for the curious minded - starts here...)
First, someone could've guessed $a_n=2a_{n-1}+2a_{n-3}$. That works too. Continuing with the BM algorithm, we may end up with different answers but both recurrences will correctly generate the last three terms using the first four terms. Try it!

Second, though our second guess correctly predicts $a_4$, it fails for $a_3$. This is not a problem as Flawed guesses, especially for the first terms in the first half, gets 'course-corrected' at later stage of the algorithm. In the final algorithm, we only be making guesses, that too trivial guesses, at the initial stage.
(... and ends here)

We try it for the next term and we see the recurrence correctly generates fifth term. We're good. We try it for the next term and we run into problem. We have a discrepancy of $8(=58-50)$.

But this time we don't have to guess anything more. Note that we re-write our 'Guess 2' as

$\text{Guess 2:}\quad a_n=2a_{n-1}+a_{n-2}+k(a_{n-2}-2_{n-3})$

for any real $k$.

Where did the terms in the bracket come from? It our first guess shifted by two terms.

Why have we shifted it by two terms? Because we are currently finding a discrepancy at the sixth term and our last guess created a discrepancy at fourth term.

Will this not affect the previous terms? That is the key of BM algorithm. The terms in the bracket will become zero when we evaluate it for earlier terms.

For example, when we evaluate Guess 2 for $a_5$, the terms in the bracket will be $a_3-2a_2$ which is actually zero. In fact, only because it was zero, we made it our first guess.

Now that we have an additional factor of $k$ in our 'Guess 2' now, we can use it to our advantage.  Choosing $k=-4$, we remove the discrepancy at the sixth position.

From now on, we can continue using the same idea. Anytime we face a discrepancy, we can choose an right multiple (based on the current and previous deviation) of our previous guess (shifted by the right amount) to correct that discrepancy. And this continues for all elements in the given sequence, thus completing the algorithm.

This type of problem comes every now and then in Competitive programming. For a given problem, we can quickly code a brute force method that can churn out the first few terms and we would be interested in finding a linear recurrence if there is one.

If we have the recurrence, finding the $n$th term (even for $n$ of order $10^{10^7}$) of the recurrence can be found in a matter of seconds. For example, we have this wonderful Codechef Editorial that explains how to do that.

Here is an (over)commented Python code of the BM algorithm based on what I wrote above.

from timeit import default_timer

start = default_timer()

def berlmassey(l0, p0):
     # Input: list of length 2n and a prime p

     # Output: Recurrence of length n that takes the first n terms as initial values and generates the next n (and possibly the further) terms

     l, p = l0[:], p0
     if len(l) <= 1:
          return l
     if l[0] != 1:          # WLOG, first element of the list should be 1
          return []
     old_rec = []          # there is no recurrence to start with
     v = 1          # deviation of the (empty) recurrence that occurs in the first element
     old_pos = 0          # trivially the list was fine until the first element
     curr_rec = [l[1] % p]          # we guess something that's trivially true
     curr_pos = 1          # starting with the second element of the list
     while curr_pos < len(l):
          temp = 0
          for i in range(len(curr_rec)):          # check what the current recurrence gives for the next element
               temp += curr_rec[i] * l[curr_pos - i - 1]
               temp %= p
          temp = (l[curr_pos] - temp) % p          # is there a deviation?
          if temp == 0:          # If no deviation, current recurrence is still good and hence moving on..
               curr_pos += 1
               continue
          temp_curr_rec = curr_rec[:]          # just a temporary storage
          v = temp * pow(v, -1, p) % p          # normalizing factor based on current and last deviation
          old_rec = [p - x for x in old_rec]          # re-strucuturing the last successful recurrence
          old_rec = [1] + old_rec
          old_rec = [0] * (curr_pos - old_pos - 1) + old_rec
          if len(old_rec) > len(curr_rec):          # making last successful and current recurrences the same length
               curr_rec = curr_rec + [0] * (len(old_rec) - len(curr_rec))
          else:
               old_rec = old_rec + [0] * (len(curr_rec) - len(old_rec))
          curr_rec = [(curr_rec[i] + v * old_rec[i]) % p for i in range(len(curr_rec))]          # updating the current recurrence with the normalizing constant and the (old) last successful recurrence
          old_rec = temp_curr_rec[:]          # the recurrence we started with becomes the last successful recurrence
          v = temp          # the deviation that occured at this position
          old_pos = curr_pos          # this position becomes the last deviated position
          curr_pos += 1          # On to the next iteration..
     return curr_rec


print(berlmassey([1, 2, 4, 10, 24, 50, 124], 1000000007))
print(default_timer() - start)

Finally, a couple of points regarding the code. As there are divisions involved, I have used a prime and made everything with modulo operations. We can also use exact divisions if there is a need. Hope you enjoyed it.


Until then
Yours Aye
Me

A Binary Grid Counting Problem - Part I

This post is inspired by the A collection of grid counting problems. The entire wordpress page is a great read and I thoroughly enjoy it every time I visit.

In Example 5 of the quoted page, the author discusses the number of ways of coloring a $n \times n$ grid in which every row and every column has exactly two black squares. In the example that's immediately next, he discusses a counting problem that uses one of favorite results in combinatorics - the Burnside's Lemma.

That got me thinking about the number of unique (upto rotations and reflections) colorings of a $2n \times 2n$ grid with exactly two black squares in each row and column.

We'll also be using Burnside's lemma for this task. We'll be considering the eight symmetries of the square and count how many colorings are fixed by each of these symmetries. We'll progressing on the increasing order of difficulty. Let's begin.

Identity - Let $t'_n$ denote total number of such colorings on an $n \times n$ grid. This is essentially what the wordpress page was doing. In fact we have A001499 that gives us a recurrence. With $t'_0=1$, we have

$t'_n=n(n-1)t'_{n-1}+\frac{1}{2}n(n-1)^2t'_{n-2}$

As we are only interested in $2n \times 2n$ grid, let's define $t_n=t'_{2n}$. With $t_0=1$, this follows the recurrence,

$\displaystyle t_n=\frac{n(2n-1)^2}{2n-3}(8n^2-16n+7)t_{n-1}-2n(n-1)^2(2n-1)^2(2n-3)t_{n-2}$

(or)

$\displaystyle \frac{t_n}{2n-1}=n(2n-1)(8n^2-16n+7)\frac{t_{n-1}}{2n-3}-2n(n-1)^2(2n-1)(2n-3)(2n-5)\frac{t_{n-2}}{2n-5}$

Diagonal (Anti-Diagonal) symmetry - Let $d'_n$ denote the required colorings that are symmetric along the main diagonal. Again, we have A000986 that gives us the following recurrence. Starting with $d'_0=1$,

$d'_n=2(n-1)d'_{n-1}-(n-1)(n-3)d'_{n-2}-\frac{1}{2}(n-1)(n-2)(n-3)d'_{n-4}$

Again, let's define $d_n=d'_{2n}$ which is the number of required colorings with diagonal (anti-diagonal) symmetry in a $2n \times 2n$ grid. With $d_0=1$,

$\displaystyle d_n=\frac{2n-1}{2n-3}(8n^2-16n+7)d_{n-1}-\frac{(n-1)(2n-1)}{2n-5}(16n^3-104n^2+230n-173)d_{n-2}-2(n-1)(n-2)(2n-1)(8n^2-38n+41)d_{n-3}-2(n-1)(n-2)(n-3)(2n-1)(2n-3)(2n-7)d_{n-4}$

(or)

$\displaystyle \frac{d_n}{2n-1}=(8n^2-16n+7)\frac{d_{n-1}}{2n-3}-(n-1)(16n^3-104n^2+230n-173)\frac{d_{n-2}}{2n-5}-2(n-1)(n-2)(2n-7)(8n^2-38n+41)\frac{d_{n-3}}{2n-7}-2(n-1)(n-2)(n-3)(2n-3)(2n-7)(2n-9)\frac{d_{n-4}}{2n-9}$

Horizontal (Vertical) reflection - Let $r_n$ denote the number of colorings that are fixed by Horizontal reflection on a $2n \times 2n$ grid. As we need the grid to be the same after reflecting it horizontally (vertically), we only have to fill up the top $n \times 2n$ grid. We have to make sure that each row has exactly two black squares and each column has exactly one black square in the $n \times 2n$ grid.

Proceeding in the same way as in the wordpress page, we can see that,
$\begin{align}
r_n &= [x_1^2x_2^2\cdots x_n^2y_1y_2\cdots y_{2n}]\prod_{i=1}^n\prod_{j=1}^{2n}(1+x_iy_j)\\
&= [x_1^2x_2^2\cdots x_n^2y_1y_2\cdots y_{2n}]\prod_{j=1}^{2n}(1+(x_1+x_2+\cdots+x_n)y_j+O(y_j^2))\\
&= [x_1^2x_2^2\cdots x_n^2](x_1+x_2+\cdots+x_n)^{2n}=2^{-n}(2n)!\\
\end{align}$

So, we have the recurrence $r_n=n(2n-1)r_{n-1}$ with $r_0=1$.

180 deg. rotation symmetry - Let $f_n$ denote the number of colorings of a $2n \times 2n$ grid that are fixed by rotating the grid by 180 degrees. Again like the previous case, we only have to fill up the top $n \times 2n$ grid. We still make sure that each row has exactly two black squares but unlike the previous case, we need to ensure that the total number of black squares in the $j$-th column and $(2n+1-j)$-th column is two.
$\begin{align}
f_n&=[x_1^2x_2^2\cdots x_n^2y_1y_2\cdots y_{2n}]\prod_{i=1}^n\prod_{j=1}^{2n}(1+x_iy_j)\\
&=[x_1^2x_2^2\cdots x_n^2t_1^2t_2^2\cdots t_n^2]\prod_{i=1}^n\prod_{j=1}^n(1+x_it_j)^2\\
&=[x_1^2x_2^2\cdots x_n^2t_1^2t_2^2\cdots t_n^2]\prod_{j=1}^n(1+(x_1+x_2+\cdots +x_n)t_j+(x_1x_2+x_2x_3+\cdots+x_{n-1}x_n)t_j^2+O(t_j^3))^2\\
&=\displaystyle[x_1^2x_2^2\cdots x_n^2]\Biggl(\biggl(\sum_{i=1}^n x_i\biggl)^2+2\sum_{1\leq i<j\leq n} x_ix_j\Biggl)^n\\
&=\displaystyle[x_1^2x_2^2\cdots x_n^2]\Biggl(2\biggl(\sum_{i=1}^n x_i\biggl)^2-\sum_{i=1}^n x_i^2\Biggl)^n\\
&=\displaystyle[x_1^2x_2^2\cdots x_n^2]\sum_{k=0}^n(-1)^k\binom{n}{k}\biggl(\sum_{i=1}^n x_i^2\biggl)^k 2^{n-k}\biggl(\sum_{i=1}^n x_i\biggl)^{2(n-k)}\\
&=\displaystyle\sum_{k=0}^n(-1)^k\binom{n}{k}\frac{n!}{(n-k)!} 2^{n-k}\frac{(2n-2k)!}{2^{n-k}}\\
&=\displaystyle\sum_{k=0}^n(-1)^k \binom{n}{k}^2k!(2n-2k)! \\
\end{align}$

where we have used $t_i=y_i=y_{2n+1-i}$ in the second equality.

Even though we have a closed form, I still wanted to find the recurrence relation. Taking a cue from A001499, I found the following generating function for $f_n$.

$\displaystyle \sum_{n=0}^\infty f_n \frac{z^n}{n!^2}=\frac{e^{-z}}{\sqrt{1-4z}}$

This helped me in finding A296618 which has a recurrence relation for the EGF coefficients of the above function which are actually $\frac{f_n}{n!}$. From there it was straightforward to see that with $f_0=1$,

$f_n=n(4n-3)f_{n-1}+4n(n-1)^2f_{n-2}$

90 deg. Clockwise (Counter-Clockwise) symmetry - Let $c_n$ denote the number of colorings of a $2n \times 2n$ grid that are fixed by rotating the grid by 90 degrees (270 degrees). Now we only have to fill the first quadrant of the grid. This time we have,

$c_n=\displaystyle[x_1^2x_2^2\cdots x_{2n}^2 y_1^2y_2^2\cdots y_{2n}^2]\prod_{i=1}^n\prod_{j=1}^n(1+x_iy_j\cdot x_jy_{2n+1-i}\cdot x_{2n+1-i}y_{2n+1-j}\cdot x_{2n+1-j}y_i)$

Now let $t_i=x_iy_ix_{2n+1-i}y_{2n+1-i}$. Then,

$c_n=\displaystyle[t_1^2t_2^2\cdots t_n^2]\prod_{i=1}^n\prod_{j=1}^n(1+t_it_j)$

I was not able to simplify this in anyway and I tried a bunch of approaches. Finally, I gave up and found the sequence with Mathematica. The values for $n=0,1,2,\cdots$ are $1, 1,2,12,90,810,8940,\cdots$. Unfortunately, this time we don't have a OEIS sequence to directly get a recurrence. After laboring for a day, I found the following recurrence with a lot of paper work. With $c_0=1$,

$c_n=(2n-1)c_{n-1}-(n-1)c_{n-2}+2(n-1)(n-2)c_{n-3}$

Searching for a generating function, we get,

$\displaystyle\sum_{n=0}^\infty c_n\frac{z^n}{n!}=\frac{e^{-z^2/2}}{\sqrt{1-2z}}$

From this we get the following explicit formula.

$\displaystyle c_n=\frac{n!}{2^n}\sum_{k=0}^n \frac{(-2)^k}{k!}\binom{2n-4k}{n-2k}$

With ideas from Algorithms forHolonomic Functions, Algorithmic Manipulations and Transformations of Univariate Holonomic Functions and Sequences and Computing the Generating Function of a series given its First few terms, it seems we can find the above recurrence from the generating function.

Back to the main question of this post, we finally have everything we need to apply Burnside's lemma. Let $A_n$ denote the number of unique (upto rotations and reflections) colorings of a $2n \times 2n$ grid with exactly two black squares in each row and column. Then,

$A_n=\displaystyle\frac{1}{8}(t'_{2n}+2d'_{2n}+2r_n+f_n+2c_n)=\frac{1}{8}(t_n+2d_n+2r_n+f_n+2c_n)$

On a side note, the number of unique (upto rotations and reflections) colorings of a $n \times n$ grid with exactly one black square in each row and column can be found completely using OEIS. The only thing we have to note is that there can be no colorings for vertical or horizontal reflections. If $B_n$ is the required result, then

$\displaystyle B_n=\frac{1}{8}(4\text{ A263685}(n)+2\text{ A000085}(n))$

Hope you enjoyed this. See ya later with another interesting topic.


Until then
Yours Aye
Me

A nice result in Geometry

Euler's theorem in Geometry is a great result that connects the distance between the circumcentre and incentre of a triangle with the circumradius and inradius of the triangle. A similar result, known as Fuss's theorem, holds for bi-centric quadrilaterals. It was surprising when I found a similar relation holds for three mutually tangent circular arcs.

Consider the three blue circles $A$, $B$ and $C$ (with radii $r_A$, $r_B$ and $r_C$ respectively) that are mutually tangent to each other. Let the points of tangency be $T_A$, $T_B$ and $T_C$ which becomes the vertices of our 'triangle of circular arcs'.

We can now draw two circles: a circumcirle (red) that passes through all the tangency points and an incircle (green) that is mutually internally tangent to all the blue circles. Let's position the configuration in the Argand plane such that center of the incircle is at the origin.

Let the centers of the blue circles be $z_A$, $z_B$ and $z_C$, and the center of the circumcircle be $z$.

Now if we let $R$ and $r$ be the circumradius and inradius of our triangle, then
$$\text{Eqn (1): }d^2=R^2-4Rr+r^2$$
where $d=\vert z \vert$ is the distance between the centers of the circumcircle and incircle.

To prove this, Consider $\triangle z_A z_B z_C$ whose side lengths are $r_A+r_B$, $r_B+r_C$ and $r_C+r_A$. It is easy to show that,

$\text{Eqn (2): }R^2=\displaystyle\frac{r_Ar_Br_C}{r_A+r_B+r_C}$ and $\text{Eqn (3): }r^{-1}={r_A}^{-1}+{r_B}^{-1}+{r_C}^{-1}+2R^{-1}$

$(2)$ is obtained by noting that the red circle is the incircle of $\triangle z_A z_B z_C$ and $(3)$ is a direct application of Descartes' theorem.

Now given the vertices of a complex triangle, we can express the incenter of the triangle using its vertices (using this). Now we have,

$\displaystyle z=\sum_{\text{cyc}}\frac{r_B+r_C}{2(r_A+r_B+r_C)}z_A \implies \text{Eqn (4): }2(r_A+r_B+r_C)z=\sum_{\text{cyc}}(r_B+r_C)z_A$

Multiplying $(4)$ with its complex conjugate, we have,

$\text{Eqn (5): }4(r_A+r_B+r_C)^2\vert z \vert^2=\displaystyle\sum_{\text{cyc}}(r_B+r_C)^2\vert z_A \vert^2+\sum_{\text{cyc}}(r_B+r_C)(r_C+r_A)(z_A\bar{z_B}+\bar{z_A}z_B)$

From the figure, we see that $\vert z_B-z_A \vert=r_A+r_B$. Using the fact that the product of a complex number and its conjugate is the square of its magnitude, we can see that,

$\text{Eqn (6): }z_A\bar{z_B}+\bar{z_A}z_B=\vert z_A \vert^2+\vert z_B \vert^2-\vert z_B-z_A\vert^2=\vert z_A \vert^2+\vert z_B \vert^2-(r_A+r_B)^2$

Using $(6)$ in $(5)$, and cancelling the common terms, we get,

$\text{Eqn (7): }2(r_A+r_B+r_C)\vert z \vert^2=\displaystyle\sum_{\text{cyc}}\vert z_A \vert^2(r_B+r_c)-\prod_{\text{cyc}}(r_A+r_B)$

From the figure above, we have $\vert z_x \vert=r+r_x$ for $x \in \{A,B,C\}$. Using this we get,

$\begin{align}
\displaystyle\text{Eqn (8): }\sum_{\text{cyc}}\vert z_A \vert^2(r_B+r_c)&=\sum_{\text{cyc}}(r^2+2rr_A+r_A^2)(r_B+r_C) \\
&= 2r^2(r_A+r_B+r_C)+4r(r_Ar_B+r_Br_C+r_Cr_A)+\sum_{\text{cyc}}r_A^2(r_B+r_C) \\
\end{align}$

Using $(8)$ in $(7)$, we get,

$\text{Eqn (9): }\displaystyle2(r_A+r_B+r_C)\vert z \vert^2=2r^2(r_A+r_B+r_C)+4r(r_Ar_B+r_Br_C+r_Cr_A)+\sum_{\text{cyc}}r_A^2(r_B+r_C)-\prod_{\text{cyc}}(r_A+r_B)$

Now the second term of $(9)$ can be simplified as
$\begin{align}
\text{Eqn (10): }r(r_Ar_B+r_Br_C+r_Cr_A)&=rr_Ar_Br_C(r_A^{-1}+r_B^{-1}+r_C^{-1})\\
&=rR^2(r_A+r_B+r_C)(r^{-1}-2R^{-1})\\
&=R(R-2r)(r_A+r_B+r_C)
\end{align}$
where we have used $(2)$ and $(3)$ in the second equality.

The last two terms of $(9)$ can be expanded and most of the terms gets cancelled giving,

$\displaystyle\text{Eqn (11): }\sum_{\text{cyc}}r_A^2(r_B+r_C)-\prod_{\text{cyc}}(r_A+r_B)=-2r_Ar_Br_C=-2R^2(r_A+r_B+r_C)$
where we have used $(2)$ in the last equality.

Using $(10)$ and $(11)$ in $(9)$ and dividing both sides by $2(r_A+r_B+r_C)$ , we finally have $(1)$.

From $(1)$, it follows that, $R\geq (2+\sqrt{3})r$ which holds with equality only when $r_A=r_B=r_C$.

This can probably be proved with geometry but I couldn't get my head around it.

Finally, before we finish this discussion, we also note that we can completely parametrize $(1)$ for non-co-prime integer solutions. For integer $a,b$ such that $\gcd(a,b)=1$, $3 \nmid b$ and $b \leq \sqrt{3}a$

(i) If $a\not\equiv b \pmod{2}$, then $R=3a^2+4ab+b^2$, $r=2ab$ and $d=3a^2-b^2$
(ii) If $a\equiv b \pmod{2}$, then $R=(3a^2+4ab+b^2)/2$, $r=ab$ and $d=(3a^2-b^2)/2$

Hope you found this interesting. See ya later.

References:
Tangencies: Three tangent Circles
Efficiently constructing tangent circles
Inversion and applications to Ptolemy and Euler


Until then
Yours Aye
Me