tag:blogger.com,1999:blog-69027612024-03-07T00:08:18.948-08:00My Journey Towards Emptiness!!Algorithms, Theory, Spirituality, Life, Technology, Food and Workout : trying to sort these deterministically in $\Theta(1)$ time (constant time).Vamsi Kundetihttp://www.blogger.com/profile/06920946350286969217noreply@blogger.comBlogger167125tag:blogger.com,1999:blog-6902761.post-33991059431804765712017-07-11T00:18:00.001-07:002017-07-11T12:34:22.288-07:00A matrix formulation of extended Euclid algorithm.I have written about <i>Extended Euclid Algorithm</i> in the past. However more recently while solving a problem I had to express it in a matrix form which I thought is something worthwhile. <br />
<br />
Recall the problem statement: Given two non-zero integers $a,b$ (assume $a\geq b$), the problem is to determine integers $x,y$ such that $\mbox{gcd}(a,b) = a\times x + b \times y$. Notice that the original Euclid's algorithm can compute the gcd but not the coefficients $x,y$ whose linear combination yields the gcd. The idea here is to use the <i>quotient</i> instead of the <i>remainder</i> as in the original algorithm. We express the operands in each iteration a linear combination of the original inputs and use the following to update as follows. The quotient in $i^{th}$ iteration is denoted by $q_i$ leading to the following formulation.<br />
<br />
$$<br />
<br />
\begin{array}{cc}<br />
\mbox{gcd}(a,b) = \mbox{ExGcd}(\left[\begin{array}{c} a \\ b \end{array}\right]) &\\<br />
\left [<br />
\begin{array}{cc}<br />
x_{0,0} & y_{0,0}\\<br />
x_{1,0} & y_{1,0}<br />
\end{array}<br />
\right ] = \left [ \begin{array}{cc}1 & 0\\ 0 & 1\end{array} \right ] & \mbox{initialization } i=0\\<br />
<br />
<br />
\mbox{ExGcd}(\left [\begin{array}{cc} a & b \end{array}\right ] \left [ \begin{array}{cc} x_{0,i} & y_{0,i} \\ x_{1.i} & y_{1,i} \end{array} \right]) = \mbox{ExGcd}(\left[\begin{array}{cc} a & b \end{array}\right] <br />
\left [<br />
\begin{array}{cc}<br />
x_{0,i-1} & y_{0,i-1} \\<br />
x_{1,i-1} & y_{1,i-1} <br />
\end{array}<br />
\right ]<br />
\left [<br />
\begin{array}{cc}<br />
0 & 1 \\<br />
1 & -q_{i-1}<br />
\end{array}<br />
\right ]) & i > 0 \\<br />
<br />
q_i = \mbox{Quotient}(x_{0,i}a+x_{1,i}b, y_{0,i}a+y_{1,i}b) & \\<br />
\end{array}<br />
$$<br />
<br />
<br />
Following is an implementation using <em>std::div</em>:<br />
// Input: a,b <br />
// Output: x, y (y0,y1)<br />
<pre class="brush:cpp;">void ExtendedEuclid(long a, long b, long& y0, long& y1)
{
long x0, x1, t0, t1;
if (a < b) { std::swap(a, b); }
// Initialize:
x0 = 1; x1 = 0; // a = 1 * a + 0 * b
y0 = 0; y1 = 1; // b = 0 * a + 1 * b
std::ldiv_t div_res = std::div(a, b);
while (div_res.rem) {
t0 = y0; t1 = y1;
y0 = x0 - (div_res.quot * y0);
y1 = x1 - (div_res.quot * y1);
x0 = t0; x1 = t1;
// Invariant: (x0*a+x1*b >= y0*a+y1*b)
div_res = std::div(x0*a+x1*b, y0*a+y1*b);
}
}
</pre><br />
Following are few examples obtained by calling <b>ExtendedEuclid(a,b,x,y)</b>:<br />
<pre>a=25375,b=667 x=1, y=-38, gcd=29
29 = 25375 * 1 + 667*38
a=17765,b=13481 x=-343, y=452, gcd=17
17 = 17765 * -343 + 13481 * 452
a=19109,b=18321 x=-23, y=24, gcd=197
197 = 19109 * -23 + 18321 * 24
</pre>Vamsi Kundetihttp://www.blogger.com/profile/06920946350286969217noreply@blogger.com0tag:blogger.com,1999:blog-6902761.post-43771954629834248392017-05-28T11:54:00.000-07:002017-05-28T11:54:31.352-07:00Efficient polygon area computation with less multiplications.<p>Linear algebra computations typically try to prefer mathematical formulations which preserve numerical precision, reduce multiplications and avoid divisions if possible. More recently I came across an elegant formulation of polygon area computation with <em>half</em> the number of multiplications relative to the standard method. I will try to briefly derive it for posterity.<br />
</p><p><em>Standard Formulation</em>: Let $P = [v_0, v_1, v_2, v_3\ldots v_{n-1}, v_{n}=v_0]$ be a polygon with $n$ vertices, $v_i = (x_i,y_i)$. Recall the standard formula for the area of the polygon $A(P) = \sum_{i=0}^{n-1} (x_{i}y_{i+1} - y_{i}x_{i+1})$. Every term $T(i) = (x_{i}y_{i+1} - y_{i}x_{i+1})$ corresponds to the edge $(v_i, v_{i+1})$ of the polygon -- more precisely this is the signed area of the triangle ${(0,0), (x_i,y_i), (x_{i+1}, y_{i+1})}$. Clearly computation of $A(P) = \sum_{i=0}^{n-1} T(i)$ needs exactly $2n$ multiplication operations.<br />
</p><br />
<p><em> Alternate Formulation:</em> Consider $T(i) = x_{i}y_{i+1} - y_{i}x_{i+1}$, now add $x_iy_i$ and subtract $x_{i+1}y_{i+1}$ resulting in the following:<br />
$$<br />
\begin{array}{lll}<br />
D(i) &=& x_iy_i + T(i) - x_{i+1}y_{i+1} \\<br />
D(i) + D(i+1) &=& x_iy_i + T(i) - x_{i+1}y_{i+1} + x_{i+1}y_{i+1} + T(i+1) - x_{i+2}y_{i+2} \\<br />
&=& x_iy_i + T(i) + T(i+1) - x_{i+2}y_{i+2}\\<br />
\sum_{i=0}^{n-1}D(i) &=& x_{0}y_{0} + \sum_{i=0}^{n-1} T(i) - x_ny_n\,\, \mbox{notice that}\, (x_0,y_0) = (x_n,y_n)\\<br />
&=& \sum_{i=0}^{n-1} T(i)<br />
\end{array}<br />
$$<br />
<br />
Whats special about $D(i) = x_iy_i + T(i) - x_{i+1}y_{i+1}$ ?. Notice that $D(i)$ can be written as the following:<br />
$$<br />
\begin{array}{lll}<br />
D(i) &=& x_iy_i + T(i) - x_{i+1}y_{i+1} \\<br />
&=& x_iy_i + (x_iy_{i+1} - x_{i+1}y_i) - x_{i+1}y_{i+1} \\<br />
&=& (x_i + x_{i+1})(y_i - y_{i+1}) \,\, \mbox{notice the only one multiplication}.<br />
\end{array}<br />
$$<br />
<br />
Finally we state that the area of the polygon is equivalent to $A(P) = \sum_{i=0}^{n-1} (x_i+x_{i+1})(y_i - y_{i+1})$ and needs exactly $n$ multiplications instead of $2n$ with standard formulation.<br />
</p>Vamsi Kundetihttp://www.blogger.com/profile/06920946350286969217noreply@blogger.com0tag:blogger.com,1999:blog-6902761.post-37634759020993285222017-01-20T16:46:00.003-08:002017-01-27T15:43:18.732-08:00An Efficient Algorithm to Determine the Order of a Permutation.Consider at set $S = \{1,2,3\ldots n\}$ and recall that $S_n$ a set of all permutations of $S$ form a <em>group</em> under <em>composition</em> (<$S_n, *$>). E.g. $S=\{1,2,3,4\}$ , $\pi_1 = [3\,4\,2\,1]$ and $\pi_2 = [4\,1\,3\,2]$. To understand composition operation it may be easier to look at a permutation $\pi \in S_n$ as <em> bijective function </em> $\pi : S \rightarrow S$. The composition $\pi_1 * \pi_2$ would be equivalent to $\pi_1(\pi_2(i)), \forall i\in S$ which results in $\pi_1 * \pi_2 = [1\,3\,2\,4] $. We use the notation $\pi^k = \pi*\pi*\pi\ldots *\pi$ to denote a composition applied $k$ times. Since <$S_n,*$> is a group it must have an identity element $e$ which is a permutation such that $e(i) = i, \forall i\in S$. <br />
<br />
Now lets focus on the following question:<br />
<b> Given any permutation $\pi \in S_n$ is there an integer $k \geq 0$ such that $\pi^k = e$ ? </b><br />
<br />
<br />
Lets try few examples with $\pi_1 = [3\, 4\, 2\, 1]$, $\pi_1^2 = [2\, 1\, 4\, 3]$, $\pi_1^3 = [4\, 3\, 1\, 2]$, $\pi_1^4 = [1\, 2\, 3\, 4]$. So for this example $\pi_1$ we were able to find $k=4$ such that $\pi_1^4 = e$. Later we will show that finding such a $k$ is always possible for any permutation $\pi$ and such a $k$ is called the <b><em> order </em></b> of the permutation. <br />
<br />
Recall that every permutation can be decomposed into one or more <a href="http://mathworld.wolfram.com/PermutationCycle.html">disjoint cycles</a>. We call a permutation <em>cyclic</em> if it has exactly one cycle in it. E.g. $\pi_c = [3\, 1\, 4\, 2]$ is a cyclic permutation in $S_4$ with cycle being $1 \rightarrow 3 \rightarrow 3 \rightarrow 4 \rightarrow 2 \rightarrow 1$.<br />
<br />
<b>Lemma:</b> If $\sigma \in S_n$ is a <em>cyclic</em> permutation then $\sigma^n = e$<br />
<br />
<b>Proof:</b> Since $\sigma$ is cyclic all elements $S$ form a directed cycle with $n$ nodes. Starting at any node $i$ we need to traverse $n$ edges to get back to $i$. This can be algebraically expressed as $\underbrace{\sigma(\sigma(\ldots \sigma(i))\ldots)}_{n\text{-times}} =i$ <br />
$\Box$. <br />
<br />
<br />
The above lemma states that if a permutation is cyclic on $n$ elements then composing it $n$ times on itself results in an identity permutation. Now we are ready to prove our main theorem:<br />
<br />
<b>Theorem:</b> $\forall \pi \in S_n$ there exists a $k \geq 0$ such that $\pi^k = e$.<br />
<br />
<b>Proof:</b> Let $\pi = \sigma_1\sigma_2\ldots \sigma_m$ be the cycle decomposition of $\pi$. Let $|\sigma_i|$ denote size of cycle. Since each of $\sigma_i$ is cyclic permutation in $S_{|\sigma_i|}$ by the previous lemma $\sigma_i^{|\sigma_i|} = e_{|\sigma_i|}$. Where $e_{|\sigma_i|}$ is an identity permutation in $S_{|\sigma_i|}$. Also notice that any multiple of $|\sigma_i|$ is also going to result in an identity permutation (i.e. for any integer $q \geq 1$ , $\sigma_i^{q|\sigma_i|} = e_{|\sigma_i|}$). Let $t = |\sigma_1|\times |\sigma_2|\times \ldots |\sigma_m|$ then $\pi^t = \sigma_1^t\sigma_2^t\ldots \sigma_m^t$ since $t$ is a multiple of $|\sigma_i|, 1\leq i \leq m$ we have $\pi^t = e_{|\sigma_1|}e_{|\sigma_2|}\ldots e_{|\sigma_m|} = e$ <br />
$\Box$.<br />
<br />
<br />
From the above proof its clear that choosing $k = |\sigma_1|\times |\sigma_2|\ldots |\sigma_m|$ which is a multiple of all cycle lengths guarantees that $\pi^k = e$. However from a computational efficiency purpose we are much better if we let that multiple be the <b> <em> least common multiple</em> </b> of $|\sigma_i|,\, 1\leq i\leq m$. Following is an efficient algorithm -- linear time and takes at most $n$ l.c.m computations -- to find the order of any permutation:<br />
<br />
<pre class="brush:cpp;">// Finds the order of a valid permutation ( 1 \leq perm[i] \leq n).
size_t FindPermutationOrder(const std::vector< unsigned int > &perm) {
size_t n = perm.size();
size_t perm_order = 1; // This overflows if the l.c.m does not fit in 64-bits.
std::vector< bool > cycle_mask(n, true);
// Find the disjoint cycle lengths in the permutation.
size_t cycle_index = 0;
while (cycle_index < n) {
// Find the start of the next disjoint cycle.
if (cycle_mask[cycle_index]) {
size_t i = cycle_index;
size_t cycle_length = 0;
do {
cycle_mask[i] = false;
i = perm[i] - 1;
cycle_length++;
} while (i != cycle_index);
perm_order = boost::math::lcm(perm_order, cycle_length);
}
cycle_index++;
}
return perm_order;
}
</pre>Following are few examples on the algorithm: <br />
<pre class="brush:cpp;">Order of [ 5 3 1 2 4 ] is 5
Order of [ 3 4 2 1 ] is 4
Order of [ 1 2 3 4 5 ] is 1
Order of [ 7 8 9 10 1 2 4 3 6 5 ] is 5
</pre>Vamsi Kundetihttp://www.blogger.com/profile/06920946350286969217noreply@blogger.com0tag:blogger.com,1999:blog-6902761.post-85707634842044704262016-12-22T12:37:00.004-08:002016-12-22T15:19:19.534-08:00Simpler Expected Analysis of Randomized Selection Algorithm.In this post I present a simpler proof to derive expected asymptotic bounds for the <em> Randomized Selection </em> algorithm -- <a href="https://en.wikipedia.org/wiki/Las_Vegas_algorithm"> Las Vegas </a> version to be more precise. Given a set $S = \{ a_1\,a_2\,a_3\ldots a_n \}$ of keys and an integer $k \in [1,n]$ the <em> selection problem </em> asks to report the $k^{th}$ largest key in the set $S$. Recall that this problem can be solved optimally using $O(n)$ comparison operations using the BFPRT (Blum, Floyd, Pratt, Rivest and Tarjan 1973) algorithm. However a correct implementation of BFPRT algorithm may take a lot of code. On the other hand randomized algorithms for many deterministic algorithms are often simpler to implement. The algorithm chooses a random pivot for partitioning the input and recursively picks on of the partition until the partition (sub problem) size is just one. Following is an implementation of <tt><b>RandSelect</b></tt> in around 25 lines. <br />
<br />
Our goal here is the analyze the <em><b> expected run-time </b></em> and <i><b> expected number of steps </b> </i> for the problem size converge to <em> unity </em> -- more precisely the expected value of the variable <tt><b>rand_walk_steps</b></tt> in the following code. I have purposefully kept the analysis very rigorous to make it accessible to everyone. <br />
<br />
TL;DR "Expected number of comparisons in <tt><b>RandSelect</b></tt> is $\Theta(n)$ and expected number of <tt><b>rand_walk_steps</b></tt> is $\Theta(\log(n))$".<br />
<br />
<pre class="brush:cpp;" >// Randomized selection algorithm to pick k^th largest element
// in a mutable array.
template < typename T >
T RandSelect(T *array, size_t n, size_t k) {
assert(k > 0 && k <= n);
size_t problem_size = n, select = k;
// What is the expected value of rand_walk_steps ??
size_t rand_walk_steps = 0;
while (problem_size > 1) {
size_t i = 0, j = problem_size - 1;
size_t ridx = RandArrayIndex(problem_size);
const T pivot = array[ridx];
// Partition the sub problem with pivot.
while (i <= j) {
if (array[i] <= pivot) {
i++;
} else if (array[j] > pivot) {
j--;
} else {
std::swap(array[i], array[j]);
}
}
problem_size = (i >= select) ? i : problem_size - i;
array = (i >= select) ? array : &array[i];
select = (i >= select) ? select : select - i;
++rand_walk_steps;
}
return array[0];
}
</pre><br />
<br />
Let $X_i$ be the random variable corresponding to problem size in $i^{th}$ step ($1 \geq i \geq n-1$) in the algorithm -- each step corresponds to one iteration of the outer most <tt><b>while</b></tt> loop. If the algorithm terminates in $j$ steps then all $X_k$ with $k>j$ are set to $0$.<br />
$$<br />
\begin{array}{lll}<br />
X_i &=& \text{Random variable for problem size in }i^{th}\text{ iteration of the outer most}\textsf{ while}\text{ loop} \\<br />
X_i && \left\{ \begin{array}{lr} \in [1,X_{i-1}-1] & \textsf{If } X_{i-1} > 1 \\<br />
= 0 & \textsf{Otherwise}<br />
\end{array} \right. \\<br />
X_0 &=& n \,\,\,\,\,\,\, \text{Initial problem size and }\,\,\, E[X_0] = n \\<br />
&& \\<br />
\end{array}<br />
$$<br />
<p><b> Lemma-1: </b> If $X_{i-1} >1$ Then $E[X_i | X_{i-1}] = \frac{X_{i-1}}{2}$ <br />
<br />
<em> Proof: </em> For a given $X_{i-1} > 1$ the random variable $X_{i}$ is uniformly distributed in the range $[1,X_{i-1}-1]$. Hence the expected value is simply the average $\frac{1 + (X_{i-1}-1)}{2} = \frac{X_{i-1}}{2}$ $\Box$.<br />
</p><p><b> Lemma-2: </b> $E[X_i] = \frac{n}{2^i}$</br><br />
<em> Proof: </em> <br />
$$<br />
\begin{array}{llll}<br />
&E[X_i|X_{i-1}] &=& \frac{X_{i-1}}{2}\,\,\,\,\,\text{From Lemma-1}\\<br />
&E[\,E[X_i|X_{i-1}]\,] &=& E[\frac{X_{i-1}}{2}] \\<br />
\Rightarrow & E[X_i] &=& E[\frac{X_{i-1}}{2}] = \frac{E[X_{i-1}]}{2}\,\,\,\,\,\text{for any two random variables } X,Y\,\,E[E[X|Y]] = E[X]\\<br />
&&=&\frac{1}{2}\times E[X_{i-1}] = \frac{1}{2^2}\times E[X_{i-2}] \ldots \,\,\,\,\, \text{Apply recursively} \\<br />
&&\ldots& \\<br />
&&=& \frac{1}{2^i}\times E[X_{i-i}] = \frac{E[X_0]}{2^i} = \frac{n}{2^i} \,\, \Box\\<br />
\end{array}<br />
$$<br />
</p><p><b> Theorem-1: </b> Expected number of comparisons done by <tt>RandSelect</tt> is $\Theta(n)$.<br />
<br />
<em> Proof: </em> The number comparisons done by the algorithm in each iteration is $\leq 3\times X_i$. We at most $2$ comparisons in the inner <tt> if..else</tt> branch and $1$ comparison to test $i <= j$, so a total of $3$ comparisons in the each iteration of the inner most <tt>while</tt> loop. The loop runs at most $X_i$ times in the $i^{th}$ iteration of the outer most <tt> while</tt> loop. Finally we need $1$ comparison to test the termination of the <tt>while</tt> loop. Let $T$ be the random variable corresponding the total number of comparisons in the algorithm. So we need to find $E[T]$. <br />
$$<br />
\begin{array}{ll}<br />
&T \leq 3(X_0 + X_1 + X_2 + \ldots X_{n-1}) + 1 & \text{Total number of comparisons in the Algorithm. }\\<br />
&&\\<br />
\Rightarrow& \sum_{i=0}^{n-1} X_i \leq T \leq 3\sum_{i=0}^{n-1} X_{i} + 1 & \\<br />
& \sum_{i=0}^{n-1} E[X_i] \leq E[T] \leq 3\sum_{i=0}^{n-1} E[X_{i}] + 1 & \text{Linearity of expectation.} \\<br />
&&\\<br />
& \sum_{i=0}^{n-1} \frac{n}{2^i} \leq E[T] \leq 3\sum_{i=0}^{n-1} \frac{n}{2^i} + n & \text{Apply Lemma-2.} \\<br />
& 2n\left(1-\frac{1}{2^n}\right) \leq E[T] \leq 6n\left(1-\frac{1}{2^n}\right) +n & \text{Use: } \sum_{i=0}^{n-1}\frac{1}{2^i} = \frac{1-1/2^n}{1-1/2}.\\<br />
&&\\<br />
\Rightarrow &n \leq E[T] < 6n + 1 = 6n +1 & \text{Use: for } n\geq 1\,\,\,,\frac{1}{2} \leq \left(1-\frac{1}{2^n}\right) < 1 \\
&&\\
\Rightarrow & E[T] = \Theta(n)\,\,\Box & \\
\end{array}
$$
</p><br />
<br />
<br />
<br />
<br />
<br />
We now return the question of the expected number of steps for the problem size to reach unity (i.e. expected number of iterations of the outer most <tt>while</tt> loop). This is also referred as the expected number of <em>recursive</em> calls in Motwani and Raghavan's book. However the following analysis is much simpler than the application of a calculus based lemma based on random walks in the book. Let $Y_{i}$ be a random variable corresponding to the number of iterations starting with a problem size of $X_{i} > 1$.<br />
$$<br />
\begin{array}{lll}<br />
Y_{i} &:=&\text{random variable for the number of steps with problem size }\,\, X_{i} >1 \\<br />
Y_i &=& 1 + Y_{i+1} \,\,\,\,\,\, \text{(Recursive formulation)}\\<br />
Y = Y_0 &=& \underbrace{1 + 1 + 1 \dots + 1}_{k \text{-times}} + Y_{k} \,\,\,\,\,\, \text{(Expand } k \text{ times)}\\<br />
&& \\<br />
\end{array}<br />
$$<br />
<p><b> Theorem-2: </b> Expected number of steps for the problem size to reach unity, $E[Y] = \Theta(\log(n))$. <br />
<br />
<em> Proof: </em><br />
$$<br />
\begin{array}{ll}<br />
&Y_i \leq X_{i}\,\,\,\, \text{ (Number of steps are at most problem size at } i^{th} \text{ step which is } X_i)\\<br />
& \\<br />
\Rightarrow &E[Y_i] \leq E[X_i] \,\,\,\, \text{( Let } Z = X_i-Y_i\,\,\,\,,\,\,\,\, E[Z] = E[X_i]-E[Y_i] \geq 0 \text{)}\\<br />
&E[Y_i] \leq E[X_i] = \frac{n}{2^i} \,\,\,\, \text{ (Apply Lemma-2)} \\<br />
&Y = Y_0 = \underbrace{1 + 1 + 1 \dots + 1}_{t \text{-times}} + Y_{t} \\<br />
\Rightarrow & E[Y] = \underbrace{E[1] + E[1] \dots + E[1]}_{t \text{-times}} + E[Y_{t}] = t + E[Y_{t}] \\<br />
&E[Y] = t + E[Y_{t}] \leq t + E[X_{t}] \leq t + \frac{n}{2^t} \\<br />
&\\<br />
\Rightarrow &t \leq E[Y] \leq t + \frac{n}{2^t} \,\,\,\, \text{After } t \text{ steps }\\<br />
& \\<br />
&\text{When } t = \log_2(n) \text{ the above inequality becomes the following } \\<br />
& \log_2(n) \leq E[Y] \leq \log_2(n) + \frac{n}{2^{log_2(n)}} = \log_2(n) +1 \\<br />
<br />
\Rightarrow & E[Y] =\Theta(log(n)) \,\, \Box \\<br />
\end{array}<br />
$$<br />
</p><br />
<br />
<br />
This completes our expected analysis of the <tt><b>RandSelect</b></tt> algorithm. The analysis just uses elementary probability theory without relying on a Calculus based random walk theorem in the <i> <a href="https://www.amazon.com/Randomized-Algorithms-Rajeev-Motwani/dp/0521474655"> "Randomized Algorithms" </a> </i> book by Motwani and Raghavan -- (<b>Theorem-1.3 </b> on page-15).Vamsi Kundetihttp://www.blogger.com/profile/06920946350286969217noreply@blogger.com0tag:blogger.com,1999:blog-6902761.post-2499205964029844782016-11-24T21:54:00.000-08:002016-11-25T12:14:46.363-08:00Integer Division Algorithm with no Division or Multiplication Operators.I was working through a number theory problem today and I needed an efficient algorithm to find a largest integer $n$ s.t $a\times n \leq b$ for any given integers $a, b$ (w.l.o.g $a < b$). Clearly $n$ is the quotient when $b$ is <i>divided</i> by $a$. Following is an algorithm which can compute the <i>quotient</i> and <i>reminder</i> in $\Theta(\log(a) + \log(b))$ operations without using multiplication or division operators. In fact if the number of bits in $a$ and $b$ are known then the algorithm runs in $\Theta(\log(a) - \log(b))$ operations. <br />
<br />
<pre class="brush:cpp;" >typedef unsigned int uint;
// Integer division of b with a without using division
// or multiplication operators.
div_t integer_div(uint a, uint b) {
div_t result;
uint bit_diff = bit_count(a) - bit_count(b);
uint prod = 0, prod_shift = b << bit_diff, shift = 1 << bit_diff;
result.quot = 0;
do {
if (prod + prod_shift <= a) {
prod += prod_shift;
result.quot += shift;
}
shift >>= 1;
prod_shift >>= 1;
} while (shift);
result.rem = a - prod;
return result;
}
// utility function to find number of bits
// required to represent the value of a.
uint bit_count(uint a) {
unsigned int i = 0;
do {
++i;
} while ((a = a >> 1));
return i;
}
</pre><br />
The above algorithm was tested on $4294967293$ pairs of $(a,b), a,b \in [1, INT\_MAX]$ and compared with <i>std::div</i> -- running all these combinations just take around $2$ minutes on a mac book pro. It makes sense because that $\log(a) \leq 32$ (word size on the machine) hence we were just running around four billion constant operations. Typically in all algorithm analysis the word size is considered a $O(1)$, so as long as the inputs bit-widths are within the word size of the machine the runtime of the algorithm can be considered $O(1)$.<br />
<br />
<br />
Notice that the algorithm is <i>greedy</i> in nature and will provide a correctness next. <br />
<br />
Let $c=\sum_{i=0}^{\log(a)-\log(b)}t_i\times 2^{i}$ and $t_{\log(a)-\log(b)}t_{\log(a)-\log(b)-1}\ldots t_0$ the bit representation of $c$. The algorithm tries to pick the powers of $2$ in a greedy manner (i.e. $2^{\log(a)-\log(b)}\geq 2^{\log(a)-\log(b)-1} \geq \ldots 2^{0}$). Consider the choice of $i^{th}$ bit of $c$, if $2^{i}\times a \leq b$ then any $i$-bit integer with a $0$ in the $i^{th}$ bit is strictly less than $2^{i}$. Since $2^{i-1} + 2^{i-2} + \ldots 2^{0} = \frac{2^{i-1+1} -1}{2-1} = 2^{i}-1 < 2^{i}$. We can use induction of $i$ to complete the proof. Hence we can claim that if we greedily pick the powers of $2$ (i.e ${2^i\,\,| \log(a)-\log(b)\leq i \leq 0 }$) it will yield the largest integer $c$ s.t $a\times c \leq b$. The integer division algorithm is just corollary of this basic fact.<br />
<br />
<br />
<br />
<br />
<br />
<br />
<br />
<br />
<br />
Vamsi Kundetihttp://www.blogger.com/profile/06920946350286969217noreply@blogger.com0tag:blogger.com,1999:blog-6902761.post-17591683487554024832016-06-11T21:59:00.002-07:002016-06-12T07:16:01.772-07:00What kind of strings have Kolmogorov complexity equal to their length ?I have been doing some reading about Information Theory and Kolmogorov complexity lately. I just thought about the following problem. Consider a language of all <b> C </b> programs (strings) which print themselves, formally Let $$L_{self} = \{ w\,\, | w\,\, \mbox{is a C program which prints itself} \}$$.<br />
Following is an example of a <b>C</b> program which prints itself (on a side-note the trick behind the construction is to use the ASCII values for <i>quotes</i> and <i>newlines</i>). Clearly the Kolmogorov complexity of each <b>w</b> is less than or equal to the length of <b>w</b> -- given that we have constructed a program of length same as the string <b>w</b>. <br />
I'm wondering if this is indeed optimal ? -- that is given a <b>w</b>, if its possible to produce a program to generate <b>w</b> and size strictly less than <b>|w|</b>.<br />
<br />
<pre class="brush:cpp;" >char qt=34;
char nl=10;
char *str="char qt=34;%cchar nl=10;%cchar *str=%c%s%c;%cint main() {%c printf(str,nl,nl,qt,str,qt,nl,nl,nl,nl);%c}%c";
int main() {
printf(str,nl,nl,qt,str,qt,nl,nl,nl,nl);
}
</pre><br />
BTW, compile the program with <pre>gcc -w </pre>to suppress warnings about header inclusion.Vamsi Kundetihttp://www.blogger.com/profile/06920946350286969217noreply@blogger.com0tag:blogger.com,1999:blog-6902761.post-25030746595894343282015-03-07T15:48:00.001-08:002016-06-11T21:27:58.403-07:00Some Special Cases Of the Sequence Assembly Problem.<p><em> Sequence Assembly </em> (SA) is a well studied problem in Bioinformatics. To quickly summarize the problem: let $\Sigma$ be some alphabet and $S =\{ s_1 , s_2 \ldots s_n\}$ be a set of strings from $\Sigma$. The goal of the SA problem is re-construct a string $s_a$ such that every $s_i \in S$ is a sub-string in $s_a$. Although its very tempting to reduce the SA problem to the <em> Shortest Super String Problem</em> (which happens to be $\mathcal{NP}$-hard), the <em> repeats</em> in genome makes the SA problem complex and ambiguous. <br />
<br />
Any way one of my friends showed me this <a href="https://www.hackerrank.com/challenges/favourite-sequence">problem</a> yesterday, which is a much simpler special case of the SA problem. Let $s_{\pi}$ be a some permutation (unknown) of $\Sigma$. Now given a set $S' = \{s'_1,s'_2\ldots s'_m\}$ of <em> sub-sequences </em> of $s_{\pi}$, we want to re-construct this un-known $s_{\pi}$. Just like the SA problem for a given $S'$ there could be several answers for $s_{\pi}$ -- in which case we want to find out the lexicographically smallest permutation. If $\Sigma_i |s'_i| = N$ the problem could be solved efficiently by reducing it to a <em>topological sorting</em> algorithm on a DAG. Following is a quick solution to this problem I wrote while watching the RSA vs PAK world cup cricket match -- I'm little disappointed that RSA lost despite a ferocious knock by AB. The construction of the DAG from $S'$ takes $O(Nlog(N)$ operations, however the cost of topological sorting to report the answer is still $O(N)$ operations.<br />
<br />
<pre class="brush:cpp;">// 03/06/2015: vamsik (at) engr.uconn.edu//
#include <cmath>
#include <cstdio>
#include <map>
#include <list>
#include <unordered_map>
#include <iostream>
#include <algorithm>
#include <cassert>
using namespace std;
struct CDirectedGraph{
map<unsigned int, map<unsigned int, bool> > m_aEdges;
unordered_map<unsigned int, unsigned int> m_aInDegree;
inline CDirectedGraph() : m_aEdges(), m_aInDegree() {}
inline void add_edge(unsigned int i, unsigned int j){
if((m_aEdges[i]).find(j) == (m_aEdges[i]).end()){
if(m_aInDegree.find(j) == m_aInDegree.end()){
m_aInDegree[j] = 0;
}
m_aInDegree[j]++;
(m_aEdges[i])[j] = true;
}
}
inline void top_sort(){
map<unsigned int, bool> snodes;
unsigned int v;
for(auto vItr=m_aEdges.begin(); vItr!=m_aEdges.end(); ++vItr){
if(!m_aInDegree[vItr->first]){
snodes[vItr->first] = true;
}
}
while(!snodes.empty()){
auto vItr = snodes.begin();
v = vItr->first;
auto eEnd = m_aEdges[v].end();
auto eBeg = m_aEdges[v].begin();
for(auto eItr=eBeg; eItr!=eEnd; ++eItr){
m_aInDegree[eItr->first]--;
if(!m_aInDegree[eItr->first]){
snodes[eItr->first] = true;
}
}
printf("%u ", v);
snodes.erase(vItr);
}
}
};
int main() {
unsigned int N;
unsigned int K, pnode, cnode;
scanf("%u", &N);
CDirectedGraph diGraph;
for(unsigned int i=0; i<N; i++){
scanf("%u", &K);
scanf("%u", &pnode);
for(unsigned int j=0; j<K-1; j++){
scanf("%u", &cnode);
diGraph.add_edge(pnode, cnode);
pnode = cnode;
}
}
diGraph.top_sort();
return 0;
}
</pre>Vamsi Kundetihttp://www.blogger.com/profile/06920946350286969217noreply@blogger.com0tag:blogger.com,1999:blog-6902761.post-21222361780252444722015-02-01T16:17:00.001-08:002015-02-03T16:38:42.388-08:00Some properties of a PinWheel Polygon<div class="separator" style="clear: both; text-align: center;">
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg_FdvCzUtuut8Yu3M8k3KwqlwyQ4XHL_CxQX2LPoQ8vH9Ha9sVO1BXtnOtvSvuJRtQ8YXZZmSSk44DDTNmHKq_DBdxA3O1GaN1Fvt2KxSIRh7-Gca_yfZADYh-Pyl3uM58e9VKOA/s1600/20150201_160819.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEg_FdvCzUtuut8Yu3M8k3KwqlwyQ4XHL_CxQX2LPoQ8vH9Ha9sVO1BXtnOtvSvuJRtQ8YXZZmSSk44DDTNmHKq_DBdxA3O1GaN1Fvt2KxSIRh7-Gca_yfZADYh-Pyl3uM58e9VKOA/s1600/20150201_160819.jpg" height="640" width="480" /></a></div>
<div class="separator" style="clear: both; text-align: center;">
</div>
Vamsi Kundetihttp://www.blogger.com/profile/06920946350286969217noreply@blogger.com0tag:blogger.com,1999:blog-6902761.post-81444455828012929072014-11-16T15:35:00.003-08:002016-06-11T21:31:56.669-07:00Pseudo Polynomial Dynamic Programming<em> Sub Set Sum </em> and <em> Integer Knapsack </em> are few examples of combinatorial optimization problems with <em> pseudo </em> polynomial running time algorithms. Often we the implementation of these pseudo polynomial dynamic programming algorithms is done using a huge array to hold intermediate sub problems (e.g. $OPT[K][n]$ would indicate a presence of a subset in the first $n$ integer elements ($x_1,x_2\ldots x_n$) which sum to $K$). If $N$ is the bound on the subset value you were looking for then the space complexity would $\Theta(Nn)$ -- which is clearly exponential on the input size which is $O(n\log(N))$ bits. To build robust algorithms -- reduce the dependence on the value of $N$ -- for these pseudo polynomial algorithms, it would be efficient if we could exploit the underlying DAG (Directed Acyclic Graph) nature of any dynamic programming formulation. But that might be a little bit more code to traverse the sub-problems in a topological order. However we can make our life a little easy by using a hash table to keep track of the sub-problems, which might add to the worst case asymptotic runtime by a factor of $O(log(nN))$. On the other hand the average case of sub-problem lookup would be a constant and would run fairly efficiently in practice. <br />
<br />
The following code is a quick illustration of using a hash tables to keep track (and also lookup) of sub-problems in pseudo polynomial dynamic programming algorithms. You can read the problem statement <a href="https://www.hackerrank.com/contests/mar14/challenges/the-indian-job"> here </a>. Notice the obvious generalization of the problem (scheduling with constant number of resources) from the way I have formulated the dynamic program.<br />
<br />
<br />
<pre class="brush:cpp;">// 11/15/2014: vamsik@engr.uconn.edu//
#include <cmath>
#include <cstdio>
#include <vector>
#include <iostream>
#include <algorithm>
#include <unordered_map>
using namespace std;
typedef std::pair<unsigned int, unsigned int> KeyType;
struct KeyTypeHasher{
inline size_t operator()(const KeyType& a) const {
return (a.first)*1729 + a.second;
}
};
typedef std::unordered_map<KeyType, bool, KeyTypeHasher> SubProbType;
typedef typename SubProbType::iterator SubProbItr;
bool IsFeasible(const std::vector<unsigned int>& A, unsigned int G){
SubProbType sub_problems_even, sub_problems_odd;
SubProbItr itr;
//init//
if(A[0] <= G){
sub_problems_even[KeyType(A[0], 0)] = true;
sub_problems_even[KeyType(0, A[0])] = true;
}
for(size_t i=1; i < A.size(); i++){
SubProbType& curr_problems =
(i%2) ? sub_problems_even : sub_problems_odd;
SubProbType& next_problems =
(i%2) ? sub_problems_odd : sub_problems_even;
if(curr_problems.empty()){
return false;
}
next_problems.clear();
//create new sub-problems in the next level//
for(itr = curr_problems.begin();
itr != curr_problems.end(); ++itr){
const KeyType &a = itr->first;
if( (a.first + A[i]) <= G){
next_problems[ KeyType((a.first+A[i]), a.second)] = true;
}
if( (a.second + A[i]) <= G){
next_problems[ KeyType(a.first, (a.second + A[i]))] = true;
}
}
curr_problems.clear();
}
return ((A.size())%2) ? !(sub_problems_even.empty()) : !(sub_problems_odd.empty());
}
int main() {
/* Enter your code here. Read input from STDIN. Print output to STDOUT */
size_t T,N;
unsigned int G, a;
scanf("%lu", &T);
for(size_t i=0; i < T; i++){
std::vector<unsigned int> A;
scanf("%lu %u",&N, &G);
for(size_t i=0; i < N; i++){
scanf("%u", &a);
A.push_back(a);
}
printf("%s\n", IsFeasible(A, G) ? "YES" : "NO");
}
return 0;
}
</pre>
<br />
Vamsi Kundetihttp://www.blogger.com/profile/06920946350286969217noreply@blogger.com0tag:blogger.com,1999:blog-6902761.post-75571429895313209842014-08-25T18:58:00.000-07:002014-08-25T18:58:05.953-07:00Planarity PuzzleI'm obsessed with the untangling of the planar graphs. I'm now able to untangle a 30-vertex planar graph in ~240sec see below (yellow one is the final solution):<br />
<br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiWVhHeceqCPRlhhfSya_yShBvGfI196egsjP4eI5a50jHqb7eYx-wBGrpXOesZRDMO0t9IgkXQgRBkB1z37zPvOKzQzfEy1Ka9zs3CTaHV8ZqPozR9EjF7tl3dYIARgZ1FkFF_4Q/s1600/30_vertex_random_input.png" imageanchor="1" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiWVhHeceqCPRlhhfSya_yShBvGfI196egsjP4eI5a50jHqb7eYx-wBGrpXOesZRDMO0t9IgkXQgRBkB1z37zPvOKzQzfEy1Ka9zs3CTaHV8ZqPozR9EjF7tl3dYIARgZ1FkFF_4Q/s400/30_vertex_random_input.png" /></a></div><br />
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhk_7858vAhqt3rs-Db2GeJwHEFXNdo1doLUkyEI9Q6bAxtETI2S20P0T9wpviTFdIRCcfs9mxuO0BtK4lDErokxDUSb40PDsMnP93xzraJVCvROforA4WlPHAxl-Lh2fiFhuks7Q/s1600/30_vertex_solution.png" imageanchor="1" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhk_7858vAhqt3rs-Db2GeJwHEFXNdo1doLUkyEI9Q6bAxtETI2S20P0T9wpviTFdIRCcfs9mxuO0BtK4lDErokxDUSb40PDsMnP93xzraJVCvROforA4WlPHAxl-Lh2fiFhuks7Q/s400/30_vertex_solution.png" /></a></div>Vamsi Kundetihttp://www.blogger.com/profile/06920946350286969217noreply@blogger.com0tag:blogger.com,1999:blog-6902761.post-14144580444806114942014-08-10T14:30:00.001-07:002014-08-10T14:32:06.076-07:00Breaking cube into five tets<p dir="ltr">Picture says it all:</p>
<div class="separator" style="clear: both; text-align: center;"> <a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEilfUdiir-bmTeODrA4VHBODTYV33pLwziV_DgKfvCsf_5_zwtDmDigm2T5WomNw2-estwPzUaWBZvrarPuPK8rX0Ri1ZnyW9x3P_x17YH9x4NkQWmJ4V0OuXOq2naNZ41tOVx-Kg/s1600/20140810_142609.jpg" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"> <img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEilfUdiir-bmTeODrA4VHBODTYV33pLwziV_DgKfvCsf_5_zwtDmDigm2T5WomNw2-estwPzUaWBZvrarPuPK8rX0Ri1ZnyW9x3P_x17YH9x4NkQWmJ4V0OuXOq2naNZ41tOVx-Kg/s640/20140810_142609.jpg"> </a> </div>Vamsi Kundetihttp://www.blogger.com/profile/06920946350286969217noreply@blogger.com0tag:blogger.com,1999:blog-6902761.post-71130654739436620872014-01-05T15:38:00.000-08:002014-01-05T15:56:02.716-08:00Useful properties of the function $$f(x) = \frac{x}{x-k}$$The function <b>f(x) = x/(x-k)</b> has some very useful properties -- especially when we restrict <b>x</b> to Integers. Let me start with a simple question about Maxima and Minima of <b> f(x)</b>. Applying elementary Calculus may not be useful since <b>f(x) ↦ +∞ </b> as <b>x ↦ k+</b> (goes to <b>-∞</b> from the other side and is discontinuous). A quick plot on WolframAlpha reveals the vertical asymptotes [<a href="http://wolfr.am/JSXDdV">http://wolfr.am/JSXDdV</a>]. <br />
<br />
<img src="http://www.wolframalpha.com/share/img?i=d41d8cd98f00b204e9800998ecf8427es17u6n7sdj&f=HBQTQYZYGY4TOM3EMI3WENDEGMYDCNBTGYYDKMTBGE3GGMJUMY4Aaaaa"/><br />
<br />
Now lets move our attention to the case where <b>x</b> is restricted to Integers. First lets consider the function <b>f(x)</b> when <b>x</b> takes on integers strictly greater than <b>k</b> (i.e. <b> x > k</b>). Notice that <b>f(x)</b> is decreasing the maximum in this case occurs when <b>x = k+1</b>. This gives us the following useful fact:<br />
$$<br />
f(n) = \frac{n}{n-k} \leq k+1 \,\,\,\,\,\,\, n \in \{ x\, | \, x > k \vee x \in \cal{I} \}<br />
$$<br />
<br />
The above inequality is quite handy in your combinatorial analysis tool book. Let me give you a simple application of this inequality by proving that <b>nlog(n) = O(log(n!))</b> (also recall that <b> log(n!) = O(nlog(n)) </b> which means <b> nlog(n) = Θ(log(n!))</b>).<br />
<br />
$$<br />
\begin{array}{lcl}<br />
&& \frac{n}{n-k} \leq k+1 \,\,\,\,\,\, k \geq 1\\<br />
&& \\<br />
\rightarrow && n \leq (n-k)\times (k+1) \\<br />
\rightarrow && n < (n-k) \times (k+k) = 2 (n-k)\times k<br />
&& \\<br />
&&\\<br />
\rightarrow && n < 2(n-1)\times 1\\<br />
\rightarrow && n < 2(n-2)\times 2\\<br />
\rightarrow && n < 2(n-3)\times 3\\<br />
&& \ldots \\<br />
\rightarrow && n < 2(n-(n-1))\times (n-1) = 2(1\times (n-1)\\<br />
&&\\<br />
&&\mbox{multiplying all the above n-1 inequalities you get the following} \\<br />
&&\\<br />
\rightarrow&& n^{n-1} < \displaystyle\Pi_{k=1}^{n-1} 2(n-k)\times k = 2^n ((n-1)!)^2 \\<br />
\rightarrow&& n^n < 2^n (n!)^2 \\<br />
\rightarrow&& \log(n^n) < \log(2^n(n!)^2) = 2\log(n!) + n \\<br />
\rightarrow&& n\log(n) < 2\log(n!) + n < 2\log(n!) + \log(n!) = 3\log(n!) \\<br />
\rightarrow&& n\log(n) = O(\log(n!)) \,\, \Box\\<br />
\end{array}<br />
$$Vamsi Kundetihttp://www.blogger.com/profile/06920946350286969217noreply@blogger.com0tag:blogger.com,1999:blog-6902761.post-66456512594674440852013-11-30T17:41:00.001-08:002013-11-30T17:41:38.399-08:00Getting around with integer overflows during Combinatorial counting<p>Many closed form expressions characterizing the number of combinatorial objects involve manipulation of raising and falling factorials (e.g. $$\displaystyle\Pi_{i=1}^{k} (n-i+1)$$). These quantities accumulate very fast we need to use extreme caution during computation. Lets take the example in which we want to compute $${n\choose k} = \frac{n\times (n-1)\ldots (n-k+1)}{k\times k-1\times \ldots 1}$$ which is the ratio of two falling factorials which is always an integer. A naive algorithm which computes numerator and denominator will easily overflow even though the actual value is well within the limits of an integer. If $$n_{max}$$ is the largest value represented by the underlying data-type, the naive algorithm would overflow when: $$n\geq \sqrt[k]{n_{max}}$$. One way to fix this issue is to maintain the ratio implicitly as two integers and knocking off the GCD when ever possible. This reduces the overflow during computation see below for a simple example. Thanks to my friend <a href="http://www.cs.duke.edu/~chittu/">Chittu</a> for discussing this over lunch.<br />
</p><br />
<pre class="java" name="code"> public static long computeNchooseK(int n, int k){
assert n >=1 && k<=n : "Invalid n and k to compute";
long num=1, denom=1, factor;
do{
factor = gcd(n,k);
num *= (n--/factor); denom *= (k--/factor);
factor = gcd(num, denom);
num /= factor; denom /= factor;
}while(k>=1);
assert denom == 1 : "Something went wrong in computing NchooseK";
return num;
}
</pre>Vamsi Kundetihttp://www.blogger.com/profile/06920946350286969217noreply@blogger.com0tag:blogger.com,1999:blog-6902761.post-59001798796887299222013-09-21T16:57:00.000-07:002013-09-21T16:57:57.823-07:00$\sum_{k=1}^{n} H_k = (n+1)(H_{n+1} - 1)$$$<br />
\frac{1}{1-z} = 1 + z + z^2 + z^3 \dots \\ <br />
\int \frac{1}{1-z}dz = \frac{z}{1} + \frac{z^2}{2} + \frac{z^3}{3} \dots = \sum_{k=1}^{\infty} \frac{z}{k} z^k = -\ln(1-z)\\ <br />
\\<br />
H_n = \sum_{k=1}^{n} \frac{1}{k} \hspace{0.5in} \mbox{By definition, also notice they are partial sums}\\<br />
\rightarrow \mbox{O.G.F of $H_n$ is}\hspace{0.02in} \frac{1}{1-z}\times (-\ln(1-z)) = \frac{1}{1-z}\ln(\frac{1}{1-z}) \hspace{0.5in}\mbox{from our previous result on partial sums}\\<br />
\rightarrow \frac{1}{1-z}\ln(\frac{1}{1-z}) = H_1\times z + H_2\times z^2 + H_3\times z^3 \ldots \hspace{0.5in} \mbox{(1)}\\<br />
\rightarrow \frac{d}{dz}(\frac{1}{1-z}\ln(\frac{1}{1-z})) = \frac{1}{(1-z)^2}\ln(\frac{1}{1-z}) + \frac{1}{(1-z)^2} = H_1 + 2H_2\times z + 3H_3\times z^2 \ldots \hspace{0.5in}\mbox{(2)}\\<br />
\mbox{Notice that}\hspace{0.02in} \frac{1}{(1-z)^2} = 1 + 2z + 3z^2 \dots \hspace{0.5in} \mbox{(3)}\\<br />
\rightarrow\mbox{(2) - (3)} = \frac{1}{(1-z)^2}\ln(\frac{1}{(1-z)}) = (H_1-1) + (2H_2-2)\times z + (3H_3 -3) \times z^2 \ldots \\<br />
\rightarrow \frac{1}{(1-z)^2}\ln(\frac{1}{(1-z)}) \hspace{0.02in} \mbox{is the O.G.F of sequence} \hspace{0.02in} (n+1)H_{n+1} - (n+1) = (n+1)(H_{n+1}-1) \hspace{0.5in} \mbox{(4)}\\<br />
\\<br />
\mbox{Applying the result of the partial sums to (1) the O.G.F of sequence}\hspace{0.02in} \sum_{k=1}^{n} H_n \hspace{0.02in} \mbox{is} \hspace{0.02in} \frac{1}{(1-z)}\times {\frac{1}{(1-z)}\ln(\frac{1}{1-z})} = \frac{1}{(1-z)^2}\ln(\frac{1}{(1-z)}) \hspace{0.5in} \mbox{(5)}\\<br />
<br />
\mbox{From (4) and (5) the generating functions of both the sequences are same, hence the corresponding terms must be the same}<br />
<br />
<br />
<br />
<br />
<br />
$$Vamsi Kundetihttp://www.blogger.com/profile/06920946350286969217noreply@blogger.com0tag:blogger.com,1999:blog-6902761.post-7369933460556560582013-09-21T15:41:00.000-07:002013-09-21T15:41:57.332-07:00Simple proof for the O.G.F of partial sumsLet $A(z)$ be the O.G.F of the sequence $S = \{a_0, a_1, a_2 \ldots \}$, then we can derive the O.G.F of the sequence $b_n = \sum_{i=0}^{n} a_i$ (i.e. $B = \{a_0, a_0+a_1, a_0+a_1+a_2 \ldots \}$) as follows.<br />
<br />
Consider the sequence $S^1 = \{0 , a_0, a_1, a_2, a_3 \ldots \}$ then O.G.F corresponding to $S^1$ is $0\times z^0 + a_0 \times z + a_1 \times z^2 + a_2 \times z^3 \ldots = z\times A(z)$. Similarly the O.G.F corresponding to the sequence $S^2 = \{0, 0, a_0, a_1, a_3 \ldots \}$ is $z^2 \times A(z)$. Now notice that if we add corresponding terms in the sequences $S^1 , S^2, S^3 \ldots$ then this<br />
results in the sequence $B$. Hence the generating function corresponding to $B$ is $A(z) + z\times A(z) + z^2\times A(z) + z^3\times A(z) \ldots = A(z)\times(1 + z + z^2 + z^3 \ldots) = \frac{A(z)}{1-z}$ $\Box$.Vamsi Kundetihttp://www.blogger.com/profile/06920946350286969217noreply@blogger.com0tag:blogger.com,1999:blog-6902761.post-81786633366914225852013-06-09T21:55:00.001-07:002013-06-11T10:47:07.897-07:00An iterative method to compute the integer linear combination of Greatest Common DivisorEuclid's iterative method to compute the GCD is very well known. You might have also heard that $GCD(a,b)$ can be expressed as integer linear combination of $a,b$ (i.e. $GCD(a,b) = \alpha\times a + \beta \times b \,\, \alpha,\beta \in \mathcal{I}$). The following is an iterative method -- I wrote on a flight between AUS -> DFW just for fun -- which computes this integer linear combination. Runs in $O(\log(\max\{a,b\}))$ time.<br />
<br />
<br />
<pre class="cpp" name="code">template<typename Unit=int>
static inline Unit gcd_integer_linear_combination(const Unit& a, const Unit& b, Unit& alpha, Unit& beta ){
Unit divd = a, divs = b, remd=0, quot;
Unit a_divd[2] = {1,0}; //linear combination of the divident//
Unit a_divs[2] = {0,1}, temp[2] = {0,0}; //linear combination of the divisor//
bool b_flip = false;
if(divd < divs){
divd = b; divs = a;
b_flip = true;
}
//Invariant's: divd = a_divd[0]*a + a_divd[1]*b //
// divs = a_divs[0]*a + a_divs[1]*b //
while(divs){
remd = divd%divs; quot = (divd-remd)/divs;
divd = divs; divs = remd;
temp[0] = a_divs[0]; temp[1] = a_divs[1];
a_divs[0] = a_divd[0] - (quot * a_divs[0]);
a_divs[1] = a_divd[1] - (quot * a_divs[1]);
a_divd[0] = temp[0]; a_divd[1] = temp[1];
}
if(!b_flip){
alpha = a_divd[0]; beta = a_divd[1];
}else{
alpha = a_divd[1]; beta = a_divd[0];
}
return divd;
}
</pre><br />
The following are some examples $GCD(a,b) = (\alpha)*a + (\beta)*b$.<br />
<pre class="cpp" name="code">PASSED: 2 = (1613)*42 + (-8)*8468
PASSED: 1 = (-1684)*6335 + (1641)*6501
PASSED: 5 = (344)*9170 + (-551)*5725
PASSED: 1 = (-1949)*1479 + (308)*9359
PASSED: 1 = (967)*6963 + (-1508)*4465
PASSED: 2 = (858)*5706 + (-601)*8146
PASSED: 6 = (181)*3282 + (-87)*6828
PASSED: 2 = (121)*9962 + (-2450)*492
PASSED: 1 = (-596)*2996 + (919)*1943
PASSED: 1 = (-2348)*4828 + (2085)*5437
PASSED: 1 = (283)*2392 + (-147)*4605
PASSED: 1 = (-61)*3903 + (1546)*154
PASSED: 1 = (122)*293 + (-15)*2383
PASSED: 1 = (-1528)*7422 + (1301)*8717
PASSED: 1 = (615)*9719 + (-604)*9896
PASSED: 1 = (-304)*5448 + (959)*1727
PASSED: 1 = (-139)*4772 + (431)*1539
PASSED: 1 = (493)*1870 + (-93)*9913
PASSED: 4 = (628)*5668 + (-565)*6300
PASSED: 1 = (3461)*7036 + (-2461)*9895
PASSED: 4 = (-60)*8704 + (137)*3812
PASSED: 1 = (77)*1323 + (-305)*334
PASSED: 3 = (662)*7674 + (-1089)*4665
PASSED: 2 = (3)*5142 + (-2)*7712
PASSED: 1 = (1840)*8254 + (-2211)*6869
PASSED: 1 = (2027)*5548 + (-1471)*7645
PASSED: 1 = (929)*2663 + (-897)*2758
PASSED: 2 = (-301)*38 + (4)*2860
PASSED: 2 = (-823)*8724 + (737)*9742
PASSED: 1 = (-3)*7530 + (29)*779
PASSED: 1 = (-1499)*2317 + (1144)*3036
PASSED: 1 = (662)*2191 + (-787)*1843
PASSED: 1 = (10)*289 + (-27)*107
PASSED: 1 = (4289)*9041 + (-4336)*8943
</pre>Vamsi Kundetihttp://www.blogger.com/profile/06920946350286969217noreply@blogger.com0tag:blogger.com,1999:blog-6902761.post-16753538982760867002012-10-05T16:54:00.000-07:002012-10-05T17:06:39.311-07:00Friday fun with a dynamic program.Its been almost an year I wrote my personal blog. I missed so many fun things I used to do as student. Well I'm a pursuit of relentless improvement. Today I just wanted to write a dynamic program -- which I have missed for a while -- so I randomly pick one of the hardest DIV-1 problems and got the dynamic program working. This is a very interesting dynamic program -- a two level one. You can read the problem -- because I cannot' post here-- from this link <a href="http://community.topcoder.com/stat?c=problem_statement&pm=1996&rd=4710">http://community.topcoder.com/stat?c=problem_statement&pm=1996&rd=4710</a> . Following is my code which passed all the 88 system tests.<br />
<br />
<pre class="cpp" name="code">/*
* Solution to problem: http://community.topcoder.com/stat?c=problem_statement&pm=1996&rd=4710 (SRM-178, DIV-1, hard)
*
* Vamsi K. Kundeti 10/05/2012
*
**/
#include<cstdio>
#include<string>
#include<cstdio>
#include<vector>
#include<limits>
using namespace std;
#define MAX_ROWS 50
#define MAX_STROKES 3000
class MinPaint {
private:
/*Optimal number of mismatches to color a row 'i' till 'j' with a
*stoke ending black using 'k' strokes : 'm_aB[j%2][i][k]'
*/
size_t m_aB[2][MAX_ROWS][MAX_STROKES];
size_t m_aW[2][MAX_ROWS][MAX_STROKES];
/*Optimal number of mismatches in the image upto row 'i' using 'k'
*strokes*/
int m_aOPT[2][3001];
//optimal value is stored in m_aB[0][i][j]
vector<string> m_aImg;
public:
MinPaint(): m_aB(), m_aW(), m_aImg() {}
void initRow(size_t idx, size_t maxStrokes){
for(size_t j=0; j<=maxStrokes; j++){
m_aB[0][idx][j] = !j ? 1 : m_aImg[idx][0] == 'W';
m_aW[0][idx][j] = !j ? 1 : m_aImg[idx][0] == 'B';
}
}
size_t findMin(size_t a, size_t b){
return a < b ? a : b;
}
void dynamicProgramPerRow(size_t idx, size_t maxStrokes){
initRow(idx, maxStrokes);
for(size_t i=1; i < m_aImg[idx].size(); i++){
for(size_t j=0; j <= maxStrokes; j++){
if(m_aImg[idx][i] == 'W'){
m_aW[i%2][idx][j] = !j ? i+1 :
findMin(m_aW[(i-1)%2][idx][j], m_aB[(i-1)%2][idx][j-1]);
m_aB[i%2][idx][j] = !j ? i+1 :
findMin(m_aB[(i-1)%2][idx][j]+1, m_aW[(i-1)%2][idx][j-1]+1);
}else{
m_aB[i%2][idx][j] = !j ? i+1 :
findMin(m_aB[(i-1)%2][idx][j], m_aW[(i-1)%2][idx][j-1]);
m_aW[i%2][idx][j] = !j ? i+1 :
findMin(m_aW[(i-1)%2][idx][j]+1, m_aB[(i-1)%2][idx][j-1]+1);
}
}
}
//choose the min between m_aB and m_aW
size_t nidx = m_aImg[idx].size()-1;
for(size_t i=0; i<=maxStrokes; i++){
m_aB[nidx%2][idx][i] = findMin(m_aB[nidx%2][idx][i],
m_aW[nidx%2][idx][i]);
}
}
int leastBad(const vector<string> &Img, const int &maxStrokes){
m_aImg = Img;
for(size_t i=0; i< Img.size(); i++){
dynamicProgramPerRow(i, (size_t)maxStrokes);
}
//now find the optimal across rows//
for(int i=0; i<=maxStrokes; i++){
m_aOPT[0][i] = m_aB[(Img[0].size()-1)%2][0][i];
}
for(size_t j=1; j<Img.size(); j++){
/*choose number of strokes from [1, maxStrokes-j] to paint row 'j'*/
size_t rsize_idx = m_aImg[j].size()-1;
for(int usedStrokes=0; usedStrokes<=maxStrokes; usedStrokes++){
int min = numeric_limits<int>::max();
int current_min;
int min_k=0;
for(int k=0; k<=usedStrokes; k++){
//use 'k' strokes to color current row 'j' and 'usedStrokes-j' for the
//array below 'j'.
current_min =
m_aOPT[(j-1)%2][usedStrokes-k] + (m_aB[(rsize_idx)%2][j][k]);
if(current_min < min){
min = current_min;
min_k = k;
}
}
m_aOPT[j%2][usedStrokes] = min;
}
}
return m_aOPT[(Img.size()-1)%2][maxStrokes];
}
size_t OptPerRow(size_t idx, size_t maxStrokes){
return m_aB[(m_aImg[idx].size()-1)%2][idx][maxStrokes-1];
}
};
/*INPUT FILE FORMAT:
*
*
* X Y
* string1
* string2
* ...
* ...
* ...
* stringX
*
*X is the number of rows
*Y is the maximum number of strokes
**/
int main(int argc, char **argv){
char buf[4096];
vector<string> input;
size_t imgSize, maxStrokes;
MinPaint X;
while(scanf("%lu%lu", &imgSize, &maxStrokes)==2){
input.clear();
while(imgSize-- && scanf("%s", &buf) == 1){
input.push_back(string(buf));
}
int opt = X.leastBad(input, maxStrokes);
printf("%d\n", opt);
}
}
</pre><br />
Have fun -- because you can afford it only on a Friday.<br />
<br />
Vamsi.Vamsi Kundetihttp://www.blogger.com/profile/06920946350286969217noreply@blogger.com0tag:blogger.com,1999:blog-6902761.post-73739760590645893862011-05-14T11:06:00.000-07:002011-05-15T23:20:25.269-07:00[TECH] Enumerating set partitions with $DLX$<p>
I have been a great fan of Dr. Knuth's paper about <a href="http://arxiv.org/abs/cs/0011047"> Dancing Links (Algorithm X) </a>. In that paper
Don suggests how the combinatorial search using backtracking can be improved in practice with the original idea from <a href="http://portal.acm.org/citation.cfm?id=1710829">Hitotumatu and Noshita</a>. One strange thing about Dancing Links is that its not Don's original idea, however Don made it popular by writing a paper about it. Don used the $DLX$ idea to solve the exact cover (matrix cover) problem and he gave several examples (e.g. Covering an rectilinear area with Pentominos) where he certain combinatorial problems to exact cover.
</p>
<p>
I have been waiting to find a problem where I get to apply $DLX$ without reducing it to exact cover problem. I'm really happy now that simultaneous hamming neighborhood problem (given a set $S=\{s_1,s_2\ldots s_n\},\, s_i\in \Sigma^{l}$ of strings compute the set $N_d(S) = \{ t | hamming(s_i,t)\leq d ,\forall s_i \in S\}$) which I have been working lately can be solved by applying $DLX$. As I mentioned in my earlier post the algorithm to compute $N_d(S)$ needs solve a set partition problem -- which further needs an algorithm to enumerate all the subsets of a given size $r$. In my previous post we have seen how to compute $n \choose r$ efficiently with backtracking. I'm now going use $DLX$ to solve the same problem. Notice that $DLX$ is an overkill to this simple problem, however it really makes sense if we look it from the perspective of the partitioning problem -- which is what we want to solve. $DLX$ just plays the role of UNDO operation see below.
</p>
<pre name="code" class="cpp">
/*Implementation of $n \choose r$ using $DLX$*/
enum_count_t EnumSubSetsCountDLX(size_t n, size_t r){
enum_count_t n_choose_r=0;
DLXState *dlx_sentinal = GetEnumSubSetsDLXState(n); /*Initialize the DLXState*/
DLXState *sidx_dat[r+1], **state_idx=sidx_dat; /*item picked in the state*/
/*total # of items picked in a given state*/
size_t scidx_dat[r+1], *state_count_idx=scidx_dat;
size_t level;
state_count_idx--; state_idx--; /*1-indexing*/
level=1; state_count_idx[level]=0; state_idx[level] = dlx_sentinal;
while(level){
if(level == r+1){
n_choose_r++;
level--;
if(level) Dance(state_idx[level]);
}else if(state_count_idx[level]+r-level > n-1){ /*state_idx[level]*/
level--;
if(level) Dance(state_idx[level]);
}else{
state_idx[level] = state_idx[level]->right;
Break(state_idx[level]);
state_count_idx[level]++;
/*track forward*/
state_count_idx[level+1] = state_count_idx[level];
state_idx[level+1] = state_idx[level];
level++;
}
}
CleanDLXState(dlx_sentinal); /*cleanup*/
return n_choose_r;
}
/*Break and Dance operations*/
void Break(DLXState *state){
state->left->right = state->right;
state->right->left = state->left;
}
void Dance(DLXState *state){
state->left->right = state;
state->right->left = state;
}
/*Initializes a new DLX state required to enumerate $n \choose r$*/
EnumSubSetState* GetNewEnumSStateWithDLX(DLXState *dlx_sentinal,
size_t n, size_t r){
EnumSubSetState *estate = malloc(sizeof(EnumSubSetState));
assert(estate);
estate->dlx_sentinal = dlx_sentinal;
estate->state_count_idx = malloc(sizeof(size_t)*(r+1));
estate->state_idx = malloc(sizeof(DLXState *)*(r+1));
assert(estate->state_count_idx && estate->state_idx);
(estate->state_count_idx)--; (estate->state_idx)--;
estate->level=1; estate->state_count_idx[1]=0;
estate->state_idx[1] = estate->dlx_sentinal;
estate->n = n; estate->r = r;
return estate;
}
</pre>
<p> Next I'll post the most interesting part which is to enumerate the partitions using $DLX$.
</p>Vamsi Kundetihttp://www.blogger.com/profile/06920946350286969217noreply@blogger.com1tag:blogger.com,1999:blog-6902761.post-67633059019552083232011-05-06T18:52:00.000-07:002011-05-07T06:36:58.673-07:00[TECH] Enumerating subsets with backtracking and some interesting observations on the backtracking enumeration tree.<p>
During my current research on enumerating all the common <em>Hamming</em> neighborhood of a set of strings. I needed an algorithm which can enumerate seven disjoint subsets from a given set. The fundamental building block in that is an algorithm which enumerates all the subsets of size $r$ from a given set of size $n$. I'm very much aware of the standard recursive algorithm which uses the fact
${n \choose r} = {n-1 \choose r} + {n-1 \choose r-1}$. However I wanted a much efficient version because I need to use this several times and cannot afford the constants in stack based implementations. So I wrote the following to backtracking based algorithm to enumerate the subsets. Once I did that I started observing some very interesting patterns in the number of leaves for each internal backtrack node at level $r-1$.
</p>
<pre><span style=' color: Green;'>/*Back tracking*/</span>
enum_count_t EnumSubSetsPrint(<span style=' color: Blue;'>char</span> *set, size_t n, size_t r){
size_t state_idx_data[r+<span style=' color: Maroon;'>1</span>], i, level;
size_t *state_idx = state_idx_data;
enum_count_t n_choose_r=<span style=' color: Maroon;'>0</span>;
set--;
state_idx--;
<span style=' color: Blue;'>for</span>(i=<span style=' color: Maroon;'>1</span>; i<=r+<span style=' color: Maroon;'>1</span>; i++){ <span style=' color: Green;'>/*initialize*/</span>
state_idx[i] = <span style=' color: Maroon;'>0</span>;
}
level = <span style=' color: Maroon;'>1</span>;
<span style=' color: Blue;'>while</span>(level){
<span style=' color: Blue;'>if</span>(level == r+<span style=' color: Maroon;'>1</span>){
<span style=' color: Green;'>/*back track*/</span>
<span style=' color: Blue;'>for</span>(i=<span style=' color: Maroon;'>1</span>; i<=r; i++){
printf(<span style=' color: Maroon;'>"%c "</span>, set[state_idx[i]]);
}
printf(<span style=' color: Maroon;'>"\n"</span>);
n_choose_r++;
level--;
}<span style=' color: Blue;'>else</span> <span style=' color: Blue;'>if</span>(state_idx[level]+r-level > n-<span style=' color: Maroon;'>1</span>){
level--;
}<span style=' color: Blue;'>else</span>{
<span style=' color: Green;'>/*track forward*/</span>
state_idx[level]++; <span style=' color: Green;'>/*pick a item*/</span>
state_idx[level+<span style=' color: Maroon;'>1</span>] = state_idx[level];
level++;
}
}
<span style=' color: Blue;'>return</span> n_choose_r;
}
</pre>
<p>
By looking at these patterns (at the end of this paragraph) I came up with the following interesting identities for $n \choose 3$ and $n \choose 4$. <br>
<bf>IDENTITY-1:</bf> ${n\choose 3} = \sum_{i=1}^{n-2} (n-i)(i-1)$ <br>
<bf>IDENTITY-2:</bf> ${n\choose 4} = \sum_{j=1}^{n-3}\sum_{i=1}^{n-2-j} (i)(j) $<br>
I'm actually printing the number of leaves under an internal node at level $r-1$ in the backtracking tree. I'm observing the following pattern.
</p>
<pre>
=================
8 3 0
6 5 4 3 2 1
5 4 3 2 1
4 3 2 1
3 2 1
2 1
1
(8 \choose 3) = 56
Please enter n and r < n and print_option
(e.g. 15 14 1 or 15 14 0)
9 3 0
7 6 5 4 3 2 1
6 5 4 3 2 1
5 4 3 2 1
4 3 2 1
3 2 1
2 1
1
(9 \choose 3) = 84
Please enter n and r < n and print_option
(e.g. 15 14 1 or 15 14 0)
10 3 0
8 7 6 5 4 3 2 1
7 6 5 4 3 2 1
6 5 4 3 2 1
5 4 3 2 1
4 3 2 1
3 2 1
2 1
1
(10 \choose 3) = 120
Please enter n and r < n and print_option
(e.g. 15 14 1 or 15 14 0)
11 3 0
9 8 7 6 5 4 3 2 1
8 7 6 5 4 3 2 1
7 6 5 4 3 2 1
6 5 4 3 2 1
5 4 3 2 1
4 3 2 1
3 2 1
2 1
1
===============================
The following are for $n \choose 4$ for various values of $n$.
===============================
8 4 0
5 4 3 2 1
4 3 2 1
3 2 1
2 1
1
4 3 2 1
3 2 1
2 1
1
3 2 1
2 1
1
2 1
1
1
(8 \choose 4) = 70
Please enter n and r < n and print_option
(e.g. 15 14 1 or 15 14 0)
9 4 0
6 5 4 3 2 1
5 4 3 2 1
4 3 2 1
3 2 1
2 1
1
5 4 3 2 1
4 3 2 1
3 2 1
2 1
1
4 3 2 1
3 2 1
2 1
1
3 2 1
2 1
1
2 1
1
1
===============================
</pre>
<p>
Play around with this program <a href="http://trinity.engr.uconn.edu/Enumerate.c">Enumerate.c</a> , <a href="http://trinity.engr.uconn.edu/Enumerate.h">Enumerate.h</a>.
</p>
<p>
Finally I enabled $LaTeX$ setting for the blog.
</p>Vamsi Kundetihttp://www.blogger.com/profile/06920946350286969217noreply@blogger.com0tag:blogger.com,1999:blog-6902761.post-36586279803725730882011-05-05T16:46:00.000-07:002011-05-05T17:10:41.918-07:00[THEORY] Is nature a giant guessing algorithm ?<p> I was reading <a href="http://rjlipton.wordpress.com/">Dr. Lipton's blog</a> on guessing and was really impressed on how he relates the power of guessing with the question whether <b>P=NP</b> or <b>P!=NP</b>. </p>
<p>
Looks like <i>nature</i> on the other had seems to be solving several combinatorial problems almost instantly (e.g. folding of the proteins to minimize its energy), does nature guess the solution to the problem ?. One way to guess a solution to a problem is to look inside the building blocks which formulated the problem. Dr. Lipton gives an example were we were supposed to guess the digits which can open a lock -- the lock we normally use in the GYM. It may be hard to guess the digits which can open a lock without any help. However if we can get an X-ray of the levers we can guess the solution to the problem more easily. So does this mean nature looks at the problems in a totally different perspective ? , is there something very intuitive to the nature which the humans fail to recognize ?.
</p>
<p>
Why did all the great mathematicians failed to find a polynomial solution 3-SAT ? -- or prove it does not exist. Is there something missing in the way we are looking at a digital circuit ? , can we design guessing algorithms to solve the NP-Complete problems ?. We don't know the answers for all these. But I'm sure that nature is a great guessing algorithm and it may be impossible for us to guess its guessing algorithm.
</p>Vamsi Kundetihttp://www.blogger.com/profile/06920946350286969217noreply@blogger.com0tag:blogger.com,1999:blog-6902761.post-5034949927580264042011-04-12T12:46:00.000-07:002011-04-12T13:14:46.486-07:00[TECH] Unbuffered message logging of interactive programs.<p>
We generally use one of the following methods to log the output messages of programs which print a lot of messages.
<ul>
<li> <pre> $ ./program >& message.log </pre>
<li> <pre> $ ./program |& tee message.log </pre>
</ul>
Since <b>stdout</b> is line buffered we expect the same to hold when we re-direct <b>stdout</b> or any other stream. However this is not true if you keep monitoring the program status with <b> tail -f message.log</b> you will not notice any output from the program this is because the redirection is not line buffered. And in certain unfortunate cases when the running program aborts the redirection may not even flush the buffers and you find that <b>message.log</b> is empty and there is NO way to figure out what has happened.
</p>
<p>One obvious way to fix this is to put <b>fflush(stdout)</b> after every <b>printf</b> statement you have, this is OK for small programs but for very large programs its really tedious. I came up with the following technique to solve this problem. This solution is generic.
</p>
<pre>
[main.c]
int main(int argc, char **argv){
ForceLogMessages("message.log");
....
...
}
void ForceLogMessages(const char *mesg_file){
int ret_setvbuf;
/*close the stdout and re-open a text file*/
fclose(stdout);
stdout = fopen(name, "w");
//make stdout unbuffered
ret_setvbuf = setvbuf(stdout, (char *) NULL, _IOLBF, 0);
if(ret_setvbuf){
fprintf(stderr, "UNABLE TO MAKE stdout UNBUFFERED %s", strerror(errno));
}
assert(stdout);
}
</pre>
<p>
Now if we run the program and we can monitor the status with <b>tail -f message.log</b> since we have made <b>stdout</b> unbuffered we can see a message in message.log every time a <b>printf</b> is executed in the original program and you can know what exactly is happening.
</p>Vamsi Kundetihttp://www.blogger.com/profile/06920946350286969217noreply@blogger.com0tag:blogger.com,1999:blog-6902761.post-27292980137911679952011-03-24T18:59:00.000-07:002011-03-24T19:11:35.555-07:00[TECH] juggernaut-asm: An out-of-core sequence assembler<p>
Development updates on the juggernaut-asm project.I'm now branching off the trunk to implement the feature which avoids using the UNIX sort.
To checkout the main branch as follows. This is the version of the code I used to obtain the results on 4million reads which the in-memory algorithm of Velvet could not handle.
</p>
<pre>
cvs -d:pserver:anonymous@juggernaut-asm.cvs.sourceforge.net:/cvsroot/juggernaut-asm co .
</pre>
<p>
The code to replace the unix sort being developed in the branch name 'replace-unix-sort'.
</p>
<pre>
cvs -d:ext:vamsi99@juggernaut-asm.cvs.sourceforge.net:/cvsroot/juggernaut-asm co -r replace-unix-sort .
</pre>Vamsi Kundetihttp://www.blogger.com/profile/06920946350286969217noreply@blogger.com0tag:blogger.com,1999:blog-6902761.post-44690394901256471242011-03-23T13:42:00.000-07:002011-03-23T14:18:48.260-07:00[TECH] External R-Way Merge sorting.<p>
http://lib-ex-sort.sourceforge.net/ is an external sorting program the following are the limitations I wish to remove before I leave in May. Most of them are engineering changes but
are very useful for several algorithms based on external sorting.
</p>
<ul>
<li> Currently the size of each key is a constant. We need to remove this by adding a key_header for each key.
<li> Use MMAP for integer sorting to avoid copy between user and kernel space.
<li> Support for sorting when the data spans multiple files.
</ul>Vamsi Kundetihttp://www.blogger.com/profile/06920946350286969217noreply@blogger.com0tag:blogger.com,1999:blog-6902761.post-8909297520162574902010-08-27T16:19:00.000-07:002010-08-27T16:24:56.574-07:00[TECH] Combinatorial counting under symmetries with Polya's Theorem<p> Recently I have started looking at the problem of enumerating Graphs. During my journey I was thrilled at the beauty of Polya's theorem. After applying the theorem on several examples I felt that I should write an article which explains the context and need for Polya's theorem.
</p>
<p>This is the first blog I'm writing after taking a break over the summer. I had few things which had happened to me recently. I'm going to talk about them some other time as they are more philosophical.
</p>Vamsi Kundetihttp://www.blogger.com/profile/06920946350286969217noreply@blogger.com0tag:blogger.com,1999:blog-6902761.post-89749020858008571792010-04-09T01:49:00.000-07:002010-04-09T02:00:53.923-07:00[TECH] Super-fast tokenizing of sequence word patterns<html xml:lang="en" xmlns="http://www.w3.org/1999/xhtml">
<head>
<meta http-equiv="content-type" content="text/html; charset=utf-8">
<title>C++ code colored by C++2HTML</title>
<meta name="generator" content="C++2HTML by Jasper Bedaux">
<!-- To generate your own colored code visit http://www.bedaux.net/cpp2html/ -->
<style type="text/css">
.comment { color: #999999; font-style: italic; }
.pre { color: #000099; }
.string { color: #009900; }
.char { color: #009900; }
.float { color: #996600; }
.int { color: #999900; }
.bool { color: #000000; font-weight: bold; }
.type { color: #FF6633; }
.flow { color: #FF0000; }
.keyword { color: #990000; }
.operator { color: #663300; font-weight: bold; }
.operator { color: #663300; font-weight: bold; }
</style>
</head><body>
<pre><span class="pre">#!/bin/perl
#
# AHO-CORASICK-TOKENIZER on words, will also work on any delemiter.
# INPUT(1): Set of patterns (sequence of words) 'P'
# INPUT(2): Text 'T'
# OUTPUT: report all the occuring patterns -- we don't even
# want to know their locations.
#
# NOTE: We don't want to worry much about the time required to
# construct the failure function. All we need is to improve the
# scan speed.
#
# PATTERN FILE FORMAT:
# <annotation> $$$ pattern
#
# vamsik@engr.uconn.edu 04/08/2010
#
</span><span class="operator">
(</span>$#ARGV<span class="operator"> ==</span><span class="int"> 1</span><span class="operator">)</span><span class="operator"> or</span> die<span class="operator">(</span><span class="string">"USAGE: perl <name> pattern_file text_file"</span><span class="operator">);</span>
open<span class="operator">(</span>PATTERN_FILE<span class="operator">,</span> $ARGV<span class="operator">[</span><span class="int">0</span><span class="operator">])</span><span class="operator"> or</span> die<span class="operator">(</span>$<span class="operator">!);</span>
open<span class="operator">(</span>TEXT_FILE<span class="operator">,</span> $ARGV<span class="operator">[</span><span class="int">1</span><span class="operator">])</span><span class="operator"> or</span> die<span class="operator">(</span>$<span class="operator">!);
%</span>AC_TREE<span class="operator">;
%</span>PATTERN_HASH<span class="operator">;</span>
$ANNOT_DELIM<span class="operator"> =</span><span class="string"> ";;a;;"</span><span class="operator">;</span>
$FAILED_DELIM<span class="operator"> =</span><span class="string"> ";;f;;"</span><span class="operator">;</span>
$HEIGHT_DELIM<span class="operator"> =</span><span class="string"> ";;h;;"</span><span class="operator">;</span> #to move within the text
sub BuildKwordTree<span class="operator">{</span><span class="pre">
# root of the Aho-Corasick tree #
</span> my @params<span class="operator"> =</span> @_<span class="operator">;</span>
my $ac_root<span class="operator"> =</span> $params<span class="operator">[</span><span class="int">0</span><span class="operator">];</span>
my $PFILE<span class="operator"> =</span> $params<span class="operator">[</span><span class="int">1</span><span class="operator">];</span>
my<span class="operator"> (</span>$line<span class="operator">,</span> $pword_list<span class="operator">,</span> $pattern<span class="operator">,</span> $curr_node<span class="operator">,</span> $pword<span class="operator">);</span><span class="flow">
while</span><span class="operator">(<</span>$PFILE<span class="operator">>){</span>
$line<span class="operator"> =</span> $_<span class="operator">;</span> chomp<span class="operator">(</span>$line<span class="operator">);</span><span class="flow">
if</span><span class="operator">(</span>$line<span class="operator"> =~ /([^\</span>$<span class="operator">]+)\</span>$<span class="operator">\</span>$<span class="operator">\</span>$<span class="operator"> (.*)/){</span>
$annotation<span class="operator"> =</span> $<span class="int">1</span><span class="operator">;</span>
$pattern<span class="operator"> =</span> $<span class="int">2</span><span class="operator">;</span><span class="pre">
# split the pattern into words delimited by any #
# non-newline character. #
</span> @pword_list<span class="operator"> =</span> split<span class="operator">(/\</span>s<span class="operator">+/,</span> $pattern<span class="operator">);</span>
$curr_node<span class="operator"> =</span> $ac_root<span class="operator">;</span>
foreach $pword<span class="operator"> (</span>@pword_list<span class="operator">){</span>
$<span class="operator">{</span>$curr_node<span class="operator">}{</span>$pword<span class="operator">} =
{}</span> unless exists $<span class="operator">{</span>$curr_node<span class="operator">}{</span>$pword<span class="operator">};</span>
$curr_node<span class="operator"> =</span> $<span class="operator">{</span>$curr_node<span class="operator">}{</span>$pword<span class="operator">};
}</span><span class="pre">
# put-the annotation into the last-node of the key-word tree
</span> $<span class="operator">{</span>$curr_node<span class="operator">}{</span>$ANNOT_DELIM<span class="operator">} =</span> $annotation<span class="operator">;
}
}
}</span><span class="pre">
#
# Build the failure function level-by level
#
#
</span>sub BuildFailureFunction<span class="operator">{</span>
my @params<span class="operator"> =</span> @_<span class="operator">;</span>
my $ac_root<span class="operator"> =</span> $params<span class="operator">[</span><span class="int">0</span><span class="operator">];</span>
my<span class="operator"> (</span>$k<span class="operator">,</span> @bfs_list<span class="operator">,</span> $keyword<span class="operator">);</span><span class="pre">
#
# Incremental algorithm; For level-1 its the $ac_root itself
#
</span> $<span class="operator">{</span>$ac_root<span class="operator">}{</span>$FAILED_DELIM<span class="operator">} =</span> $ac_root<span class="operator">;</span>
$<span class="operator">{</span>$ac_root<span class="operator">}{</span>$HEIGHT_DELIM<span class="operator">} =</span><span class="int"> 0</span><span class="operator">;</span>
foreach $k<span class="operator"> (</span>keys<span class="operator"> %</span>$ac_root<span class="operator">){</span><span class="flow">
if</span><span class="operator">((</span>$k eq $FAILED_DELIM<span class="operator">)</span><span class="operator"> or</span><span class="operator"> (</span>$k eq $ANNOT_DELIM<span class="operator">)</span><span class="operator">
or</span><span class="operator"> (</span>$k eq $HEIGHT_DELIM<span class="operator">)){</span>
next<span class="operator">;
}</span>
$<span class="operator">{</span>$<span class="operator">{</span>$ac_root<span class="operator">}{</span>$k<span class="operator">}}{</span>$FAILED_DELIM<span class="operator">} =</span> $ac_root<span class="operator">;</span>
$<span class="operator">{</span>$<span class="operator">{</span>$ac_root<span class="operator">}{</span>$k<span class="operator">}}{</span>$HEIGHT_DELIM<span class="operator">} =</span><span class="int"> 1</span><span class="operator">;</span>
push<span class="operator">(</span>@bfs_list<span class="operator">,</span> $<span class="operator">{</span>$ac_root<span class="operator">}{</span>$k<span class="operator">});
}</span><span class="flow">
while</span><span class="operator">(</span>$#bfs_list<span class="operator"> >=</span><span class="int"> 0</span><span class="operator">){</span>
$nref<span class="operator"> =</span> shift<span class="operator">(</span>@bfs_list<span class="operator">);</span> #node
foreach $keyword<span class="operator"> (</span>keys<span class="operator"> %</span>$nref<span class="operator">){</span><span class="flow">
if</span><span class="operator">((</span>$keyword eq $FAILED_DELIM<span class="operator">)</span><span class="operator"> or</span><span class="operator">
(</span>$keyword eq $ANNOT_DELIM<span class="operator">)</span><span class="operator"> or</span><span class="operator">
(</span>$keyword eq $HEIGHT_DELIM<span class="operator">)){</span>
next<span class="operator">;
}</span>
$nnref<span class="operator"> =</span> $<span class="operator">{</span>$nref<span class="operator">}{</span>$keyword<span class="operator">};</span>
$fpref<span class="operator"> =</span> $<span class="operator">{</span>$nref<span class="operator">}{</span>$FAILED_DELIM<span class="operator">};</span><span class="flow">
while</span><span class="operator">(!(</span>exists $<span class="operator">{</span>$fpref<span class="operator">}{</span>$keyword<span class="operator">})</span><span class="operator">
and</span> $fpref<span class="operator"> !=</span> $ac_root<span class="operator">){</span>
$fpref<span class="operator"> =</span> $<span class="operator">{</span>$fpref<span class="operator">}{</span>$FAILED_DELIM<span class="operator">};
}</span><span class="flow">
if</span><span class="operator">(</span>exists $<span class="operator">{</span>$fpref<span class="operator">}{</span>$keyword<span class="operator">}){</span>
$<span class="operator">{</span>$nnref<span class="operator">}{</span>$FAILED_DELIM<span class="operator">} =</span> $<span class="operator">{</span>$fpref<span class="operator">}{</span>$keyword<span class="operator">};
}</span><span class="flow">else</span><span class="operator">{</span>
$<span class="operator">{</span>$nnref<span class="operator">}{</span>$FAILED_DELIM<span class="operator">} =</span> $ac_root<span class="operator">;
}</span>
$<span class="operator">{</span>$nnref<span class="operator">}{</span>$HEIGHT_DELIM<span class="operator">} =</span> $<span class="operator">{</span>$nref<span class="operator">}{</span>$HEIGHT_DELIM<span class="operator">}+</span><span class="int">1</span><span class="operator">;</span>
push<span class="operator">(</span>@bfs_list<span class="operator">,</span> $nnref<span class="operator">);
}
}
}</span><span class="pre">
#
# Reads word delimited text from the stream or a file
#
</span>sub MatchWordDelimText<span class="operator">{</span>
my @params<span class="operator"> =</span> @_<span class="operator">;</span>
my $ac_root<span class="operator"> =</span> $params<span class="operator">[</span><span class="int">0</span><span class="operator">];</span>
my $stream<span class="operator"> =</span> $params<span class="operator">[</span><span class="int">1</span><span class="operator">];</span>
my<span class="operator"> (</span>$curr_state<span class="operator">,</span> $annot<span class="operator">,</span> $sword<span class="operator">,</span> $c<span class="operator">,</span> $height<span class="operator">);</span>
my $line<span class="operator">;</span><span class="flow">
while</span><span class="operator">(<</span>$stream<span class="operator">>){</span>
$line<span class="operator"> =</span> $_<span class="operator">;</span> chomp<span class="operator">(</span>$line<span class="operator">);</span>
$curr_state<span class="operator"> =</span> $ac_root<span class="operator">;</span>
@swords<span class="operator"> =</span> split<span class="operator">(/\</span>s<span class="operator">+/,</span> $_<span class="operator">);</span>
$c<span class="operator"> =</span><span class="int"> 0</span><span class="operator">;</span><span class="flow">
while</span><span class="operator">(</span>$c<span class="operator"> <=</span> $#swords<span class="operator">){</span>
$sword<span class="operator"> =</span> $swords<span class="operator">[</span>$c<span class="operator">];</span><span class="pre">
# Change the state of the automaton by reading the input
# print "I:$sword C:$c 1:$curr_state ";
</span><span class="flow"> if</span><span class="operator">(</span>exists $<span class="operator">{</span>$curr_state<span class="operator">}{</span>$sword<span class="operator">}){</span>
$curr_state<span class="operator"> =</span> $<span class="operator">{</span>$curr_state<span class="operator">}{</span>$sword<span class="operator">};</span>
$c<span class="operator">++;
}</span><span class="flow">else</span><span class="operator">{</span><span class="flow">
if</span><span class="operator">(</span>$curr_state<span class="operator"> ==</span> $ac_root<span class="operator">){</span>
$c<span class="operator">++;
}</span>
$curr_state<span class="operator"> =</span> $<span class="operator">{</span>$curr_state<span class="operator">}{</span>$FAILED_DELIM<span class="operator">};
}</span><span class="pre">
# print "2:$curr_state \n";
</span><span class="flow"> if</span><span class="operator">(</span>exists $<span class="operator">{</span>$curr_state<span class="operator">}{</span>$ANNOT_DELIM<span class="operator">}){</span>
$annot<span class="operator"> =</span> $<span class="operator">{</span>$curr_state<span class="operator">}{</span>$ANNOT_DELIM<span class="operator">};</span>
print<span class="string"> "$annot \n"</span><span class="operator">;
}
}
}
}</span>
sub StressTestMatchWordDelimText<span class="operator">{</span>
my @params<span class="operator"> =</span> @_<span class="operator">;</span>
my $ac_root<span class="operator"> =</span> $params<span class="operator">[</span><span class="int">0</span><span class="operator">];</span>
my $stream<span class="operator"> =</span> $params<span class="operator">[</span><span class="int">1</span><span class="operator">];</span>
my<span class="operator"> (</span>$curr_state<span class="operator">,</span> $annot<span class="operator">,</span> $sword<span class="operator">,</span> $c<span class="operator">,</span> $height<span class="operator">);</span>
my $line<span class="operator">;</span>
my<span class="operator"> (</span>$stress_me<span class="operator">,</span> $ss<span class="operator">);</span>
$line<span class="operator"> = <</span>$stream<span class="operator">>;</span>
print<span class="string"> "ABSTRACT: $line\n"</span><span class="operator">;</span>
chomp<span class="operator">(</span>$line<span class="operator">);</span>
$stress_me<span class="operator"> =</span><span class="int"> 1000000</span><span class="operator">;</span><span class="flow">
for</span><span class="operator">(</span>$s<span class="operator">=</span><span class="int">0</span><span class="operator">;</span> $s<span class="operator"><</span> $stress_me<span class="operator">;</span> $s<span class="operator">++){</span><span class="flow">
if</span><span class="operator">(</span>$s<span class="operator">%</span><span class="int">100000</span><span class="operator"> ==</span><span class="int"> 0</span><span class="operator">){</span>
print<span class="string"> "iteration $s \n"</span><span class="operator">;
}</span>
$curr_state<span class="operator"> =</span> $ac_root<span class="operator">;</span>
@swords<span class="operator"> =</span> split<span class="operator">(/\</span>s<span class="operator">+/,</span> $_<span class="operator">);</span>
$c<span class="operator"> =</span><span class="int"> 0</span><span class="operator">;</span><span class="flow">
while</span><span class="operator">(</span>$c<span class="operator"> <=</span> $#swords<span class="operator">){</span>
$sword<span class="operator"> =</span> $swords<span class="operator">[</span>$c<span class="operator">];</span><span class="pre">
# Change the state of the automaton by reading the input
# print "I:$sword C:$c 1:$curr_state ";
</span><span class="flow"> if</span><span class="operator">(</span>exists $<span class="operator">{</span>$curr_state<span class="operator">}{</span>$sword<span class="operator">}){</span>
$curr_state<span class="operator"> =</span> $<span class="operator">{</span>$curr_state<span class="operator">}{</span>$sword<span class="operator">};</span>
$c<span class="operator">++;
}</span><span class="flow">else</span><span class="operator">{</span><span class="flow">
if</span><span class="operator">(</span>$curr_state<span class="operator"> ==</span> $ac_root<span class="operator">){</span>
$c<span class="operator">++;
}</span>
$curr_state<span class="operator"> =</span> $<span class="operator">{</span>$curr_state<span class="operator">}{</span>$FAILED_DELIM<span class="operator">};
}</span><span class="pre">
# print "2:$curr_state \n";
</span><span class="flow"> if</span><span class="operator">(</span>exists $<span class="operator">{</span>$curr_state<span class="operator">}{</span>$ANNOT_DELIM<span class="operator">}){</span>
$annot<span class="operator"> =</span> $<span class="operator">{</span>$curr_state<span class="operator">}{</span>$ANNOT_DELIM<span class="operator">};</span><span class="pre">
# print "$annot \n";
</span><span class="operator"> }
}
}
}</span><span class="pre">
#
# Just to see how the tree looks -- for viewing pleasure :)
#
</span>sub PrintKwordTree<span class="operator">{</span>
my @params<span class="operator"> =</span> @_<span class="operator">;</span>
my $ac_root<span class="operator"> =</span> $params<span class="operator">[</span><span class="int">0</span><span class="operator">];</span>
my $level<span class="operator"> =</span> $params<span class="operator">[</span><span class="int">1</span><span class="operator">];</span>
my<span class="operator"> (</span>$k<span class="operator">,</span> $aref<span class="operator">,</span> $annot<span class="operator">,</span> $j<span class="operator">,</span> $nref<span class="operator">,</span> $fref<span class="operator">,</span> $height<span class="operator">);</span>
foreach $k<span class="operator"> (</span>keys<span class="operator"> %</span>$ac_root<span class="operator">){</span><span class="flow">
if</span><span class="operator">((</span>$k eq $ANNOT_DELIM<span class="operator">)</span><span class="operator">
or</span><span class="operator"> (</span>$k eq $FAILED_DELIM<span class="operator">)</span><span class="operator"> or</span><span class="operator"> (</span>$k eq $HEIGHT_DELIM<span class="operator">)){</span>
next<span class="operator">;
}</span><span class="flow">
for</span><span class="operator">(</span>$j<span class="operator">=</span><span class="int">0</span><span class="operator">;</span> $j<span class="operator"><</span>$level<span class="operator">;</span> $j<span class="operator">++){</span>
print<span class="string"> " "</span><span class="operator">;
}</span>
$fref<span class="operator"> =</span> $<span class="operator">{</span>$<span class="operator">{</span>$ac_root<span class="operator">}{</span>$k<span class="operator">}}{</span>$FAILED_DELIM<span class="operator">};</span>
$nref<span class="operator"> =</span> $<span class="operator">{</span>$ac_root<span class="operator">}{</span>$k<span class="operator">};</span>
$height<span class="operator"> =</span> $<span class="operator">{</span>$nref<span class="operator">}{</span>$HEIGHT_DELIM<span class="operator">};</span>
print<span class="string"> "-$k-> N=$nref F=$fref h=$height\n"</span><span class="operator">;</span>
PrintKwordTree<span class="operator">(</span>$<span class="operator">{</span>$ac_root<span class="operator">}{</span>$k<span class="operator">},</span> $level<span class="operator">+</span><span class="int">1</span><span class="operator">);
}</span><span class="flow">
if</span><span class="operator">(</span>exists $<span class="operator">{</span>$ac_root<span class="operator">}{</span>$ANNOT_DELIM<span class="operator">}){</span>
$annot<span class="operator"> =</span> $<span class="operator">{</span>$ac_root<span class="operator">}{</span>$ANNOT_DELIM<span class="operator">};</span><span class="flow">
for</span><span class="operator">(</span>$j<span class="operator">=</span><span class="int">0</span><span class="operator">;</span> $j<span class="operator"><</span>$level<span class="operator">;</span> $j<span class="operator">++){</span>
print<span class="string"> " "</span><span class="operator">;
}</span>
print<span class="string"> "; $annot ;\n\n"</span><span class="operator">;</span><span class="flow">
return</span><span class="operator">;
}
}</span><span class="pre">
#
# A simple unit-test
#
</span>sub<span class="keyword"> main</span><span class="operator">{</span>
BuildKwordTree<span class="operator">(\%</span>AC_TREE<span class="operator">,</span> PATTERN_FILE<span class="operator">);</span><span class="pre">
# PrintKwordTree(\%AC_TREE, 0);
</span> BuildFailureFunction<span class="operator">(\%</span>AC_TREE<span class="operator">);</span><span class="pre">
# PrintKwordTree(\%AC_TREE, 0);
# Now build the failure function.
</span> print<span class="string"> "please enter the words\n"</span><span class="operator">;</span><span class="pre">
# MatchWordDelimText(\%AC_TREE, STDIN);
</span> StressTestMatchWordDelimText<span class="operator">(\%</span>AC_TREE<span class="operator">,</span> TEXT_FILE<span class="operator">);</span>
print<span class="string"> "\n....THANK YOU......\n"</span><span class="operator">;
}</span><span class="keyword">
main</span><span class="operator">();</span></pre>
</body></html>Vamsi Kundetihttp://www.blogger.com/profile/06920946350286969217noreply@blogger.com0