Algorithms « Antimatroid, The

Archive for the ‘Algorithms’ Category

Minesweeper Agent

Introduction

Lately I’ve been brushing up on probability, statistics and machine learning and thought I’d play around with writing a Minesweeper agent based solely on these fields. The following is an overview of the game’s mechanics, verification of an implementation, some different approaches to writing the agent and some thoughts on the efficacy of each approach.

Minesweeper

Background

Minesweeper was created by Curt Johnson in the late eighties and later ported to Windows by Robert Donner while at Microsoft. With the release of Windows 3.1 in 1992, the game became a staple of the operating system and has since found its way onto multiple platforms and spawned several variants. The game has been shown to be NP-Complete, but in practice, algorithms can be developed to solve a board in a reasonable amount of time for the most common board sizes.

Specification

Gameplay
An agent, $\mathcal{A}$ , is presented a $n \times m$ grid containing $M$ uniformly distributed mines. The agent’s objective is to expose all the empty grid locations and none of the mines. Information about the mines’ grid locations is gained by exposing empty grid locations which will indicate how many mines exist within a unit (Chebyshev) distance of the grid location. If the exposed grid location is a mine, then the player loses the game. Otherwise, once all empty locations are exposed, the player wins.	$\textbf{PLAY-GAME}(\mathcal{A}, n, m, M) \newline \indent H \leftarrow \textbf{INITIALIZE-HIDDEN}(n,m,M) \newline \indent V \leftarrow \textbf{INITIALIZE-VISIBLE}(n,m,M) \newline \newline \indent \textbf{do} \newline \indent \indent (i,j) \leftarrow \mathcal{A}(V) \newline \indent \indent \textbf{EXPOSE}(H, V, i, j) \newline \indent \indent \textbf{if} \quad V_{i,j} = \text{*} \newline \indent \indent \indent \textbf{return} \quad \textit{Failure} \newline \indent \textbf{while} \quad M \neq \lvert \textbf{U-LOCATIONS}(V) \rvert \newline \newline \indent \textbf{return} \quad \textit{Success}$
Initialization
The board consists of hidden and visible states. To represent the hidden, $H$ , and visible state, $V$ , of the board, two character matrices of dimension $n \times m$ are used. Characters ’0′-’8′ represent the number of neighboring mines, character ‘U’ to represent an unexposed grid location and character ‘*’ for a mine. Neighbors of a grid location $(i,j)$ is the set of grid locations such that $1 \le \lVert (u,v) - (i,j) \rVert_{\infty} \le r$ . By default, $r = 1$ .	$\textbf{INITIALIZE-VISIBLE}(n, m) \newline \indent \textbf{return} \quad \textbf{INITIALIZE}(n,m, \text{U})$ $\textbf{INITIALIZE-HIDDEN}(n, m, M) \newline \indent H \leftarrow \textbf{INITIALIZE}(n,m,\text{0}) \newline \newline \indent S \leftarrow \emptyset \newline \indent \textbf{while} \quad \lvert S \lvert < M \newline \indent \indent S \leftarrow S \cup \textbf{RANDOM}(\textbf{LOCATIONS}(H)) \newline \newline \indent \textbf{foreach} \quad (i, j) \in S \newline \indent \indent H_{i, j} = \text{*} \newline \indent \indent \textbf{foreach} \quad (u,v) \in \textbf{NEIGHBORS}(n, m, i, j) \setminus S \newline \indent \indent \indent H_{u,v} \leftarrow H_{u,v} + 1 \newline \newline \indent \textbf{return} \quad H$
Exposing Cells
The expose behavior can be thought of as a flood fill on the grid, exposing any empty region bordered by grid locations containing mine counts and the boundaries of the grid. A matrix, $F \in \mathbb{Z}^{n \times m}$ , represents the topography of the board. A value of zero is reserved for sections of the board that have yet to be visited, a value of one for those that have, two for those that are boundaries and three for mines. A stack, $S$ , keeps track of locations that should be inspected. If a cell location can be exposed, then each of its neighbors will be added to the stack to be inspected. Those neighbors that have already been inspected will be skipped. Once all the reachable grid locations have been inspected, the process terminates.	$\textbf{EXPOSE}(H, V, i, j) \newline \indent \textbf{foreach} \quad (u,v) \in \textbf{LOCATIONS}(H) \newline \indent \indent \textbf{if} \quad H_{u,v} = \text{0} \newline \indent \indent \indent F_{u,v} \leftarrow 0 \newline \indent \indent \textbf{if} \quad H_{i,j} = \text{*} \newline \indent \indent \indent F_{u,v} \leftarrow 3 \newline \indent \indent \textbf{else} \newline \indent \indent \indent F_{u,v} \leftarrow 2 \newline \newline \indent \textbf{PUSH}(S, (i, j)) \newline \indent \textbf{do} \newline \indent \indent (u, v) \leftarrow \textbf{POP}(S) \newline \indent \indent \textbf{if} \quad F_{u,v} = 0 \newline \indent \indent \indent F_{u,v} \leftarrow 1 \newline \indent \indent \indent \textbf{foreach} \quad (r,s) \in \textbf{NEIGHBORS}(H, u, v) \newline \indent \indent \indent \indent \textbf{PUSH}(S, (r,s)) \newline \indent \indent \textbf{elseif} \quad F_{u,v} \in (1, 2) \newline \indent \indent \indent V_{u,v} \leftarrow H_{u, v} \newline \indent \textbf{while} \quad \textbf{COUNT}(S) > 0$

Verification

Methodology

Statistical tests are used to verify the random aspects of the game’s implementation. I will skip the verification of the game’s logic as it requires use of a number of different methods that are better suited for their own post.

There are two random aspects worth thinking about: the distribution of mines and the distribution of success (i.e., not clicking a mine) for random trials. In both scenarios it made since to conduct Pearson’s chi-squared test. Under this approach there are two hypotheses:

$H_{0}$ : The distribution of experimental data follows the theoretical distribution
$H_{a}$ : The distribution experimental data does not follow the theoretical distribution

$H_{0}$ is accepted when the test statistic, $\chi^2$ , is less than the critical value, $\chi_{\alpha}^2$ . The critical value is determined by deciding on a p-value (e.g., 0.05, 0.01, 0.001), $\alpha$ , that results in the tail area beneath the chi-squared distribution $\chi_{k}^2(x)$ equal to $\alpha$ . $k$ is the degrees of freedom in the observation.

Mine distribution

The first aspect to verify was that mines were being uniformly placed on the board. For a standard $9 \times 9$ board with $10$ mines, the expectation is that each grid location should be assigned $\frac{10}{81} N$ times for $N$ trials. $k = 80$ for this experiment.

In the above experiment, $\chi^2 = 71.78$ and $\chi_{0.05}^2 = 101.87$ . Since $\chi^2 < \chi_{0.05}^2$ , this affirms $H_{0}$ and that the implemented distribution of mines is indeed uniform with a statistical significance of $0.05$ .

Distribution of successful clicks

The second aspect to verify is that the number of random clicks before exposing a mine follows a hypergeometric distribution. The hypergeometric distribution is appropriate since we are sampling (exposing) without replacement (the grid location remains exposed after clicking). This hypothesis relies on a non-flood-fill exposure.

The distribution has four parameters. The first is the number of samples drawn (number of exposures), the second the number of successes in the sample (number of empty exposures), the third the number of successes in the population (empty grid locations) and the last the size of the population (grid locations): $h(nm-M;nm-M,nm-M, nm)$ .

The expected frequencies for the hypergeometric distribution is given by $N h(nm - M; nm - M, nm - M, nm)$ for $N$ trials. $k = 70$ in this case.

In the above experiment $\chi^2 = 47.24$ and $\chi_{0.05}^2 = 90.53$ . Since $\chi^2 < \chi_{0.05}^2$ , this affirms $H_{0}$ and that the number of locations exposed prior to exposing a mine follows a hypergeometric distribution with a statistical significance of $0.05$ .

Also included in the plot is the observed distribution for a flood based exposure. As one might expect, the observed frequency of more exposures decreases more rapidly than that of the non-flood based exposure.

Agents

Methodology

Much like how a human player would learn to play the game, I decided that each model would have knowledge of game’s mechanics and no prior experience with the game. An alternative class of agents would have prior experience with the game as the case would be in a human player who had studied other player’s strategies.

To evaluate the effectiveness of the models, each played against a series of randomly generated grids and their respective success rates were captured. Each game was played on a standard beginner’s $9 \times 9$ grid containing between $[1, 10]$ mines.

For those models that refer to a probability measure, $\mathbb{P}$ , it is assumed that the measure is determined empirically and treated as an estimate of the probability of an event and not as an a priori measure.

Marginal Model

Development

The first model to consider is the Marginal Model. It is designed to simulate the behavior of a naive player who believes that if he observes a mine at a grid location that the location should be avoid in future trials.

The model treats the visible board, $V$ , as a matrix of discrete random variables where each grid location is interpreted as either $\textit{Empty}$ or (a) $\textit{Mine}$ . This model picks the grid location with the greatest empirical probability of being empty:

$\underset{(i,j)}{\arg\max} \quad \mathbb{P}(X_{i,j} = \textit{Empty})$

$\textbf{MARGINAL-MODEL}(V) \newline \indent p_{\textit{max}} \leftarrow \perp \newline \indent (i, j)_{\textit{max}} \leftarrow \perp \newline \newline \indent \textbf{foreach} \quad (i,j) \in \textbf{U-LOCATIONS}(V) \newline \indent \indent p \leftarrow \mathbb{P}(V_{i,j} = \textit{Empty}) \newline \indent \indent \textbf{if} \quad p > p_{\textit{max}} \newline \indent \indent \indent p_{\textit{max}} \leftarrow p \newline \indent \indent \indent (i,j)_{\textit{max}} \leftarrow (i,j) \newline \newline \indent \textbf{return} \quad (i,j)_{\textit{max}}$

Test Results

Since the mine distribution is uniform, the model should be equivalent to selecting locations at random. The expected result is that avoiding previously occupied grid locations is an ineffective strategy as the number of mines increases. This does however, provide an indication of what the success rate should look like for chance alone.

Conditional Model

Development

One improvement over the Marginal Model is to take into account the visual clues made visible when an empty grid location is exposed. Since an empty grid location represents the number of neighboring mines, the Conditional Model can look at these clues to determine whether or not an unexposed grid location contains a mine.

This boils down to determining the probability of $\mathbb{P}(\textit{Mine} \lvert \textrm{Evidence})$ . A simplification in calculating the probability is to assume that each piece of evidence is independent. Under this assumption the result is a Naïve Bayes Classifier:

$\mathbb{P}( C = c \vert X = x ) = \dfrac{\mathbb{P}(C = c) \prod \mathbb{P}( X_{i} = x_{i} \vert C = c)}{\prod\mathbb{P}(X_{i} = x_{i})}$

As in the case of the Marginal Model, the Conditional Model returns the grid location that it has determined has the greatest probability of being empty given its neighbors:

$\underset{(i,j)}{\arg\max} \quad \mathbb{P}(C_{i,j} = \textit{Empty} | N(V_{i,j}) )$

$\textbf{CONDITIONAL-MODEL}(V, r) \newline \indent C \leftarrow \lbrace \textit{Empty}, \textit{Mine} \rbrace \newline \indent S \leftarrow \textbf{U-LOCATIONS}(V) \newline \newline \indent T \leftarrow \emptyset \newline \indent \textbf{foreach} \quad (i,j) \in S \newline \indent \indent N \leftarrow \textbf{NEIGHBORS}(V, i, j, r) \newline \indent \indent p_{\textit{max}} \leftarrow \perp \newline \indent \indent c_{\textit{max}} \leftarrow \perp \newline \indent \indent \textbf{foreach} \quad c \in C \newline \indent \indent \indent p \leftarrow \mathbb{P}(C = c) \newline \indent \indent \indent \textbf{foreach} \quad (u,v) \in N \newline \indent \indent \indent \indent p \leftarrow p * \mathbb{P}( X_{i, j} = V_{i, j} \vert C = c ) \newline \indent \indent \indent \textbf{if} \quad p > p_{\textit{max}} \newline \indent \indent \indent \indent p_{\textit{max}} \leftarrow p \newline \indent \indent \indent \indent c_{\textit{max}} \leftarrow c \newline \indent \indent \textbf{if} \quad c_{\textit{max}} = \textit {Empty} \newline \indent \indent \indent T \leftarrow T \cup (i, j) \newline \newline \indent \textbf{return} \quad \textbf{RANDOM}(T)$

Test Results

The Naïve Bayes Classifier is regarded as being an effective approach to classifying situations for a number of different tasks. In this case, it doesn’t look like it is effective at classifying mines from non-mines. The results are only slightly better than the Marginal Model.

Graphical Model

Development

One shortfall of the Conditional Model is that it takes a greedy approach in determining which action to take. A more sophisticated approach is to not just consider the next action, but the possible sequence of actions that will minimize the possibility of exposing a mine.

Each of the possible observable grids, $S$ , can be thought of as a set of vertices in a graph whose corresponding set of edges represent the transition between a state, $s$ , to the next observable state, $s^{\prime}$ . Each transition was achieved by performing an action, $a \in A$ , on the state. The specific action, $\alpha : S \times S \to A$ , is chosen from a subset of permitted actions given the state. Each transition has a probability, $\mathbb{P}(s^{\prime} \vert s)$ , of taking place.

It is possible to pick a path, $\pi$ , through this graph that minimizes the risk by assigning a reward, $\rho : S \to \mathbb{R}$ , to each state and attempting to identify an optimal path, $\pi^{*}_{s}$ , from the present state that yields the greatest aggregate reward,

$\displaystyle \varrho(\pi) = \sum_{(s, s^{\prime}) \in \pi} \rho(s^{\prime}) \mathbb{P}(s^{\prime} \vert s)$

Solving for $\pi^{*}_{s}$ is equivalent to solving the Longest Path Problem and can be computed efficiently using a dynamic programming solution.

$\displaystyle \pi_{s}^{*} \gets \underset{\pi_{s}}{\arg\max} \ \sum_{(\sigma, \sigma^{\prime}) \in \pi_{s}} \rho(\sigma^{\prime}) \mathbb{P}(\sigma^{\prime} \vert \sigma) \ \pi_{s} \in \Pi_{s}$

$\textbf{GRAPHICAL-MODEL}(V) \newline \indent (a, r)_{\textit{max}} \gets (\bot, \bot) \newline \indent T \gets \emptyset \newline \newline \indent \textbf{foreach} \quad U \in \textbf{SUB-GRIDS}(V) \newline \indent \indent (a, r)_{U} \gets \textbf{OPTIMAL-ACTION}(U, \bot) \newline \indent \indent \textbf{if} \quad r_{U} > r_{\textit{max}} \newline \indent \indent \indent (a,r)_{\textit{max}} \gets (a,r)_{U} \newline \indent \indent \indent T \gets \emptyset \newline \indent \indent \textbf{if} \quad r_{U} = r_{\textit{max}} \newline \indent \indent \indent T \gets T \cup (a,r)_{U} \newline \newline \indent (a, r)_{\textit{max}} \gets \textbf{RANDOM}(T) \newline \newline \indent \textbf{return} \quad a_{\textit{max}} \newline \newline \textbf{OPTIMAL-ACTION}(V) \newline \indent (a, r)_{\textit{max}} \gets (\bot, \bot) \newline \newline \indent \textbf{foreach} \quad (i,j) \in \textbf{U-LOCATIONS}(V) \newline \indent \indent \textbf{foreach} \quad V^{\prime} \in \textbf{OUTCOMES}(V, (i,j)) \newline \indent \indent \indent (a, r)_{V^{\prime}} \gets \textbf{OPTIMAL-ACTION}(V^{\prime}) \newline \indent \indent \indent r \gets r_{V^{\prime}} + \mathbb{P}(V^{\prime} \vert V) * \textbf{REWARD}(V^{\prime}) \newline \indent \indent \indent \textbf{if} \quad r \ge r_{\textit{max}} \newline \indent \indent \indent \indent (a, r)_{\textit{max}} \gets ((i,j), r) \newline \newline \indent \textbf{return} \quad (a,r)_{\textit{max}} \newline \newline \textbf{REWARD}(V) \newline \indent \textbf{if} \quad \textbf{SUCCESS}(V) \newline \indent \indent \indent \textbf{return} \quad +100 \newline \indent \textbf{if} \quad \textbf{FAILURE}(V) \newline \indent \indent \indent \textbf{return} \quad -100 \newline \newline \indent \textbf{return} +1$

From the optimal walk, a sequence of optimal actions is determined by mapping $\alpha$ over the path. Taking the first action gives the optimal grid location to expose given the current visible state of the board.

This description constitutes a Markov Decision Process. As is the case for most stochastic processes, it is assumed that the process holds the Markov Property; that future states only depend upon the current states and none of the prior states. In addition to being a Markov Decision Process, this is also an example of Reinforcement Learning.

First thing to observe is that the game state space is astronomical. For a standard beginner’s grid there is at most a sesvigintillion $(10^{81})$ possible grids that a player can encounter. Which as an aside, is on the order of the number of atoms in the observable universe! The set of actions at each state is slightly more manageable with at most eighty-one actions.

To simplify the state space, I chose to only consider $3 \times 3$ boards and when evaluating a full grid, consider the possible sub-grids and evaluate the optimal sequence of actions for each sub-grid and pick the maximum reward associated for each sub-grid that was evaluated as the action to take on the full grid.

Test Results

The Graphical Model produces results that are only a margin better than those of the Conditional Model.

Semi-deterministic Model

Development

The last model I’m going to talk about is a semi-deterministic model. It works by using the visible grid to infer the topology of the hidden grid and from the hidden grid, the topology that the visible grid can become.

The grid can be viewed as a graph. Each grid location is a vertex and an edge is an unexposed grid location’s influence on another grid location’s neighbor mine count.

For each of the exposed grid locations on the board, $v_{i,j}$ , it’s neighbors, $N(v_{i,j})$ , are all mines when the number of inbound edges $E_{i,j} = d(v_{i,j})$ , matches the visible mine count $V_{i,j}$ .

The model produces its inferred version, $F$ , of the influence graph $E$ by using the determined mine locations $M$ .

For each of the grid locations that are exposed and the inferred influence matches the visible count, then each of the neighbors about that location can be exposed provided they are not already exposed and not an inferred mine.

From this set of possibilities, a mine location is chosen. When no mine locations can be determined, then an alternative model can be used.

$\textbf{SEMI-DETERMINISTIC-MODEL}(V) \newline \indent E \leftarrow 0_{n, m} \newline \indent \textbf{foreach} \quad (i,j) \in \textbf{LOCATIONS}(V) \newline \indent \indent \textbf{if} \quad V_{i,j} \ne \textit{U} \newline \indent \indent \indent \textbf{continue} \newline \indent \indent \textbf{foreach} \quad (u,v) \in \textbf{NEIGHBORS}(V, i, j) \newline \indent \indent \indent E_{u,v} \leftarrow E_{u,v} + 1 \newline \newline \indent M \leftarrow \textit{False}_{n,m} \newline \indent \textbf{foreach} \quad (i,j) \in \textbf{LOCATIONS}(V) \newline \indent \indent \textbf{if} \quad V_{i,j} \in \lbrace \textit{U}, \textit{0} \rbrace \lor V_{i,j} \neq E_{i,j} \newline \indent \indent \indent \textbf{continue} \newline \indent \indent \textbf{foreach} \quad (u,v) \in \textbf{NEIGHBORS}(V, i, j) \newline \indent \indent \indent \textbf{if} \quad V_{i,j} = \textit{U} \newline \indent \indent \indent \indent M_{u,v} \leftarrow \textit{True} \newline \newline \indent F \leftarrow 0_{n,m} \newline \indent \textbf{foreach} \quad (i,j) \in \textbf{LOCATIONS}(V) \newline \indent \indent \textbf{foreach} \quad (u,v) \in \textbf{NEIGHBORS}(V, i, j) \newline \indent \indent \indent \textbf{if} \quad M_{u,v} \newline \indent \indent \indent \indent F_{i,j} \leftarrow F_{i,j} + 1 \newline \newline \indent S \leftarrow \emptyset \newline \indent \textbf{foreach} \quad (i,j) \in \textbf{LOCATIONS}(V) \newline \indent \indent \textbf{if} \quad V_{i,j} = \textit{U} \lor F_{i,j} \ne V_{i,j} \newline \indent \indent \indent \textbf{continue} \newline \indent \indent \textbf{foreach} \quad (u,v) \in \textbf{NEIGHBORS}(V, i, j) \newline \indent \indent \indent \textbf{if} \quad V_{i,j} \ne \textit{U} \lor M_{u,v} \newline \indent \indent \indent \indent \textbf{continue} \newline \indent \indent \indent S \leftarrow S \cup (u, v) \newline \newline \indent \textbf{return} \quad \textbf{RANDOM}(S)$

Test Results

Since the model is a more direct attempt at solving the board, its results are superior to the previously presented models. As the number of mines increases, it is more likely that it has to rely on a more probabilistic approach.

Summary

Each of the models evaluated offered incremental improvements over their predecessors. Randomly selecting locations to expose is on par with choosing a location based on previously observed mine locations. The Conditional Model and Graphical Model yield similar results since they both make decisions based on conditioned probabilities. The Semi-deterministic Model stands alone as the only one model that produced reliable results.

The success rate point improvement between the Condition and Marginal models is most notable for boards consisting of three mines and the improvement between Graphical and Semi-deterministic models for seven mines. Improvements between Random and Marginal models is negligible and between Conditional and Graphical is minor for all mine counts fewer than seven.

Given the mathematical complexity and nondeterministic nature of the machine learning approaches, (in addition the the complexity and time involved in implementing those approaches) they don’t seem justified when more deterministic and simpler approaches exist. In particular, it seems like most people have implemented their agents using heuristics and algorithms designed to solve constraint satisfaction problems. Nonetheless, this was a good refresher to some of the elementary aspects of probability, statistics and machine learning.

References

“Classification – Naïve Bayes.” Data Mining Algorithms in R. Wikibooks. 3 Nov. 2010. Web. 30 Oct. 2011.

“Windows Minesweeper.” MinesweeperWiki. 8 Sept. 2011. Web. 30 Oct. 2011.

Kaye, Richard. “Minesweeper Is NP-complete.” [pdf] Mathematical Intelligencer 22.2 (2000): 9-15. Web. 30 Oct. 2011.

Nakov, Preslav, and Zile Wei. “MINESWEEPER, #MINESWEEPER.” 14 May 2003. Web. 14 Apr. 2012.

Richard, Sutton, and Andrew G. Barto. “3.6 Markov Decision Processes.” Reinforcement Learning: An Introduction. Cambridge, Massachusetts: Bradford Book, 1998. 4 Jan. 2005. Web. 30 Oct. 2011.

Rish, Irene “An Empirical Study of the Naive Bayes Classifer.” [pdf] IJCAI-01 Workshop on Empirical Methods in AI (2001). Web. 30 Oct. 2011.

Russell, Stuart J., and Peter Norvig. Artificial Intelligence: A Modern Approach. Upper Saddle River, NJ: Prentice Hall/PearsonEducation., 2003. Print.

Sun, Yijun, and Jian Li. “Adaptive Learning Approach to Landmine Detection.” [pdf] IEEE Transactions of Aerospace and Electronic Systems 41.3 (2005): 1-9. 10 Jan. 2006. Web. 30 Oct. 2011.

Taylor, John R. An introduction to error analysis: the study of uncertainties in physical measurements. Sausalito, CA: University Science Books, 1997. Print.

Written by lewellen

2012-05-01 at 8:00 am

Posted in Algorithms, Artificial Intelligence, Graph Theory, Probability and Statistics

Tagged with Chi-squared Test, Flood Fill, Hypergeometric Distribution, Longest Path Problem, Markov Decision Process, Minesweeper, Naive Bayes Classifier, Reinforcement Learning

Tree Drawing: Force-based Algorithm

with 18 comments

I’ve written a couple times on different data visualization techniques for viewing hierarchical information. One method that I haven’t discussed is the force-based algorithm approach. After digging around a bit, it looks like this method was originally developed in Birbil’s and Fang’s 2003 publication [pdf], “An Electromagnetism-like Mechanism for Global Optimization”. In terms of tree drawing, the idea is fairly simple: treat the tree as a dynamical system composed of springs and point charges, apply the appropriate physical laws and over time, you will get a nicely laid out tree. My implementation in the clip below demonstrates this approach.

Let’s think about why we want a dynamical system. An idealized tree drawing is aesthetically pleasing, compact, nodes are uniformly distributed and edges between nodes do not cross. Under this proposed system, we can’t guarantee these criteria, but we can get “close enough” to produce something that is both nice on paper and sensible in practice. Treating the nodes as point charges will make the nodes repel from one another. The magnitude of the repulsion is determined by Column’s Law which states that the force applied to a point charge by another point charge follows an inverse square law. If all we had was Coulomb’s Law, then all the nodes will fly away from one another and finally come to rest, but the nodes would be too distantly spread across the screen. To alleviate this, we can treat the edges between nodes as springs. Following Hooke’s Law we can tether the nodes together to preserve the compactness criteria. The magnitude of this force is directly linear with respect to the distance between the nodes.

Let’s jump into some quick physics. To determine the net force applied on a given node, we need to calculate the spring force applied to the node, and then we need to calculate all the point charge forces applied to the node. The key to getting this right is making sure that we get the direction and magnitude of the forces correct. To make sure that we have the direction correct, let $\bold{d}_{i,j} = \bold{x}_{i} - \bold{x}_{j}$ be the vector representing the distance between the i’th and j’th nodes ( $\bold{x}_{i}$ representing the i’th node’s location). Let $\hat{\bold{d}}_{i,j} = \frac{\bold{d}_{i,j}}{\lVert \bold{d}_{i,j} \rVert}$ be the normalized vector representing the direction that the force that needs to be applied. For the spring force we are looking at $\bold{F}_{i, j} = -k (\lVert \bold{d}_{i,j} \rVert - r) \hat{\bold{d}}_{i,j}$ where $k$ is the spring constant, $r$ is the length of the spring at rest and $\bold{F}_{i,j}$ is the force between the i’th and j’th node. For the point charge we are interested in $\bold{F}_{i,j} = k \frac{ Q_{i} Q_{j} }{ \lVert \bold{d}_{i,j} \rVert ^{2} } \hat{\bold{d}}_{i,j}$ . $Q_{i}$ is the charge of the i’th node and $k$ is the coulomb constant. (I default all constants to 1.0 in my implementation).

We’ve got the net force on each node figured out, so now it is necessary to figure out the velocity of the nodes as well as their location. To do this, we recall that by Newton’s Laws, a force is the product of a mass and acceleration $F = ma$ . If we know an acceleration, we can derive (integrate mathematically) the velocity and location each iteration. Thus, if we know a node’s present velocity $v_{n}$ , then we know that $v_{n+1} = a h + v_{n}$ . If the node has current location $p_{n}$ , then its next location will be $p_{n+1} = \frac{h^{2}}{2} a + h v_{n} + p_n$ . Once we’ve determined these values between time steps $h$ , we will need to remember to zero out the net force on each node in our implementation. In my implementation I chose a value of 0.1 for my time step) One important aspect to note for this algorithm, is the need to apply dampening to the velocity over time, otherwise the structure continue to move in space and we run into the possibility of running out of numbers on the machine resulting in an overflow error.

To go about implementing all of this you’ll find the time complexity of the algorithm should be on the order $O(\lVert{V}\rVert^2)$ . The square order comes from apply Coulomb’s law. By Newton’s third law (equal and opposite forces), you can reduce this down to $\frac{ \lVert V \lVert^2 - \lVert V \lVert }{2}$ , but for significantly sized $V$ you are still going to be looking at $O(\lVert{V}\rVert^{2})$ . You could improve this by applying a binary space partitioning (in particular the R-tree and its variants) data structure to the tree and only use those nodes within a certain radius of the node when calculating the Coulomb forces. Testing my implementation (using Lists), I found things get a little sluggish for a real-time render around 100 nodes on my 2.4GHz laptop. Depending on your application (domain and implementation), you millage may vary. Following are some test cases I used in writing my solution.

Fan	Flat

Each node a number of nodes equivalent to its number of siblings minus one	Tree containing one level of nodes
Complete Binary	Random

Every node has exactly two children	Randomly generated tree of depth four

I chose to implement my solution in C# using the Microsoft .net 3.5 Framework (as a side note, 4.0 was recently released on the 12th of April). I went with my custom mathematics library and the WinForms library of the Framework (in the future I would like to migrate this over to WPF) for displaying and end-user manipulation of the algorithm’s variables. I went with a vanilla UserControl that contains a setter for the tree data structure, the algorithm to apply to the structure and the settings for that algorithm. Upon setting the tree, a Timer instance kicks off (non-reentrant) to invoke the algorithm below every 50ms. To get the tree to display correctly on the control, the algorithm calculates the boundary of the nodes and then the control performs a linear map from the model space to the view space. The first implementation used the DoubleBuffered property of the UserControl, but I found that it was ineffective in reducing flickering so I implemented custom double buffering using the BufferedGraphicsContext class. It’s worth noting that most implementations track the kinetic energy of the system to determine when to terminate the algorithm. I chose not to do this, as I didn’t find value in adding in the additional calculation.

using System;
using System.Collections.Generic;
using System.Drawing;
using Library.Mathematics;

namespace ForceBasedTreeLayout {
    public class ForceBasedTreeLayout {
        private LayoutSettings settings;

        public ForceBasedTreeLayout(LayoutSettings settings) {
            this.settings = settings;
        }

        public RectangleF Evaluate(Node root) {
            List<Node> nodes = new List<Node>();
            foreachNode(root, (x) => nodes.Add(x));

            // Apply Hooke's law
            // F = -k x
            foreachNode(root, (parent) => {
                parent.Children.ForEach((child) => {
                    Vector dist = parent.Location - child.Location;
                    Vector restingDist = settings.MinimumSpringLength * dist.Normalize();
                    Vector force = -settings.SpringConstant * (dist - restingDist);

                    parent.Acceleration += force;
                    child.Acceleration -= force;
                });
            });

            // Apply Coulomb's Law
            // F = Q1 Q1 / d^2
            for (int n = 0; n < nodes.Count; n++) {
                for (int m = n + 1; m < nodes.Count; m++) {
                    Vector dist = nodes[n].Location - nodes[m].Location;
                    Vector norm = dist.Normalize();
                    Vector force = new Vector(2, (i) => norm[i] / Math.Pow(dist.Norm() + 1.0, 2.0));

                    nodes[n].Acceleration += force;
                    nodes[m].Acceleration -= force;
                }
            }

            Vector xExtrema = new Vector(2);
            xExtrema[0] = double.MaxValue;
            xExtrema[1] = double.MinValue;

            Vector yExtrema = new Vector(2);
            yExtrema[0] = double.MaxValue;
            yExtrema[1] = double.MinValue;

            // update the locations & velocity && figure out new bounding region
            foreach (Node node in nodes) {
                // p = a0t^2/2 + v0t + p0
                node.Location = (settings.TimeStep * settings.TimeStep * 0.5) * node.Acceleration + (settings.TimeStep) * node.Velocity + node.Location;

                // v = at + v0
                node.Velocity = (settings.TimeStep) * node.Acceleration + node.Velocity;
                node.Velocity = (settings.VelocityDampening) * node.Velocity;
                
                node.Acceleration = new Vector(2, (i) => 0.0);

                xExtrema[0] = Math.Min(xExtrema[0], node.Location[0]);
                xExtrema[1] = Math.Max(xExtrema[1], node.Location[0]);

                yExtrema[0] = Math.Min(yExtrema[0], node.Location[1]);
                yExtrema[1] = Math.Max(yExtrema[1], node.Location[1]);
            }

            RectangleF R = new RectangleF();
            R.X = (float)xExtrema[0];
            R.Y = (float)yExtrema[0];
            R.Width = (float)(xExtrema[1] - xExtrema[0]);
            R.Height = (float)(yExtrema[1] - yExtrema[0]);

            R.X -= R.Width / 2;
            R.Y -= R.Height / 2;
            R.Width *= 2;
            R.Height *= 2;

            return R;
        }

        private void foreachNode(Node root, Action<Node> act) {
            Stack<Node> stack = new Stack<Node>();
            stack.Push(root);
            while (stack.Count > 0) {
                Node node = stack.Pop();
                act(node);
                node.Children.ForEach((x) => stack.Push(x));
            }
        }
    }
}

Edit: 2010-10-21

By popular demand, here is the vector class:

public class Vector {
	private double[] V;

	public int Dimension { get { return V.Length; } }

	public double this[int n] {
		get {
			if (n < 0 || n >= Dimension)
				throw new Exception(string.Format("{0} must be between 0 and {1}", n, Dimension));
			return V[n];
		}
		set {
			if (n < 0 || n >= Dimension)
				throw new Exception(string.Format("{0} must be between 0 and {1}", n, Dimension));
			V[n] = value;
		}
	}

	public Vector() : this(0) { }

	public Vector(int n) { V = new double[n]; }

	public Vector(int n, VectorInitializer initializer) : this(n) {
		for (int i = 0; i < Dimension; i++)
			V[i] = initializer(i);
	}

	public double dot(Vector y) {
		if (Dimension != y.Dimension)
			throw new Exception();
		double d = 0.0;
		for (int n = 0; n < Dimension; n++)
			d += this[n] * y[n];
		return d;
	}

	public override bool Equals(object obj) {
		Vector x = obj as Vector;
		if (x == null || x.Dimension != Dimension)
			return false;

		for (int n = 0; n < Dimension; n++)
			if (this[n] != x[n])
				return false;
		return true;
	}

	static public Vector operator +(Vector x, Vector y) {
		if (x.Dimension != y.Dimension)
			throw new Exception();

		Vector z = new Vector(x.Dimension);
		for (int n = 0; n < z.Dimension; n++)
			z[n] = x[n] + y[n];
		return z;
	}

	static public Vector operator -(Vector x, Vector y) {
		if (x.Dimension != y.Dimension)
			throw new Exception();

		Vector z = new Vector(x.Dimension);
		for (int n = 0; n < z.Dimension; n++)
			z[n] = x[n] - y[n];
		return z;
	}

	static public double operator *(Vector x, Vector y) {
		return x.dot(y);
	}

	static public Vector operator *(double c, Vector x) {
		Vector y = new Vector(x.Dimension);
		for (int n = 0; n < y.Dimension; n++)
			y[n] = c*x[n];
		return y;
	}
}

Written by lewellen

2010-05-01 at 8:00 am

Posted in Algorithms, Computer Graphics

Tagged with C-Sharp, Dynamical System, Force-based Algorithm, Tree Drawing

Sudoku Solver in Haskell

with 6 comments

This month will be a bit of short article since I haven’t had a whole lot of time on my hands lately. Haskell is a wonderful little language that has begun to pick up a bit of moment in the past year that I’ve been playing with on-and-off now for several years. Since I don’t post enough on Haskell, I figured I’d post my bare-bones Haskell Sudoku Solver.

import Data.List
import Data.Maybe

toRowColumn :: Int -> (Int, Int)
toRowColumn index = (r, c)
	where	r = div index 9
		c = mod index 9

toIndex :: (Int, Int) -> Int
toIndex (r, c) = r * 9 + c

toRegion :: (Int, Int) -> Int
toRegion (r, c) = (div r 3) * 3 + (div c 3)

columnIndicies :: Int -> [Int]
columnIndicies c = [c, c + 9..80]

regionIndicies :: Int -> [Int]
regionIndicies g = [toIndex(r + x, c + y) |
		x <- [0..2],
		y <- [0..2]
	]
	where	r = (div g 3) * 3
		c = (mod g 3) * 3

rowIndicies :: Int -> [Int]
rowIndicies r = [9 * r..9 * (r + 1) - 1]

values :: [Int] -> [Int] -> [Int]
values board indicies = filter (>0) (map (board!!) indicies)

possibleValues :: [Int] -> (Int, Int) -> [Int]
possibleValues board rowColumn = 
	foldl (\\) [1..9] (
		map (values board) (
			map (\f -> (fst f . snd f) rowColumn) (
				zip [rowIndicies, columnIndicies, regionIndicies]
					[fst, snd, toRegion]
			)
		)
	)

validBoard :: [Int] -> Bool
validBoard board = (length board == 81) && (and $ map (==0) l)
	where	l = map length s
		s = map (\\[1..9]) v
		v = [values board (xIndicies x) |
			x <- [0..8],
			xIndicies <- [rowIndicies, columnIndicies, regionIndicies]]

solvedBoard :: [Int] -> Bool
solvedBoard board = and $ map (>0) board

hasUnassigned :: [Int] -> Bool
hasUnassigned board = isJust $ elemIndex 0 board

assignFirstUnassigned :: [Int] -> Int -> [Int]
assignFirstUnassigned (b:bs) value
	| b == 0 = value : bs
	| otherwise = b : (assignFirstUnassigned bs value)

possibleBoards :: [Int] -> [Int] -> [[Int]]
possibleBoards board possibleAssignments = map (assignFirstUnassigned board)
	possibleAssignments

solve :: [Int] -> [[Int]]
solve board
	| not (validBoard board) = [[]]
	| solvedBoard board = [board]
	| not (hasUnassigned board) = [[]]
	| otherwise = concated
	where	concated = concat mapped
		mapped = map solve filtered
		filtered = filter (not . null) possibleSolved
		possibleSolved = possibleBoards board possibleAssignments
		possibleAssignments = possibleValues board unassignedRowColumn
		unassignedRowColumn = toRowColumn unassignedIndex
		unassignedIndex = fromJust $ elemIndex 0 board

demo :: [Int]
demo = [2,0,0,0,8,0,3,0,0,
	0,6,0,0,7,0,0,8,4,
	0,3,0,5,0,0,2,0,9,
	0,0,0,1,0,5,4,0,8,
	0,0,0,0,0,0,0,0,0,
	4,0,2,7,0,6,0,0,0,
	3,0,1,0,0,7,0,4,0,
	7,2,0,0,4,0,0,6,0,
	0,0,4,0,1,0,0,0,3]

Written by lewellen

2009-06-01 at 12:00 am

Posted in Algorithms

Tagged with Haskell, Sudoku

Sudoku Solver in C# using Lambda Expressions

leave a comment »

Seems that everywhere you look someone has a Sudoku Solver that they want to showcase, well, I’m no different so I figured I’d post my take on the subject. Microsoft has introduced/included/borrowed a number of functional programming features into the latest version of C# (3.0) that have made it easier for developers to write better, cleaner code. One of those features continues the trend of improving anonymous methods, which extend delegates which extend interfaces now known as Lambda expressions. E.g., the following are all the same for the expression $\lambda x . x * x$ :

Func<int, int> square = (x) => x * x;

Func<int, int> square = new Func<int, int>(delegate(int x) { return x * x; });

public interface Func<T, R> {
    R Evaluate(T x);
}

public class Square : Func<int, int> {
    public int Evaluate(int x) {
        return x*x;
    }
}

...

Func<int,int> square = new Square();

Given that level of expressive power, I thought I would approach this implementation using as many lambda expressions as possible to see how concise and easy to follow an implementation I could create.

To start off, the solver will assume that the board will be passed in as a 81 character array containing digits 0-9. Zero shall represent an unassigned cell.

The core loop is pretty simple: Start with the initial board on a stack, in a loop take the first board off the stack and check if it is valid. A board is said to be valid if for every row, column and region each structure contains one and only one instance of the digits 1-9. If the board is not valid, no further work should be done.

Next we need to check if the board is solved. If it is not, then we need to explore the possible boards that can be derived from that board. To do so, we will need to find the first possible cell that is unassigned and push on to the the stack a derived board using the possible values that can be assigned to the identified cell.

Once the loop finally exits, print out the solved board.

using System;
using System.Collections.Generic;

namespace Sudoku {
    public class Program {
        static public void Main(string[] args) {
            string input = "200080300060070084030500209000105408000000000402706000301007040720040060004010003";
            Board inputBoard = null;
            try {
                inputBoard = new Board(input);
            } catch (InvalidInputException iie) {
                Console.WriteLine(iie.InvalidInput);
                return;
            }

            Stack boards = new Stack();
            boards.Push(inputBoard);
            Board board = null;
            do {
                board = boards.Pop();
                if (!board.Valid)
                    continue;
                if (!board.Solved)
                    board.FirstAvailable((r, c) => 
                        board.PossibleValuesAt(r, c, (v) => 
                            boards.Push(board.DeriveUsing(r, c, v))
                        )
                    );
            } while (boards.Count > 0 && !board.Solved);
            Console.WriteLine(board);
        }
    }
}

I’m going to start with the private member methods since they form the basic grammar that I will use to implement the public member methods and properties.

The first thing you’ll notice is the private class Structure; it is a simple pair class that contains two member properties for accessing a function that iterates over all possible structures and a function that iterates over all the cells in a specific structure.

The private constructor instantiates a private member array containing the function pointers that enumerate over Rows, Column and Regions.

assertStructure method which iterates over every instance of a structure in the table, and verifies that one and only one instance of the digits 1-9 exist in that structure instance.

indexInStructure iterates over all of the indices of a structure- in this case, 0-8.

using System;
using System.Text;

namespace Sudoku {
    public partial class Board {
        private class Struct {
            public Action<Action<int>> In { get; set; }
            public Action<int, Action<int>> ForValues { get; set; }
            public Struct(Action<Action<int>> _in, Action<int, Action<int>> forValues) {
                In = _in; ForValues = forValues;
            }
        }

        private Struct[] structures;

        private Board(int[] board) {
            this.board = board;
            structures = new Struct[] {
                new Struct(Rows, ValuesInRow), new Struct(Columns, ValuesInColumn),
                new Struct(Regions, ValuesInRegion)
            };
        }

        private bool assertStructure(Action<Action<int>> structure, Action<int, Action<int>> 
            valuesInStructure) {
            bool asserted = true;
            structure((x) => {
                int[] used = new int[10];
                valuesInStructure(x, (y) => used[y]++);
                used[0] = 0;
                asserted &= forAllValues((v) => v < 2, used);
            });
            return asserted;
        }

        private void indexInStructure(Action<int> actOnIndex) {
            indexInStructure((n) => true, actOnIndex);
        }

        private void indexInStructure(Predicate<int> p, Action<int> actOnIndex) {
            indexFromToWhere(0, 9, p, actOnIndex, false);
        }

The forAllValues method is simply a way of writing the predicate $(\forall x) Px$ applied to members of the universe of discourse (the array that was passed in).

indexFromToWhere is a simple wrapper around a common for loop with filtering and the option to break after the first object to pass through the filter is found.

        private bool forAllValues(Predicate f, T[] A) {
            bool held = true;
            for (int n = 0; n < A.Length && held; n++)
                held &= f(A[n]);
            return held;
        }

        private void indexFromToWhere(int min, int max, Predicate<int> p, Action<int> actOnIndex, 
            bool breakAfterFirst) {
            for (int n = min; n < max; n++)
                if (p(n)) {
                    actOnIndex(n);
                    if (breakAfterFirst)
                        break;
                }
        }

The following methods are all related to working with the board representation. I choose to implement the board as an integer array. The at method maps a logical row and column to a row order array value within board. Each of the different indexInBoard methods allows for iterating over the indices of the board.

        private int[] board;

        private int at(int row, int col) {
            return row * 9 + col;
        }

        private void indexInBoard(Action<int> actOnIndex) {
            indexInBoard((n) => true, actOnIndex);
        }

        private void indexInBoard(Predicate<int> p, Action<int> actOnIndex) {
            indexInBoard(p, actOnIndex, false);
        }

        private void indexInBoard(Predicate<int> p, Action<int> actOnIndex, bool breakAfterFirst) {
            indexFromToWhere(0, board.Length, p, actOnIndex, breakAfterFirst);
        }
    }
}

The first set of public member methods and properties we can look at are for managing the state of the board. The board is only Solved if each index is assigned. The board is only Valid if every structure in the board is asserted to be true. The constructor checks the input to make sure it is valid, delegates the some work to the private constructor and loads the input string in to the integer array.

using System;
using System.Text;

namespace Sudoku {
    public partial class Board {
        public bool Solved {
            get { 
                return forAllValues((x) => x > 0, board);
            }
        }
        public bool Valid { 
            get { 
                return forAllValues((s) => assertStructure(s.In, s.ForValues), structures);
            }
        }
        
        public Board(string input) : this (new int[81]) {
            if (string.IsNullOrEmpty(input))
                throw new InvalidInputException(InvalidInput.Empty);
            if (input.Length != 81)
                throw new InvalidInputException(InvalidInput.Length);
            for (int n = 0; n  board[n] = input[n] - '0');
        }

        override public string ToString() {
            StringBuilder S = new StringBuilder(board.Length);
            indexInBoard((n) => S.Append((char)('0' + board[n])));
            return S.ToString();
        }

Next, we need a way to operate on each of the structures in the board. Each row from top to bottom, each column from left to right and each region from top left to bottom right (zig-zag) will be assigned an index from 0-8 respectively.

        public void Columns(Action<int> actOnColumn) {
            indexInStructure(actOnColumn);
        }

        public void Regions(Action<int> actOnRegion) {
            indexInStructure(actOnRegion);
        }

        public void Rows(Action<int> actOnRow) {
            indexInStructure(actOnRow);
        }

It is easy to then iterate the cells in a given structure. The values in a column are simply a trip down the Rows and we only want to act when the value is defined. The values in a row are just as easy by traveling across the Columns and acting when the value is defined. Iterating over the values in the region is a little more involved, but nonetheless just as easy to follow- map the region to the appropriate reference row and column and then iterate over the 3×3 grid and act only when the value is defined.

        public void ValuesInColumn(int column, Action<int> actOnValue) {
            Rows((r) => {
                int value = board[at(r, column)];
                if (value > 0)
                    actOnValue(value);
            });
        }

        public void ValuesInRegion(int region, Action<int> actOnValue) {
            int row = (region / 3) * 3;
            int column = (region % 3) * 3;
            int value = 0;
            for (int r = 0; r < 3; r++)
                for (int c = 0; c  0)
                        actOnValue(value);
                }
        }

        public void ValuesInRow(int row, Action<int> actOnValue) {
            Columns((c) => {
                int value = board[at(row, c)];
                if (value > 0)
                    actOnValue(value);
            });
        }

Finally, we have the interesting methods used in the main loop. DeriveUsing will copy the board into a cloned integer array and then assign the derived at row and column with the value supplied. The newly derived board is then returned.

FirstAvailable iterates over all the indices of the board until it finds an unassigned value and then it acts upon the reference row and column.

PossibleValuesAt goes and collects a list of possible values by first collecting the values used in the row, column and region that the reference row and column reside within, then each value not found is acted upon.

        public Board DeriveUsing(int row, int colum, int withValue) {
            int[] derived = new int[81];
            indexInBoard((n) => derived[n] = board[n]);
            derived[at(row, colum)] = withValue;
            return new Board(derived);
        }

        public void FirstAvailable(Action<int, int> actOnRowAndColumn) {
            indexInBoard((n) => 
                board[n] == 0, (n) => actOnRowAndColumn(n / 9, n % 9), true);
        }

        public void PossibleValuesAt(int row, int column, Action<int> actOnPossibleValue) {
            bool[] used = new bool[10];
            ValuesInRow(row, (x) => used[x] = true);
            ValuesInColumn(column, (x) => used[x] = true);
            ValuesInRegion((row / 3) * 3 + (column / 3), (x) => used[x] = true);
            indexFromToWhere(1, used.Length, (n) => !used[n], actOnPossibleValue, false);
        }
    }
}

Having approached this implementation as I did, I found some interesting bugs that I hadn’t come across before and I figure I’ll close with one that caught me off guard. I had spent an hour writing all my code and figured I’d run it to see what kind of output I got. To my surprise I got an immediate StackOverflowExeception. So I spent and an additional 10 minutes debugging and found the following offending code. Take a look at it and see if you can see what’s wrong with it before reading on.

public void act(Func<int,int> assign){ 
    act((n) => board[n] = assign(n));
}

public void act(Action<int> actOn) {
    for(int n = 0; n < board.Length; n++)
      actOn(n);
}

After a minute I had my ah-ha moment (as I’m sure have as well) and remembered that the assignment operator returns the value of the assignment so the compiler was turning (n) => board[n] = assign(n), into Func<int, int> instead of Action<int> as I had hoped. To fix the bug, I had to do the following to get the compiler to pick the statement up as Action<int>.

public void act(Func<int,int> assign) {
    for(int n = 0; n < board.Length; n++)
        act((n) => { board[n] = assign(n); });
}

Having made the change, I tested my input string and out printed the solution immediately. I decided against writing similar function overloading in the final implementation to prevent any unintended bugs.

Written by lewellen

2009-05-01 at 12:00 am

Posted in Algorithms

Tagged with Lambda Expressions, Sudoku

Integer Factorization by Dynamic Programming with Number Theoretic Applications

with 6 comments

Having been a participant of a number of mathematical programming competitions over the years, I’ve had to find a number of efficient ways of implementing many common Number Theoretic Functions. In this write up, I’m going to go over a method I’ve found useful for easily factoring numbers using a sieving method, go over some the implementation of a few number theory functions along with time complexity analysis of each. The cornerstone to many of these implementations relies on the ability to quickly factor integers and find primes.

One of the most popular methods of for finding primes is the Sieve of Eratosthenes. The algorithm starts by populating a table with every positive integer from 2 to a ceiling value. Then find the first integer not yet crossed off, in the case 2, and eliminate every multiple of 2 from the table then return to 2 and find the next positive integer not yet crossed off and repeat the procedure until the end of the table is reached. The method is fine and all, but a lot of really great information is lost in that computation. Here is a sample implementation:

bool[] isPrime = new bool[400];
for (uint n = 2; n < isPrime.Length; n++)
    isPrime[n] = true;
for (uint n = 2; n < isPrime.Length; n++)
    if (isPrime[n])
        for (uint m = 2, c = 0; (c = m * n) < isPrime.Length; m++)
            isPrime[c] = false;

On the other hand, say we approach sieve a little differently. Create an empty table as large as the ceiling value. Start at 2 and for every multiple of 2, create a record that has two parts: 2 and half of the multiple. Return to 2 and find the next integer in the table that has yet to be recorded, in this case 3. For every multiple of 3, create a record that has two parts: 3 and third of the multiple (only if the multiple was not previously recorded. E.g., 6 because 2 previously recorded the record). Return to 3 and find the next integer in the table that has yet to be recorded so on and so forth until every integer in the table has been recorded.

The following graphic demonstrates this process for a ceiling values of 25. If we wish to factor 16, we go to 16′s record (2, 8), follow to 8′s record (2,4), again follow 4′s record (2, 2) and finally 2′s record (2, λ). Thus the prime factorization of 16 is 2, 2, 2, 2.

It should be apparent that this algorithm is a simple dynamic programming solution that yields two major results:

We have factored every positive integer up to a ceiling value.
We have found every positive prime integer up to a ceiling value.

And one major draw back

Uses a lot of memory as a trade off for speed.

Let’s get into the C# implementation. To start off, we need a record class that’ll store the information about the first prime that divides an entry and the composite to jump to if the record corresponds to a composite.

public class Record {
    public uint Prime {get; set;}
    public uint? JumpTo {get; set;}
}

We’ll have a class called NumberTheory and assume that it has the following structure. If you want, you could make this a Singleton class but I felt it was unnecessary for the scope of this write up.

public class NumberTheory {
    private Record[] table;

    ...
}

It makes sense to put the core algorithm in the constructor and then have member methods for each of the functions we’d like to have. It should be assumed that for the lifetime of the class that the largest value ever called on the methods will be N otherwise an exception should be thrown by the methods (omitted here for brevity).

public NumberTheory(uint N){
    uint c = 0;
    table = new Record[N+1];
    for(uint n = 2; n < table.Length; n++) {
        if(table[n] != null)
            continue;
        table[n] = new Record() { Prime = n };
        for(uint m = 2; (c = n * m) < table.Length; m++)
            if(table[c] == null)
                table[c] = new Record() { JumpTo = m, Prime = n };
    }
}

The time complexity of the implementation can be derived using some analysis and by having some knowledge of certain identities. Starting from 2 there are N – 2 numbers to check of which 1/2 will be visited by the interior loop, starting from 3 there are N – 3 numbers to check of which 1/3 will be visited visited by the interior loop, so on and so forth leading to the following summation:
$\displaystyle T(n) = \sum_{p \le n} \left( \frac{1}{p}(n - p) \right)$
$\displaystyle T(n) = \sum_{p \le n} \left( \frac{n}{p} - 1 \right)$

If we separate the summation into the harmonic series of primes (HSP) and prime counting function (aggregate 1 for every prime less than n) we get:
$\displaystyle T(n) = n \sum_{p \le n} \left( \frac{1}{p} \right) - \pi(n)$

Asymptotically, the HSP tends towards $\ln \ln n$ and $\pi(n)$ towards $\frac{n}{\ln n}$ .
$\displaystyle T(n) = n \left( \ln \ln n \right) - \pi(n)$
$\displaystyle T(n) = n \left( \ln \ln n \right) - \frac{n}{\ln n}$
$\displaystyle T(n) = n \frac{ \ln n \cdot \ln \ln(n) - 1}{\ln n}$

Giving us our final asymptotic time complexity of $O(n \ln \ln n)$ .

The first member method we’ll implement is a trivial check to see if a given number is prime by checking the table’s Record’s JumpTo property for null.

public bool IsPrime(uint N){
    return !table[N].JumpTo.HasValue;
}

The time complexity here is a simple $O(1)$ .

From IsPrime, we can easily implement a function that will get every single prime up to a given value by iterating over the table.

public void PrimesLessThan(uint value, Action<uint> actOnPrime) {
    for(int n = 2; n < value; n++)
	if(IsPrime(n))
		actOnPrime(n);
}

We can get the prime counting function $\pi(n)$ but utilizing the PrimeLessThan function.

public uint CountPrimes(uint n) {
	uint count = 0;
	PrimesLessThan(n, (p) => {count++;});
	return count;
}

Here the time complexity is the same as PrimesLessThan: $O(n)$ .

We can get the prime factorization of a composite easily. To do so, we simply iterate over the records until we reach a Record with no JumpTo value.

public void PrimeFactorsOf(uint composite, Action actOnFactor) {
    Record temp = table[composite];
    while(temp != null) {
        actOnFactor(temp.Prime);
        if(temp.JumpTo.HasValue)
            temp = table[temp.JumpTo.Value];
        else
            temp = null;
    }
}

The time complexity of this implementation relies on the the prime omega function $\Omega(n)$ which is the number of prime factors (not necessarily distinct) of n. The function tends to $O( \ln \ln n)$ .

From PrimeFactorsOf, we can also easily implement Euler’s Totient Function $\phi(n)$ - the function tells us how many positive integers less than n are coprime to n. It is defined as:
$\displaystyle \phi(n) = n \prod_{p|n} \left( 1 - \frac{1}{p} \right)$
Which essentially states that if you multiply all of the repeated prime factors of n together by all of the non-repeat prime factors – 1 of n together, you will have the result of $\phi(n)$ .

public uint EulerTotient(uint n){
    uint phi = 1, last = 0;
    PrimeFactorsOf(n, (p) => {
        if(p != last) {
            phi *= p - 1;
            last = p;
        } else {
            phi *= p;
        }
    });
    return phi;
}

A similar function known as the Dedekind $\psi(n)$ Function defined as
$\displaystyle \psi(n) = n \prod_{p|n} \left ( 1 + \frac{1}{p} \right )$
can be implemented in a similar way as $\phi(n)$ :

public uint DedekindPsi(uint n){
    uint phi = 1, last = 0;
    PrimeFactorsOf(n, (p) => {
        if(p != last) {
            phi *= p + 1;
            last = p;
        } else {
            phi *= p;
        }
    });
    return phi;
}

The Von Mangoldt Function $\Lambda(n)$ is another interesting function, unfortunately I haven’t had a chance to use it, but it is trivial to implement so I will include it here for completeness. It is defined as
$\displaystyle \Lambda (n) = \begin{cases} \ln{p} & \text {if n = prime to some positive integer power} \\ 0 & \text{otherwise} \end{cases}$

public double VonMangoldt(uint n){
    uint P = 0;
    PrimeFactorsOf(n, (p) => {
        if(P == 0) {
            P = p;
        } else if (P != p) {
            P = 1;
        }
    });
    return Math.Log(P);
}

The Möbius Function $\mu(n)$ is a handy little function for determining if a number is square free or not (among many more interesting things). It is defined as:
$\displaystyle \mu(n) = \begin{cases} 0 & \text{if n has one or more replicated factors} \\ 1 & \text{if n = 1} \\ (-1)^{k} & \text{if n is a product of k distinct primes} \end{cases}$

public int MoebiusFunction(uint N){
    if(N == 1)
        return 1;
    bool distinct = true;
    uint last = 0, k = 0;
    PrimeFactorsOf(N, (p) => {
        if(p == last) {
            distinct = false;
        } else {
            k++;
            last = p;
        }
    });

    if(distinct) 
        return ((k & 1) == 0) ? 1 : -1; 
    return 0;
}

Since EulerTotient, DedekindPsi, VonMangoldt and Moebius each use PrimeFactorsOf without any additional lifting, their time complexities are the same as PrimeFactorsOf – $O(\ln \ln n)$ .

The last function I’ll implement is the Mertens Function $M(n)$ which is simply defined as $\displaystyle M(n) = \sum_{k = 1}^{n} \mu(k)$ :

public int Mertens(uint n){
    int m = 0;
    for(uint k = 1; k <= n; k++)
        m += MoebiusFunction(k);
    return m;
}

The implementation’s time complexity is simply $O(n \ln \ln n)$ .

Written by lewellen

2009-04-01 at 12:00 am

Posted in Algorithms, Number Theory

Tagged with C-Sharp, Dynamic Programming, Integer Factorization

Converting a Floating-point Number to a Fraction

with 2 comments

The problem

How do we go about converting a floating point number $q$ to a quotient of natural numbers $n, m$ to a given precision $p$ such that $\gcd(n, m) = 1$ ?

There are two ways to go about this, one of which is a subset of the other. We could use the Farey Sequence, but that is limited to mapping floating point values to fractions which fall into the range $\lvert q \rvert \in [0, 1]$ or we can use a Stern-Brocot tree, which allows for values to fall into the range $\vert q \rvert \in [0, \infty]$ .

If we think of the Stern-Brocot tree as a list whose composition is determined by iteration, we will want to define the initial list as the domain $\frac{0}{1}, \frac{1}{0}$ . We want to know a value that lies between these two values such that the numerator and denominator are coprime. (So that the tree we construct contains only fractions that are in minimal terms.) The mediant achieves that goal and is defined as $\frac{a}{c} < \frac{a + b}{c + d} < \frac{b}{d}$ . If we insert the mediant between every pair in the list until the last pair is accessed we’ve completed an iteration. If we take this process ad infinitum we get the Stern-Brocot tree. For example:

$\displaystyle T_{0} = \{ \frac{0}{1}, \frac{1}{0} \}$
$\displaystyle T_{1} = \{ \frac{0}{1}, \frac{1}{1}, \frac{1}{0} \}$
$\displaystyle T_{2} = \{ \frac{0}{1}, \frac{1}{2}, \frac{1}{1}, \frac{2}{1}, \frac{1}{0} \}$
$\displaystyle T_{3} = \{ \frac{0}{1}, \frac{1}{3}, \frac{1}{2}, \frac{2}{3}, \frac{1}{1}, \frac{3}{2}, \frac{2}{1}, \frac{3}{1}, \frac{1}{0} \}$

It is worth noting that if we had instead started with $\frac{0}{1}, \frac{1}{1}$ that we would instead be operating over the Farey Sequence where every iteration is denoted as $\mathfrak{F}_{n}$ . Arguably, we could have gone with this sequence by inverting all $q > 1$ and then inverting the resolved fraction- which would have added unnecessary complexity to the end implementation.

The implementation

Using the Stern-Brocot structure, we can implement a simple binary search over the list to map a floating point number to a fraction within a given precision. The following C# implementation illustrates how this can be done but it is not complete nor is it optimized in any way.

public static Fraction toFraction(this double q, int precision) {
    if (double.IsNaN(q))
        return null;
    if (double.IsPositiveInfinity(q))
        return new Fraction(1, 0);
    else if (double.IsNegativeInfinity(q))
        return new Fraction(-1, 0);
    else if (q == 0.0)
        return new Fraction(0, 0);

    bool negative = q < 0;
    double val;

    Fraction L = new Fraction(0, 1);
    Fraction R = new Fraction(1, 0);
    Fraction M = null;

    q = Math.Round(Math.Abs(q), precision);

    do {
        M = Fraction.mediant(L, R);
        val = M.getValue(precision);

        if (val < q)
            L = M;
        else if(val > q)
            R = M;
    } while (val != q);

    return new Fraction(negative ? -1 : 1 * M.N, M.D);
}

public class Fraction {
       public long N, D;

       public Fraction(long n, long d) {
           N = n; D = d;
       }

       public double value(int precision) {
           return Math.Round(N / (double)D, precision);
       }

       public static Fraction mediant(Fraction l, Fraction r) {
           return new Fraction(l.N + r.N, l.D + r.D);
       }
}

Written by lewellen

2008-07-27 at 9:30 pm

Posted in Algorithms

Tagged with C-Sharp, Stern-Brocot tree

Antimatroid, The

Archive for the ‘Algorithms’ Category

Minesweeper Agent

Introduction

Minesweeper

Background

Specification

Gameplay

Initialization

Exposing Cells

Verification

Methodology

Mine distribution

Distribution of successful clicks

Agents

Methodology

Marginal Model

Development

Test Results

Conditional Model

Development

Test Results

Graphical Model

Development

Test Results

Semi-deterministic Model

Development

Test Results

Summary

References

Tree Drawing: Force-based Algorithm

Sudoku Solver in Haskell

Sudoku Solver in C# using Lambda Expressions

Integer Factorization by Dynamic Programming with Number Theoretic Applications

Converting a Floating-point Number to a Fraction

The problem

The implementation

Pages

Categories

Tags

Archives

Meta

Follow “Antimatroid, The”