Mathematics « Antimatroid, The

Archive for the ‘Mathematics’ 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

Abstract Algebra in C#

leave a comment »

Motivation

In C++ it is easy to define arbitrary template methods for computations involving primitive numeric types because all types inherently have arithmetic operations defined. Thus, a programmer need only implement one method for all numeric types. The compiler will infer the use and substitute the type at compile time and emit the appropriate machine instructions. This is C++’s approach to parametric polymorphism.

With the release of C# 2.0 in the fall of 2005, the language finally got a taste of parametric polymorphism in the form of generics. Unfortunately, types in C# do not inherently have arithmetic operations defined, so methods involving computations must use ad-hoc polymorphism to achieve the same result as in C++. The consequence is a greater bloat in code and an increased maintenance liability.

To get around this design limitation, I’ve decided to leverage C#’s approach to subtype polymorphism and to draw from Abstract Algebra to implement a collection of interfaces allowing for C++ like template functionality in C#. The following is an overview of the mathematical theory used to support and guide the design of my solution. In addition, I will present example problems from mathematics and computer science that can be represented in this solution along with examples how type agnostic computations that can be performed using this solution.

Abstract Algebra

Abstract Algebra is focused on how different algebraic structures behave in the presence of different axioms, operations and sets. In the following three sections, I will go over the fundamental sub-fields and how they are represented under the solution.

In all three sections, I will represent the distinction between algebraic structures using C# interfaces. The type parameters on these interfaces represent the sets being acted upon by each algebraic structure. This convention is consistent with intuitionistic (i.e., Chruch-style) type theory embraced by C#’s Common Type System (CTS). Use of parameter constraints will be used when type parameters are intended to be of a specific type. Functions on the set and elements of the set will be represented by methods and properties respectively.

Group Theory

Group Theory is the simplest of sub-fields of Abstract Algebra dealing with the study of a single binary operation, $(\cdot)$ , acting on a set $a, b, c \in S$ . There are five axioms used to describe the structures studied under Group Theory:

Closure: $(\cdot) : S \times S \to S$
Associativity: $(a \cdot b) \cdot c = a \cdot (b \cdot c)$
Commutativity : $a \cdot b = b \cdot a$
Identity: $a \cdot e = e \cdot a$
Inverse: $a \cdot b = e$

The simplest of these structures is the Groupoid satisfying only axiom (1). Any Groupoid also satisfying axiom (2) is known as a Semi-group. Any Semi-group satisfying axiom (4) is a Monoid. Monoid’s also satisfying axiom (5) are known as Groups. Any Group satisfying axiom (3) is an Abelian Group.

public interface IGroupoid<T> {
    T Operation(T a, T b);
}

public interface ISemigroup<T> : IGroupoid<T> {

}

public interface IMonoid<T> : ISemigroup<T> {
    T Identity { get; }
}

public interface IGroup<T> : IMonoid<T> {
    T Inverse(T t);
}

public interface IAbelianGroup<T> : IGroup<T> {

}

Ring Theory

The next logical sub-field of Abstract Algebra to study is Ring Theory which is the study of two operations, $(\cdot)$ and $(+)$ , on a single set. In addition to the axioms outlined above, there is an addition axiom for describing how one operations distributes over the other.

Distributivity: $a \cdot (b + c) = (a \cdot b) + (a \cdot c), (a + b) \cdot c = (a \cdot c) + (b \cdot c)$

All of the following ring structures satisfy axiom (6). Rings are distinguished by the properties of their operands. The simplest of these structures is the Ringoid where both operands are given by Groupoids. Any Ringoid whose operands are Semi-groups is a Semi-ring. Any Semi-ring whose first operand is a Group is a Ring. Any Ring whose second operand is a Monoid is a Ring with Unity. Any Ring with Unity whose second operand is a Group is Division Ring. Any Division Ring whose operands are both Abelian Groups is a Field.

public interface IRingoid<T, A, M>
    where A : IGroupoid<T>
    where M : IGroupoid<T> {
    A Addition { get; }
    M Multiplication { get; }

    T Distribute(T a, T b);
}

public interface ISemiring<T, A, M> : IRingoid<T, A, M>
    where A : ISemigroup<T>
    where M : ISemigroup<T> {

}

public interface IRing<T, A, M> : ISemiring<T, A, M>
    where A : IGroup<T>
    where M : ISemigroup<T> {

}

public interface IRingWithUnity<T, A, M> : IRing<T, A, M>
    where A : IGroup<T>
    where M : IMonoid<T> {

}

public interface IDivisionRing<T, A, M> : IRingWithUnity<T, A, M>
    where A : IGroup<T>
    where M : IGroup<T> {

}

public interface IField<T, A, M> : IDivisionRing<T, A, M>
    where A : IAbelianGroup<T>
    where M : IAbelianGroup<T> {

}

Module Theory

The last, and more familiar, sub-field of Abstract Algebra is Module Theory which deals with structures with an operation, $(\circ) : S \times R \to R$ , over two separate sets: $a,b \in S$ and $x,y \in R$ that satisfy the following axioms.

Distributivity of $S$ : $a \circ (x + y) = (a \circ x) + (a \circ y)$
Distributivity of $R$ : $(a + b) \circ x = (a \circ x) + (b \circ x)$
Associativity of $S$ : $a \circ (b \circ x) = (a \cdot b) \circ x$

All of the following module structures satisfy axioms (7)-(9). A Module consists of a scalar Ring and an vector Abelian Group. Any Module whose Ring is a Ring with Unity is a Unitary Module. Any Unitary Module whose Ring with Unity is a Abelian Group is a Vector Space.

public interface IModule<
    TScalar, 
    TVector, 
    TScalarRing, 
    TScalarAddativeGroup, 
    TScalarMultiplicativeSemigroup,
    TVectorAddativeAbelianGroup
>
    where TScalarRing : IRing<TScalar, TScalarAddativeGroup, TScalarMultiplicativeSemigroup>
    where TScalarAddativeGroup : IGroup<TScalar>
    where TScalarMultiplicativeSemigroup : ISemigroup<TScalar>
    where TVectorAddativeAbelianGroup : IAbelianGroup<TVector> 
{

    TScalarRing Scalar { get; }
    TVectorAddativeAbelianGroup Vector { get; }

    TVector Distribute(TScalar t, TVector r);
}

public interface IUnitaryModule<
    TScalar, 
    TVector, 
    TScalarRingWithUnity, 
    TScalarAddativeGroup, 
    TScalarMultiplicativeMonoid,
    TVectorAddativeAbelianGroup
> 
    : IModule<
        TScalar, 
        TVector, 
        TScalarRingWithUnity, 
        TScalarAddativeGroup, 
        TScalarMultiplicativeMonoid,
        TVectorAddativeAbelianGroup
    >
    where TScalarRingWithUnity : IRingWithUnity<TScalar, TScalarAddativeGroup, TScalarMultiplicativeMonoid>
    where TScalarAddativeGroup : IGroup<TScalar>
    where TScalarMultiplicativeMonoid : IMonoid<TScalar>
    where TVectorAddativeAbelianGroup : IAbelianGroup<TVector>
{

}

public interface IVectorSpace<
    TScalar,
    TVector,
    TScalarField,
    TScalarAddativeAbelianGroup,
    TScalarMultiplicativeAbelianGroup,
    TVectorAddativeAbelianGroup
>
    : IUnitaryModule<
        TScalar,
        TVector,
        TScalarField,
        TScalarAddativeAbelianGroup,
        TScalarMultiplicativeAbelianGroup,
        TVectorAddativeAbelianGroup
    >
    where TScalarField : IField<TScalar, TScalarAddativeAbelianGroup, TScalarMultiplicativeAbelianGroup>
    where TScalarAddativeAbelianGroup : IAbelianGroup<TScalar>
    where TScalarMultiplicativeAbelianGroup : IAbelianGroup<TScalar>
    where TVectorAddativeAbelianGroup : IAbelianGroup<TVector> 
{

}

Representation of Value Types

The CTS allows for both value and reference types on the .NET Common Language Infrastructure (CLI). The following are examples of how each theory presented above can leverage value types found in the C# language to represent concepts drawn from mathematics.

Enum Value Types and the Dihedral Group $D_8$

One of the simplest finite groups is the Dihedral Group of order eight, $D_{8}$ , representing the different orientations of a square, $e$ , obtained by reflecting the square about the vertical axis, $b$ , and rotating the square by ninety degrees, $a$ . The generating set is given by $\lbrace a, b \rbrace$ and gives rise to the set $\lbrace e, a, a^2, a^3, b, ba, ba^2, ba^3 \rbrace$ These elements are assigned names as follows: $\text{Rot}(0) = e$ , $\text{Rot}(90) = a$ , $\text{Rot}(180) = a^2$ , $\text{Rot}(270) = a^3$ , $\text{Ref}(\text{Ver}) = b$ , $\text{Ref}(\text{Desc}) = ba$ , $\text{Ref}(\text{Hoz}) = ba^2$ and $\text{Ref}(\text{Asc}) = ba^3$ . The relationship between these elements is visualized below:

The easiest way to represent this group as a value type is with an enum.

enum Symmetry { Rot000, Rot090, Rot180, Rot270, RefVer, RefDes, RefHoz, RefAsc }

From this enum we can define the basic Group Theory algebraic structures to take us to $D_8$ .

public class SymmetryGroupoid : IGroupoid<Symmetry> {
    public Symmetry Operation(Symmetry a, Symmetry b) {
        // 64 cases
    }
}

public class SymmetrySemigroup : SymmetryGroupoid, ISemigroup<Symmetry> {

}

public class SymmetryMonoid : SymmetrySemigroup, IMonoid<Symmetry> {
    public Symmetry Identity {
        get { return Symmetry.Rot000; }
    }
}

public class SymmetryGroup : SymmetryMonoid, IGroup<Symmetry> {
    public Symmetry Inverse(Symmetry a) {
        switch (a) { 
            case Symmetry.Rot000:
                return Symmetry.Rot000;
            case Symmetry.Rot090:
                return Symmetry.Rot270;
            case Symmetry.Rot180:
                return Symmetry.Rot270;
            case Symmetry.Rot270:
                return Symmetry.Rot090;
            case Symmetry.RefVer:
                return Symmetry.RefVer;
            case Symmetry.RefDes:
                return Symmetry.RefAsc;
            case Symmetry.RefHoz:
                return Symmetry.RefHoz;
            case Symmetry.RefAsc:
                return Symmetry.RefDes;
        }

        throw new NotImplementedException();
    }

}

Integral Value Types and the Commutative Ring with Unity over $\mathbb{Z} / 2^n \mathbb{Z}$

C# exposes a number of fixed bit integral value types that allow a programmer to pick an integral value type suitable for the scenario at hand. Operations over these integral value types form a commutative ring with unity whose set is the congruence class $\mathbb{Z} / 2^n \mathbb{Z} = \lbrace \overline{0}, \overline{1}, \ldots, \overline{2^n-1} \rbrace$ where $n$ is the number of bits used to represent the integer and $\overline{m}$ is the equivalance class $\overline{m} = \lbrace m + k \cdot 2^n\rbrace$ with $k \in \mathbb{Z}$ .

Addition is given by $\overline{a} + \overline{b} = \overline{(a + b)}$ just as multiplication is given by $\overline{a} \cdot \overline{b} = \overline{(a \cdot b)}$ . Both statements are equivalent to the following congruence statements: $(a + b) \equiv c \pmod{2^n}$ and $(a \cdot b) \equiv c \pmod{2^n}$ respectively.

Under the binary numeral system, modulo $2^n$ is equivalent to ignoring the bits exceeding $n-1$ , or equivalently, $\displaystyle \sum_{i = 0}^{\infty} c_i 2^i \equiv \sum_{i = 0}^{n-1} c_i 2^i \pmod{2^n}$ where $c \in \lbrace 0, 1 \rbrace$ . As a result only the first (right most) bits need to be considered when computing the sum or product of two congruence classes, or in this case, integer values in C#. Thus, in the following implementation, it is not necessary to write any extra code to represent these operations other than writing them in their native form.

The reason why we are limited to a commutative ring with unity instead of a full field is that multiplicative inverses do not exist for all elements. A multiplicative inverse only exists when $ax \equiv 1 \pmod{2^n}$ where $x$ is the multiplicative inverse of $a$ . For a solution to exist, $\gcd(a, 2^n) = 1$ . Immediately, any even value of $a$ will not have a multiplicative inverse in $\mathbb{Z} \ 2^n \mathbb{Z}$ . However, all odd numbers will.

public class AddativeIntegerGroupoid : IGroupoid<long> {
    public long Operation(long a, long b) {
        return a + b;
    }
}

public class AddativeIntegerSemigroup : AddativeIntegerGroupoid, ISemigroup<long> {

}

public class AddativeIntegerMonoid : AddativeIntegerSemigroup, IMonoid<long> {
    public long Identity {
        get { return 0L; }
    }
}

public class AddativeIntegerGroup : AddativeIntegerMonoid, IGroup<long> {
    public long Inverse(long a) {
        return -a;
    }
}

public class AddativeIntegerAbelianGroup : AddativeIntegerGroup, IAbelianGroup<long> {

}

public class MultiplicativeIntegerGroupoid : IGroupoid<long> {
    public long Operation(long a, long b) {
        return a * b;
    }
}

public class MultiplicativeIntegerSemigroup : MultiplicativeIntegerGroupoid, ISemigroup<long> {

}

public class MultiplicativeIntegerMonoid : MultiplicativeIntegerSemigroup, IMonoid<long> {
    public long Identity {
        get { return 1L; }
    }
}

public class IntegerRingoid : IRingoid<long, AddativeIntegerGroupoid, MultiplicativeIntegerGroupoid> {
    public AddativeIntegerGroupoid Addition { get; private set;}
    public MultiplicativeIntegerGroupoid Multiplication { get; private set;}

    public IntegerRingoid() {
        Addition = new AddativeIntegerGroupoid();
        Multiplication = new MultiplicativeIntegerGroupoid();
    }

    public long Distribute(long a, long b) {
        return Multiplication.Operation(a, b);
    }
}

public class IntegerSemiring : IntegerRingoid, ISemiring<long, AddativeIntegerSemiring, MultiplicativeIntegerSemiring> {
    public AddativeIntegerSemiring Addition { get; private set;}
    public MultiplicativeIntegerSemiring Multiplication { get; private set;}

    public IntegerSemiring() : base() {
        Addition = new AddativeIntegerSemiring();
        Multiplication = new MultiplicativeIntegerSemiring();
    }
}

public class IntegerRing : IntegerSemiring, IRing<long, AddativeIntegerGroup, MultiplicativeIntegerSemigroup>{
    public new AddativeIntegerGroup Addition { get; private set; }

    public IntegerRing() : base() {
        Addition = new AddativeIntegerGroup();
    }
}

public class IntegerRingWithUnity : IntegerRing, IRingWithUnity<long, AddativeIntegerGroup, MultiplicativeIntegerMonoid> {
    public MultiplicativeIntegerMonoid Multiplication { get; private set; }

    public IntegerRingWithUnity() : base() {
        Multiplication = new MultiplicativeIntegerMonoid();
    }
}

Floating-point Value Types and the Real Vector Space $\mathbb{R}^n$

C# offers three types that approximate the set of Reals: floats, doubles and decimals. Floats are the least representative followed by doubles and decimals. These types are obviously not continuous, but the error involved in rounding calculations with respect to the calculations in question are negligible and for most intensive purposes can be treated as continuous.

As in the previous discussion on the integers, additive and multiplicative classes are defined over the algebraic structures defined in the Group and Ring Theory sections presented above. In addition to these implementations, an additional class is defined to describe a vector.

public class Vector<T> {
    private T[] vector;

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

    public T this[int n] {
        get { return vector[n]; }
        set { vector[n] = value; }
    }

    public Vector() {
        vector = new T[2];
    }
}

With these classes, it is now possible to implement the algebraic structures presented in the Module Theory section from above.

public class RealVectorModule : IModule<double, Vector<double>, RealRing, AddativeRealGroup, MultiplicativeRealSemigroup, VectorAbelianGroup<double>> {
    public RealRing Scalar {
        get;
        private set;
    }

    public VectorAbelianGroup<double> Vector {
        get;
        private set;
    }

    public RealVectorModule() {
        Scalar = new RealRing();
        Vector = new VectorAbelianGroup<double>(new AddativeRealAbelianGroup());
    }

    public Vector<double> Distribute(double t, Vector<double> r) {
        Vector<double> c = new Vector<double>();
        for (int i = 0; i < c.Dimension; i++)
            c[i] = Scalar.Multiplication.Operation(t, r[i]);
        return c;
    }
}

public class RealVectorUnitaryModule : RealVectorModule, IUnitaryModule<double, Vector<double>, RealRingWithUnity, AddativeRealGroup, MultiplicativeRealMonoid, VectorAbelianGroup<double>> {
    public new RealRingWithUnity Scalar {
        get;
        private set;
    }

    public RealVectorUnitaryModule()
        : base() {
        Scalar = new RealRingWithUnity();
    }
}

public class RealVectorVectorSpace : RealVectorUnitaryModule, IVectorSpace<double, Vector<double>, RealField, AddativeRealAbelianGroup, MultiplicativeRealAbelianGroup, VectorAbelianGroup<double>> {
    public new RealField Scalar {
        get;
        private set;
    }

    public RealVectorVectorSpace()
        : base() {
        Scalar = new RealField();
    }
}

Representation of Reference Types

The following are examples of how each theory presented above can leverage reference types found in the C# language to represent concepts drawn from computer science.

Strings, Computability and Monoids

Strings are the simplest of reference types in C#. From an algebraic structure point of view, the set of possible strings, $\Sigma^{*}$ , generated by an alphabet, $\Sigma$ , and paired with a concatenation operation, $(+)$ , forms a monoid.

public class StringGroupoid : IGroupoid<string> {
    public string Operation(String a, String b) {
        return string.Format("{0}{1}", a, b);
    }
}

public class StringSemigroup : StringGroupoid, ISemigroup<string> {

}

public class StringMonoid : StringSemigroup, IMonoid<string> {
    public string Identity {
        get { return string.Empty; }
    }
}

Monoids over strings have a volley of applications in the theory of computation. Syntactic Monoids describe the smallest set that recognizes a formal language. Trace Monoids describe concurrent programming by allowing different characters of an alphabet to represent different types of locks and synchronization points, while the remaining characters represent processes.

Classes, Interfaces, Type Theory and Semi-rings

Consider the set of types $\mathcal{T}^{*}$ , that are primitive and constructed in C#. The generating set of $\mathcal{T}^{*}$ is the set of primitive reference and value types, $\mathcal{T}$ , consisting of the types discussed thus far. New types can be defined by defining classes and interfaces.

A simple operation $(\oplus)$ over $\mathcal{T}^{*}$ takes two types, $\alpha, \beta$ , and yields a third type, $\gamma$ , known as a sum type. In type theory, this means that an instance of $\gamma$ can be either an instance of $\alpha$ or $\beta$ . A second operation $(\otimes)$ over $\mathcal{T}^{*}$ takes two types and yields a third type representing a tuple of the first two types. In other words, $\alpha \otimes \beta = (\alpha, \beta)$ .

Both operations form a semi-group $(\mathcal{T}^{*}, \otimes)$ and $(\mathcal{T}^{*}, \otimes)$ and in conjunction the two form a semi-ring.

To implement this semi-ring is a little involved. The .NET library supports emitting dynamic type definitions at runtime. For sum types, this would lead to an inheritance view of the operation. Types $\alpha$ and $\beta$ would end up deriving from $\gamma$ which means that any sequence of sums would yield an inheritance tree. A product type would result in composition of types with projection operations, $\pi_{n} : \prod_{i = 0} \tau_{i} \to \tau_{n}$ , to access and assign the $n\text{'th}$ element of the composite. Both type operation implementations are outside the scope of this write-up and I’ll likely revisit this topic in a future write-up.

Delegates and Process Algebras

The third type of reference type to mention is the delegate type which is C#’s approach to creating first-class functions. The simplest of delegates is the built-in Action delegate which represents a single procedure taking no inputs and returning no value.

Given actions $a, b \in \mathcal{A}$ , we can define a possible execution operator, $(\Vert) : \mathcal{A} \times \mathcal{A} \to \mathcal {A}$ , where either $a$ or $b$ is executed denoted as $a \lVert b$ . The choice operation forms a commutative semigroup $(\mathcal{A}, \lVert)$ since operations are associative $a \lVert (b \lVert c) = (a \lVert b) \lVert c$ and the operation is commutative $a \lVert b = b \lVert a$ .

A product operation, $(\to) : \mathcal{A} \times \mathcal{A} \to \mathcal {A}$ , representing the sequential execution of $a$ and then $b$ is given by $a \to b$ . The sequence operator forms a groupoid with unity since the operation is not associative $a \to (b \to c) \neq (a \to b) \to c$ and there is an identity action $e$ representing a void operation resulting in $e \to a = a$ .

Both operations together form a ringoid, $(\mathcal{A}, \to, \lVert)$ since the sequence operation distributes over the choice operation $a \to (b \lVert c) = (a \to b) \lVert (a \to c)$ . Meaning that $a$ takes place and then $b$ or $c$ takes places is equivalent to $a$ and then $b$ takes place or $a$ and then $c$ takes place.

public class SequenceGroupoidWithUnity<Action> : IGroupoid<Action> {
    public Action Identity { 
        get { return () => {}; }
    }

    public Action Operation(Action a, Action b) {
        return () => { a(); b(); }
    }
}

public class ChoiceGroupoid<Action> : IGroupoid<Action> {
    public Action Operation(Action a, Action b) {
        if(DateTime.Now.Ticks % 2 == 0)
            return a;
        return b;
    }
}

The process algebra an be extended further to describe parallel computations with an additional operation. The operations given thus far enable one to derive the possible execution paths in a process. This enables one to comprehensively test each execution path to achieve complete test coverage.

Examples

The motivation of this work was to achieve C++’s approach to parametric polymorphism by utilizing C# subtype polymorphism to define the algebraic structure required by a method (akin to the built-in operations on types in C++). To illustrate how these interfaces are to be used, the following example extension methods operate over a collection of a given type and accept the minimal algebraic structure to complete the computation. The result is a single implementation of the calculation that one would expect in C++.

static public class GroupExtensions {
    static public T Sum<T>(this IEnumerable<T> E, IMonoid<T> m) {
        return E
            .FoldL(m.Identity, m.Operation);
    }
}

static public class RingoidExtensions {
    static public T Count<T, R, A, M>(this IEnumerable<R> E, IRingWithUnity<T, A, M> r)
        where A : IGroup<T>
        where M : IMonoid<T> {

        return E
            .Map((x) => r.Multiplication.Identity)
            .Sum(r.Addition);
    }

    static public T Mean<T, A, M>(this IEnumerable<T> E, IDivisionRing<T, A, M> r)
        where A : IGroup<T>
        where M : IGroup<T> {

        return r.Multiplication.Operation(
            r.Multiplication.Inverse(
                E.Count(r)
            ),
            E.Sum(r.Addition)
        );
    }

    static public T Variance<T, A, M>(this IEnumerable<T> E, IDivisionRing<T, A, M> r)
        where A : IGroup<T>
        where M : IGroup<T> {

        T average = E.Mean(r);

        return r.Multiplication.Operation(
            r.Multiplication.Inverse(
                E.Count(r)
            ),
            E
                .Map((x) => r.Addition.Operation(x, r.Addition.Inverse(average)))
                .Map((x) => r.Multiplication.Operation(x, x) )
                .Sum(r.Addition)
        );
    }
}

static public class ModuleExtensions {
    static public TV Mean<TS, TV, TSR, TSRA, TSRM, TVA>(this IEnumerable<TV> E, IVectorField<TS, TV, TSR, TSRA, TSRM, TVA> m)
        where TSR : IField<TS, TSRA, TSRM>
        where TSRA : IAbelianGroup<TS>
        where TSRM : IAbelianGroup<TS>
        where TVA : IAbelianGroup<TV> {

        return m.Distribute(
            m.Scalar.Multiplication.Inverse(
                E.Count(m.Scalar)
            ),
            E.FoldL(
                m.Vector.Identity,
                m.Vector.Operation
            )
        );
    }
}

Conclusion

Abstract Algebra comes with a rich history and theory for dealing with different algebraic structures that are easily represented and used in the C# language to perform type agnostic computations. Several examples drawn from mathematics and computer science illustrated how the solution can be used for both value and reference types in C# and be leveraged in the context of a few example type agnostic computations. The main benefit of this approach is that it minimizes the repetitious coding of computations required under the ad-hoc polymorphism approach adopted by the designers of C# language. The downside is that several structures must be defined for the types being computed over and working with C# parameter constraint system can be unwieldy. While an interesting study, this solution would not be practical in a production setting under the current capabilities of the C# language.

References

Baeten, J.C.M. A Brief History of Process Algebra [pdf]. Department of Computer Science, Technische Universiteit Eindhoven. 31 Mar. 2012.

ECMA International. Standard ECMA-335 Common Language Infrastructure [pdf]. 2006.

Fokkink, Wan. Introduction of Process Algebra [pdf]. 2nd ed. Berlin: Springer-Verlang, 2007. 10 Apr. 2007. 31 Mar. 2012.

Goodman, Joshua. Semiring Parsing [pdf]. Computational Linguistics 25 (1999): 573-605. Microsoft Research. 31 Mar. 2012.

Hungerford, Thomas. Algebra. New York: Holt, Rinehart and Winston, 1974.

Ireland, Kenneth. A classical introduction to modern number theory. New York: Springer-Verlag, 1990.

Litvinov, G. I., V. P. Maslov, and A. YA Rodionov. Universal Algorithms, Mathematics of Semirings and Parallel Computations [pdf]. Spring Lectures Notes in Computational Science and Engineering. 7 May 2010. 31 Mar. 2012 .

Mazurkiewicz, Antoni. Introduction to Trace Theory [pdf]. Rep. 19 Nov. 1996. Institute of Computer Science, Polish Academy of Sciences. 31 Mar. 2012.

Pierce, Benjamin. Types and programming languages. Cambridge, Mass: MIT Press, 2002.

Stroustrup, Bjarne. The C++ Programming Language. Reading, Mass: Addison-Wesley, 1986.

Written by lewellen

2012-04-01 at 8:00 am

Posted in Abstract Algebra, Theoretical Computer Science

Tagged with C++, C-Sharp, Group Theory, Module Theory, Ring Theory

Menger Sponge in C++ using OpenGL

with one comment

This past summer I was going through some old projects and came across a Menger Sponge visualizer that I wrote back in college. A Menger Sponge is simple fractal that has infinite surface area and encloses zero volume. The sponge is constructed in successive iterations and the first four iterations are rendered in the video below.

The sponge starts as a single cube that is segmented into twenty-seven equally sized cubes. The center cubes of each face and that of the parent cube are then discarded and the process is applied again to each of the remaining cubes. Visually, the process looks like the following:

The geometry of the process is straight forward. Starting with a cube’s origin, $\vec{o}$ , and edge length, $e$ , each of the children’s attributes can be calculated. Each child’s edge length is given by $e_{Child} = \frac{1}{3} e_{Parent}$ . Each child’s origin given by $\vec{o}_{Child} = \vec{o}_{Parent} + e_{Child} \vec{c}_{Child}$ . The constant represents a child’s relative position (e.g., $(1, -1, 0)$ ) to its parent.

The following implementation isn’t particularly well written, but it accomplishes the desired end result. The point and Cube classes achieve the logic that I’ve outlined above. Cube can be thought of as a tree structure that is generated upon instantiation. The visualize() method pumps out the desired OpenGL commands to produce the faces of the cubes.

#include <GL\glut.h>

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>

//=================================================================================
//=================================================================================
class point
{
public:
	point(GLfloat x, GLfloat y, GLfloat z, point* ref = NULL);
	void visualize();

	GLfloat x,y,z;
};

point::point(GLfloat x, GLfloat y, GLfloat z, point* ref)
{
	this->x = x;
	this->y = y;
	this->z = z;

	if(ref != NULL)
	{
		this->x += ref->x;
		this->y += ref->y;
		this->z += ref->z;
	}
}

//=================================================================================
//=================================================================================

class Cube
{
public:
	Cube(point* origin, GLfloat edgelength, GLfloat depth);
	~Cube();

	void visualize();

private:
	void MakeFace(int i, int j, int k, int l);
	void ActAsContainer(point* o, GLfloat e, GLfloat d);
	void ActAsCube(point* o, GLfloat e);

	point** PointCloud;
	Cube** SubCubes;
};

Cube::Cube(point* origin, GLfloat edgelength, GLfloat depth)
{
	if(depth <= 1.0)
	{
		ActAsCube(origin, edgelength);
	} else {
		ActAsContainer(origin, edgelength, depth);
	}
}

Cube::~Cube()
{
	int i;

	if(PointCloud != NULL)
	{
		for(i = 0; i < 8; i++)
			delete PointCloud[i];
		delete[] PointCloud;
	}

	if(SubCubes != NULL)
	{
		for(i = 0; i < 20; i++)
			delete SubCubes[i];
		delete[] SubCubes;
	}
}

void Cube::ActAsCube(point* o, GLfloat e)
{
	GLfloat ne = e / 2.0;

	PointCloud = new point*[8];		// This is the actual physical cube coordinates;
	SubCubes = NULL;

	PointCloud[0] = new point( ne,  ne,  ne, o);	// net
	PointCloud[1] = new point( ne, -ne,  ne, o);	// set
	PointCloud[2] = new point(-ne,  ne,  ne, o);	// nwt
	PointCloud[3] = new point(-ne, -ne,  ne, o);	// swt
	PointCloud[4] = new point( ne,  ne, -ne, o);	// neb
	PointCloud[5] = new point( ne, -ne, -ne, o);	// seb
	PointCloud[6] = new point(-ne,  ne, -ne, o);	// nwb
	PointCloud[7] = new point(-ne, -ne, -ne, o);	// swb
}

void Cube::ActAsContainer(point* o, GLfloat e, GLfloat d)
{
	GLfloat ne = e / 3.0;

	SubCubes = new Cube*[20];	// These are the centers of each sub cube structure
	PointCloud = NULL;

	SubCubes[0] = new Cube(new point(-ne,  ne,  ne, o), ne, d-1.0);
	SubCubes[1] = new Cube(new point(0.0,  ne,  ne, o), ne, d-1.0);
	SubCubes[2] = new Cube(new point( ne,  ne,  ne, o), ne, d-1.0);
	SubCubes[3] = new Cube(new point( ne, 0.0,  ne, o), ne, d-1.0);
	SubCubes[4] = new Cube(new point( ne, -ne,  ne, o), ne, d-1.0);
	SubCubes[5] = new Cube(new point(0.0, -ne,  ne, o), ne, d-1.0);
	SubCubes[6] = new Cube(new point(-ne, -ne,  ne, o), ne, d-1.0);
	SubCubes[7] = new Cube(new point(-ne, 0.0,  ne, o), ne, d-1.0);
	
	SubCubes[8] = new Cube(new point( ne,  ne,  0.0, o), ne, d-1.0);
	SubCubes[9] = new Cube(new point( ne, -ne,  0.0, o), ne, d-1.0);
	SubCubes[10] = new Cube(new point(-ne, ne,  0.0, o), ne, d-1.0);
	SubCubes[11] = new Cube(new point(-ne, -ne,  0.0, o), ne, d-1.0);
	
	SubCubes[12] = new Cube(new point(-ne,  ne, -ne, o), ne, d-1.0);
	SubCubes[13] = new Cube(new point(0.0,  ne, -ne, o), ne, d-1.0);
	SubCubes[14] = new Cube(new point( ne,  ne, -ne, o), ne, d-1.0);
	SubCubes[15] = new Cube(new point( ne, 0.0, -ne, o), ne, d-1.0);
	SubCubes[16] = new Cube(new point( ne, -ne, -ne, o), ne, d-1.0);
	SubCubes[17] = new Cube(new point(0.0, -ne, -ne, o), ne, d-1.0);
	SubCubes[18] = new Cube(new point(-ne, -ne, -ne, o), ne, d-1.0);
	SubCubes[19] = new Cube(new point(-ne, 0.0, -ne, o), ne, d-1.0);
}

void Cube::MakeFace(int i, int j, int k, int l)
{
		glVertex3f(PointCloud[i]->x, PointCloud[i]->y, PointCloud[i]->z);
		glVertex3f(PointCloud[j]->x, PointCloud[j]->y, PointCloud[j]->z);
		glVertex3f(PointCloud[k]->x, PointCloud[k]->y, PointCloud[k]->z);
		glVertex3f(PointCloud[l]->x, PointCloud[l]->y, PointCloud[l]->z);

}

void Cube::visualize()
{
	int i;

	if(PointCloud != NULL)
	{
		glBegin(GL_QUADS);
			glColor3f(1.0,0.0,0.0);// top
			MakeFace(0,2,3,1);
			glColor3f(0.0,1.0,1.0);//bottom
			MakeFace(4,6,7,5);
			
			glColor3f(0.0,1.0,0.0);// north
			MakeFace(0,2,6,4);
			glColor3f(1.0,0.0,1.0);// south
			MakeFace(1,3,7,5);
			
			glColor3f(0.0,0.0,1.0);//east
			MakeFace(0,4,5,1);
			glColor3f(1.0,1.0,0.0);// west
			MakeFace(2,6,7,3);
		glEnd();
	}

	if(SubCubes != NULL)
	{
		for(i = 0; i < 20; i++)
		{
			SubCubes[i]->visualize();
		}
	}
}

The implementation of the program is your run-of-the-mill OpenGL boilerplate. The application takes in an argument dictating what order of sponge it should produce. It sets up the camera and positions the sponge at the origin. The sponge is left stationary, while the camera is made to orbit upon each display(). On idle(), a redisplay message is sent back to the OpenGL system in order to achieve the effect that the sponge is spinning.

//=================================================================================
//=================================================================================
Cube* MengerCube;

void idle()
{
	glutPostRedisplay();
}

void display()
{
	static GLfloat rtri = 0.0;
	
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
	gluLookAt(1.0,1.0,1.0, 0.0,0.0,0.0,0.0,1.0,0.0);
	glRotatef((rtri+=0.932), 1.0, 0.5, -1.0);

	MengerCube->visualize();

	glutSwapBuffers();
}

void reshape(int w, int h)
{
	glViewport(0,0,w,h);
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	glOrtho(-8.0, 8.0,-8.0, 8.0,-8.0, 8.0);
}

void init()
{
	glShadeModel(GL_SMOOTH);
	glClearColor(0.0, 0.0, 0.0, 0.0);
	glClearDepth(1.0f);
	glEnable(GL_DEPTH_TEST);
	glColor3f(1.0, 1.0, 1.0);
}

GLfloat getDepth(char* depth)
{
	int k = atoi(depth);
	if(k <= 1) return 1.0;
	else if (k>= 5) return 5.0;
	else return (GLfloat) k;
}

int main(int argc, char* argv[])
{
	GLfloat depth;
	bool viewAsPointCloud = false;
	point origin(0.0, 0.0, 0.0);

	printf("%d\n",argc);

	switch(argc)
	{
	case 2:
		depth = getDepth(argv[1]);
		break;
	default:
		depth = 2.0;
		break;
	}

	MengerCube = new Cube(&origin, 8.0, depth);

	glutInit(&argc, argv);
	glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
	glEnable(GL_DEPTH_TEST);
	glutInitWindowSize(500,500);
	glutInitWindowPosition(0,0);
	glutCreateWindow(*argv);
	glutReshapeFunc(reshape);
	glutDisplayFunc(display);
	glutIdleFunc(idle);
	init();
	glutMainLoop();

	delete MengerCube;
}

Written by lewellen

2012-02-01 at 8:00 am

Posted in Computer Graphics, Geometry and Topology

Tagged with C++, Menger Sponge, OpenGL

Water and Wine Problem

leave a comment »

Take two containers, $A$ and $B$ . Initially, $A$ is full of water and $B$ is full of wine. Remove a teaspoon of water from $A$ and put it in $B$ , then a teaspoon of the mixture in $B$ and put it in $A$ . In what proportions are the water and wine now mixed in $A$ and $B$ ? What are the proportions of water and wine in each container after $n$ iterations?

First some assumptions: containers $A$ and $B$ are identical in volume and both capable of containing the combined fluid of each. The initial amount of fluid in each container is identical. When transferring fluids between the containers, none of the fluid is lost. Rather than focusing on what amount a teaspoon represents, we’ll state that we are going to transfer a percentage, $\alpha$ , of $A$ to $B$ . Likewise, when transferring fluid from $B$ to $A$ , we’ll use a percentage $\beta$ . Note that the amount that is transferred, $\alpha \in (0, 1]$ , from $A$ to $B$ has to be equal to the amount that is later transferred, $\beta(\alpha + 1)$ , from $B$ to $A$ . From that fact we can conclude that $\beta = \frac{\alpha}{\alpha + 1}$ . Following these conventions, we arrive at the following formulated events in the problem statement:

	Container A		Container B
	Water	Wine	Water	Wine
Initially, A is full of water and B is full of wine.	$1$	$0$	$0$	$1$
Remove a teaspoon of water from A and put it in B,	$1-\alpha$	$0$	$\alpha$	$1$
then a teaspoon of the mixture in B and put it in A.	$1 - \alpha + \alpha \frac{\alpha}{\alpha+1}$	$\frac{\alpha}{\alpha+1}$	$\alpha(1-\frac{\alpha}{\alpha+1})$	$1-\frac{\alpha}{\alpha+1}$

One way of capturing these events is to model the system, $\mathcal{S}$ , as a matrix. Each row of a matrix represents a container and each column of the matrix represents the amount of water or wine in the container. The initial state of the system, $\mathcal{S}_{0} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}$ , is the amount of water and wine in each container. The act of transferring fluid from one container to another and back again is given by $\mathcal{T} = \begin{pmatrix} 1 & \frac{\alpha}{\alpha+1} \\ 0 & \frac{1}{\alpha+1} \end{pmatrix} \begin{pmatrix} 1-\alpha & 0 \\ \alpha & 1 \end{pmatrix}$ which simplifies to $\frac{1}{\alpha+1} \begin{pmatrix} 1 & \alpha \\ \alpha & 1 \end{pmatrix}$ . To determine the quantities of water and wine in each container after $n$ iterations, we can looks at the recurrence relation $\mathcal{S}_{n} = \mathcal{T} \mathcal{S}_{n-1}$ which simplifies to $\mathcal{S}_{n} = \mathcal{T}^{n}\mathcal{S}_{0}$ .

To figure out a closed form solution to $\mathcal{S}_{n}$ we’ll note that the system is a system of difference equations. The general solution being $f(n) = \sum{ c_{i} \mathcal{Q}_{(i)} \lambda_{i}^n }$ where $c_{i}$ is a coefficient based on initial conditions, $\mathcal{Q}_{(i)}$ is the i’th eigenvector and $\lambda_{i}$ is the i’th eigenvalue, i.e., the i’th solution to $\lvert \mathcal{A} - \lambda \mathcal{I} \rvert = 0$ .

Following this practice, we arrive at $f(n) = \displaystyle c_{1} \begin{pmatrix} 1 \\ 1 \end{pmatrix} 1^{n} + c_{2} \begin{pmatrix} -1 \\ 1 \end{pmatrix} \frac{1-\alpha}{1+\alpha}^n$ . To find $c_{1}$ and $c_{2}$ we solve for $n = 0$ and $\displaystyle c_{1} \begin{pmatrix} 1 \\ 1 \end{pmatrix} + c_{2} \begin{pmatrix} -1 \\ 1 \end{pmatrix} = \begin{pmatrix} 1 \\ 0 \end{pmatrix}$ . We go with $\begin{pmatrix} 1 \\ 0 \end{pmatrix}$ because it represents a container starting with one type of fluid and none of the other. The solution is $\displaystyle f(n) = \frac{1}{2} \begin{pmatrix} 1 \\ 1 \end{pmatrix} - \frac{1}{2} \begin{pmatrix} -1 \\ 1 \end{pmatrix} \frac{1-\alpha}{1+\alpha}^n$ . Important to note that the solution represents both containers and one can swap the water and wine labels.

One last thing of interest is to look at the limiting behavior of $f(n)$ . As $n$ increases, the first term remains constant and the second term tends towards zero given the domain of $\alpha$ . Thus, $\displaystyle \lim_{n \to \infty} f(n) = \frac{1}{2} \begin{pmatrix} 1 \\ 1 \end{pmatrix}$ . Thus, no mater what size scope we take, we’ll always arrive at a container containing equal parts water and wine.

An alternative approach to figuring out the limiting behavior is to look directly at $\mathcal{S}_{n}$ . Matrix computations are time consuming, so we’ll want to have a way to compute $\mathcal{S}_{n}$ that minimizes the number of computations. One way to do this is to apply an Eigen Decomposition of a matrix. The idea is that we can decompose a matrix into a product of three matrices: first is the eigenvectors of the matrix, the second has the matrix’s eigenvalues along the diagonal and the third is an inverse of the first. The decomposition looks like the following: $\mathcal{A} = \mathcal{Q} \Lambda \mathcal{Q}^{-}$ with $\Lambda_{i,i} = \lambda_{i}$ . One important consequence is that we can compute $\mathcal{A}^{n}$ easily now by $\mathcal{Q} \Lambda \mathcal{Q}^{-} \mathcal{Q} \Lambda \mathcal{Q}^{-} \cdots \mathcal{Q} \Lambda \mathcal{Q}^{-} = \mathcal{Q} \Lambda^n \mathcal{Q}^{-}$ and as a direct consequence of how $\Lambda$ is defined, $\Lambda^{n}_{i,i} = \lambda_{i}^{n}$ .

Given that knowledge, we can decompose $\mathcal{T}^{n} = \mathcal{Q}\Lambda^n\mathcal{Q}^{-} = \frac{1}{2} \begin{pmatrix} 1 & -1 \\ 1 & 1 \end{pmatrix} \begin{pmatrix} 1 & 0 \\ 0 & \frac{1-\alpha}{1+\alpha}^{n} \end{pmatrix} \begin{pmatrix} 1 & 1 \\ -1 & 1 \end{pmatrix}$ .

Taking the limit $\displaystyle \lim_{n \rightarrow \infty} \mathcal{T}^n = \frac{1}{2} \begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix}$ leads to the steady state solution of $\mathcal{S}_{\infty} = \frac{1}{2} \begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix}$ .

Written by lewellen

2011-03-01 at 8:00 am

Posted in Dynamical Systems, Linear Algebra

Tagged with Difference Equations, Eigen Decomposition

One tree many possibilities: Part 3

with 2 comments

Introduction

When I first began publishing content on this blog, I wrote a couple of posts entitled One tree many possibilities that covered how to enumerate the sets produced by different counting techniques. You can catch up on the idea in part one and the C# implementation in part two. In this post I’m going to cover how to take the knowledge from the first two posts and derive closed form solutions to each of the counting techniques through analysis of their respective recurrence relations.

Preliminaries

The basic idea is that that the set $S$ and each of its elements, $e$ , are supplied to a binary operation, $\oplus$ , that maps to a subset $S^{\prime}$ of $S$ . This process is performed iteratively on $S^{\prime}$ until the process has been applied to a depth of $k$ or the resulting subset becomes the empty set.

Since we can have a variety of conceivable sets, we need to map sets of a given cardinality to a common set of the same cardinality. Let $b(S) = P_{\lvert S \rvert}$ be a bijective function that maps a set to the set $P_{n}$ of positive integers from 1 to $n$ .

k-ary Cartesian Product

The k-ary cartesian product is the way of selecting $k$ items from a set where the order of the selection is important and items can be placed back into the set after selection.

The process uses $\oplus = id$ where $id(S, e) = S$ is an identity function. As an example, here is a graph of the process applied to $P_{3}$ with $k = 3$ :

To count the number of elements, we need to count the number of nodes at depth $k$ . We will write this using the following recurrence relation:

$\displaystyle f(P_{n}, 0) = 1$
$\displaystyle f(P_{n}, k) = \sum_{i=1}^{n} {f(P_{n}, k - 1) }$

If we think of the tree having a height of $k$ , then we can label each layer from $k$ (at the root) to zero. As a consequence, we can interpret the first statement as a node at depth zero will be considered a leaf node and should be counted once. The second statement says that for each element in $P_{n}$ we will add up whatever number of leaf nodes were counted in the layer beneath the current layer $k - 1$ .

Rearranging the second statement leads to $f(P_{n}, k) = n f(P_{n}, k - 1)$ . Reducing further $k$ times leads to $f(P_{n}, k) = n (n \cdots (n(1))) = n^{k}$ .

k-Permutations

The permutation and k-permutation is a way of selecting $k$ items ( $k \le n$ ) from a list where the order of the selection is important but items cannot be placed back into the set after selection.

The process uses $\oplus = \setminus$ where $\setminus(S, e) = \lbrace s \colon s \ne e, s \in S \rbrace$ is the set difference operator. We use the set difference because the item is removed from the set upon each selection. As before, here is a graph of the process applied to $P_{3}$ with $k=3$ :

To count the number of elements, we need to count the number of nodes at depth $k$ . (The regular permutation can be counted when $k = n$ .) We will write this using the following recurrence relation:

$\displaystyle f(P_{n}, 0) = 1$
$\displaystyle f(P_{n}, k) = \sum_{i = 1}^{n} {f(P_{n-1}, k - 1)}$

Using the same labeling scheme as in the k-ary cartesian product, we can interpret the first statement as any node at depth zero will be considered a leaf node and counted once. The second statement says that for each element in $P_{n}$ we will add up whatever number of leaf nodes were counted in the layer beneath the current layer $k - 1$ for the set $P_{n-1}$ . We use $P_{n-1}$ since we removed an element from $P_{n}$ .

Rearranging the second statement leads to $f(P_{n}, k) = n f(P_{n-1}, k - 1)$ . Taking $k$ times leads to $f(P_{n}, k) = n (n - 1 \cdots (n - k + 1(1))) = \frac{n!}{(n-k)!}$ .

Power Set

The power set is the way of selecting zero items to the cardinality of the set items from a set where the order of the selection is not important and items are not placed back into the set after selection.

The process uses $\oplus = <$ where $<(S, e) = \lbrace s \colon s < e, s \in S \rbrace$ . As before, here is a graph of the process applied to $P_{3}$ :

To count the number of elements, we need to count the number of nodes in the tree. We will write this using the following recurrence relation:

$\displaystyle f(P_{0}) = 1$
$\displaystyle f(P_{n}) = \sum_{i=1}^{n} {f(P_{i-1})}$

The first statement says for any node a depth zero we will count it once. The second statement states that for each element in $P_{n}$ we will add up the result from $P_{0}$ to $P_{n-1}$ . Starting at $i = 1$ is a way of encoding that all nodes should be counted. We go up to $P_{n-1}$ because there are $n - 1$ elements in $P_{n}$ less than $n$ .

Rearranging the second statement leads to $f(P_{n}) = \sum_{i=1}^{n-1} {f(P_{i-1})} + f(P_{n-1}) = 2 f(P_{n-1})$ . Taking the expression to zero leads to $f(P_{n}) = 2(2 \cdots(2(1))) = 2^{n}$ .

k-Combinations

The k-combination is the way of selecting a fixed number of items from a set where the order of the selection is not important and items are not placed back into the set after selection.

The process is identical to the power set.

To count the number of elements, we need to count the number of nodes at depth $k$ in the tree. We will write this using the following recurrence relation:

$\displaystyle f(P_{n}, 0) = 1$
$\displaystyle f(P_{n}, k) = \sum_{i=1}^{n} {f(P_{i-1}, k-1) }$

The first statement says for depth zero we will count once, the second statement states that we will add up whatever sum was produced by iterating over $P_{i-1}$ from one to $n$ and reducing the depth $k$ by one. To find the closed form solution, we’ll approach the recurrence relation inductively:

$\displaystyle f(P_{n}, 0) = 1$
$\displaystyle f(P_{n}, 1) = \sum_{i=1}^{n} {f(P_{i-1}, 0) } = \sum_{i=1}^{n} {1} = n$
$\displaystyle f(P_{n}, 2) = \sum_{i=1}^{n} {f(P_{i-1}, 1) } = \sum_{i=1}^{n} {i-1} = \frac{n(n-1)}{2}$
$\displaystyle f(P_{n}, 3) = \sum_{i=1}^{n} {f(P_{i-1}, 2) } = \sum_{i=1}^{n} {\frac{(i-1)(i-2)}{2}} = \frac{n(n-1)(n-2)}{6}$
$\displaystyle f(P_{n}, k) = \sum_{i=1}^{n} {f(P_{i-1}, k-1) } = \frac{n(n-1)(n-2)\cdots(n-k+1)}{k(k-1)(k-2)\cdots1} = \frac{n!}{k!(n-k)!}$

Important to notice that for $k = 0$ that we have the unit, for $k = 1$ we have the natural numbers, for $k = 2$ we have the triangle numbers, for $k = 3$ we have the tetrahedral numbers and so on. The generalized sequence is the Figurate number- these numbers constitute Pascal’s Triangle which is one method of calculating combinations.

Written by lewellen

2011-02-01 at 8:00 am

Posted in Combinatorics

Tagged with k-ary Cartesian product, k-Combinations, k-Permutations, Permutations, Power Set, Proofs

Mathematics Test Practice Questions

leave a comment »

I’ve had a fairly busy month, so I haven’t had a lot of time to work on any projects. I’ve got a few in the pipeline, but none really mature enough yet to get a write-up. I’ve been playing with the idea of graduate school for a couple years now and have always felt that a strong mathematics background is a good basis for pursuing more interesting (computer science) topics. So, this past week I’ve been working my way through the Mathematics Test Practice Book to see how well I’d fare. The examination seems fairly aggressive allowing for roughly 2.5 minutes per problem for 66 problems. I’m a little more methodical in my ways, so I would certainly take more time than what’s allotted. Nonetheless, I thought I’d share my take on a handful of problems from different subject matters that I found of particular interest:

Let $P_{1}$ be the set of all primes, $\{2, 3, 5, 7,\ldots \}$ , and for each integer $n$ , let $P_{n}$ be the set of all prime multiples of n, $\{2n, 3n, 5n, 7n, \ldots \}$ . Which of the following intersections is nonempty?

$P_{1} \cap P_{23}$

$P_{7} \cap P_{21}$

$P_{12} \cap P_{20}$

$P_{20} \cap P_{24}$

$P_{5} \cap P_{25}$

Given $P_{n} \cap P_{m}$ , we’ll pick two elements (one from each set) and solve for $n p_{1}^{(i)} = m p_{1}^{(j)}$ . We then assume that $n, m$ are of the form $a p_{1}^{(k)}$ where $a$ is a composite. Thus, $a p_{1}^{(j)} p_{1}^{(i)} = a p_{1}^{(i)} p_{1}^{(j)}$ , we can conclude that $p_{1}^{(j)}p_{1}^{(i)} \in P_{n} \cap P_{m}$ . Thus, D is the correct answer.

For what value of $b$ is the line $g(x) = 10x$ tangent to the curve $h(x) = e^{bx}$ at some point in the xy-plane?

Using operator notation, a line tangent to a curve is given by $T_{x_{0}} f = \frac{df}{dx}(x_{0})(x - x_{0}) + f(x_{0})$ . To solve for $b$ we solve for $T_{x_{0}} h = g \implies 10x = be^{b x_0}(x - x_0) + e^{b x_0} \implies 10x = be^{b x_0}x + e^{b x_0}(1 - bx_0)$

Which yields the following two equations: $10 = be^{b x_0}$ and $0 = e^{b x_0}(1 - bx_0)$ since $e^{b x_0}$ will never equal zero, we can conclude $0 = 1 - b x_0 \implies 1 = bx_0 \implies x_0 = \frac{1}{b}$ . Substituting for $x_{0}$ yields $10 = b e^{b \frac{1}{b}} \implies 10 = be \implies b = \frac{10}{e}$ .

Suppose that two binary operations, denoted by $\oplus$ and $\odot$ , are defined on a nonempty set $S$ , and that the following conditions are satisfied for all $x$ , $y$ and $z$ in $S$ :

(1) (Closure) $x \oplus y, x \odot y \in S$
(2) (Associativity) $x \oplus (y \oplus z) = (x \oplus y) \oplus z$ and $x \odot (y \odot z) = (x \odot y) \odot z$
(3) (Commutativity) $x \oplus y = y \oplus x$

Also, for each $x \in S$ and for each positive integer $n$ , the elements $nx$ and $x^n$ are defined recursively as follows:

(4) (Identity Element) $1x = x^1 = x$
(5) and if $kx$ and $x^k$ have been defined, then $(k + 1)x = kx \oplus x$ and $x^{k + 1} = x^k \odot x$

Which of the following must be true?

I only

II only

III only

II and III only

I, II and III

I.) $(x \odot y)^n = x^n \odot y^n$ s.t. $x, y \in S, n \in \mathbb{Z}^{+}$

Assume $x \odot y \in S$ , from which we assume $(x \odot y) \odot (x \odot y) \in S$ . If (3) were defined for $\odot$ , we could swap either side then apply (2) reducing it down to the recursive definition (5) and have a proof. Since that isn’t the case, we can’t conclude that $(x \odot y)^n = x^n \odot y^n$ .

II.) $n (x \oplus y) = nx \oplus ny$ s.t. $x,y \in S, n \in \mathbb{Z}^{+}$

Given (1) we can construct $(x \oplus y) = 1 (x \oplus y)$ by using (4). We can state $2 (x \oplus y) = 1 (x \oplus y) \oplus (x \oplus y)$ from (5), leading to $2 (x \oplus y) = (x \oplus y) \oplus (y \oplus x)$ by way of (3). $2 (x \oplus y) = x \oplus (y \oplus (y \oplus x))$ and $2 (x \oplus y) = x \oplus ((y \oplus y) \oplus x)$ by (2). From (5) $2 (x \oplus y) = x \oplus (2y \oplus x)$ . Unfolding that set of steps leads to the following: $2 (x \oplus y) = x \oplus (x \oplus 2y)$ , $2 (x \oplus y) = (x \oplus x) \oplus 2y$ and $2 (x \oplus y) = 2x \oplus 2y$ . If we were to continue down this path for $n = \{3, 4, 5,\dotsc\}$ then we would find that $n (x \oplus y) = nx \oplus ny$ .

III.) $x^m \odot x^n = x^{m+n}$ s.t. $x \in S, m,n \in \mathbb{Z}^{+}$ .

Starting with (5): $x^{m + 1} = x^m \odot x$ we introduce an additional $x$ , $x^{m + 1} \odot x = (x^m \odot x) \odot x$ and then apply (2) $x^{m + 2}= x^m \odot (x \odot x)$ which leads to $x^{m + 2} = x^m \odot x^2$ . The same set of steps once again leads to $x^{m + 2} \odot x = (x^m \odot x^2) \odot x$ $x^{m + 3} = x^m \odot (x^2 \odot x)$ $x^{m + 3} = x^m \odot x^3$ . By induction, we can conclude $x^{m + n} = x^m \odot x^n$ .

At a banquet, 9 women and 6 men are to be seated in a row of 15 chairs. If the entire seating arrangement is to be chosen at random, what is the probability that all of the men will be seated next to each other in 6 consecutive positions?

I’m going to use the notation $y_{M}$ to denote a number of men ( $Y_{M}$ for the total number of men), $x_{W}$ to denote a number of women ( $X_{W}$ for the total) and $S$ to denote the number of seats. A sequence of this notation such as $3_{W} 6_{M} 6_{W}$ indicates a sequence of three women, six men and six women seated amongst the fifteen seats. Given the definition, we are interested in the probability of $P[ \bigcup_{i = 0}^{X_{W}} {i_{W}Y_{M}(X_{W} - i)_{W}}]$ . Since there are a fixed number of seats and no two people can occupy the same seat, we know that the number of permutations of $S = S!$ is the total number of possible arrangements (our sample space). There are $Y_{M}!$ ways in which the men can be arranged and there are $X_{W}!$ ways in which the women can be arranged in the remaining seats. Since there are $X_{W} + 1$ ways to “slide” the group of men along the seats we can conclude that:

$\displaystyle P[ \bigcup_{i = 0}^{X_{W}} {i_{W}Y_{M}(X_{W} - i)_{W}}] = \frac{(X_{W}+1) X_{W}! Y_{M}!}{S!}$

Given the question, the solution is

$\displaystyle P[ 0_{W}6_{M}9_W \cup 1_{W}6_{M}8_W \cup 2_{W}6_{M}7_W \cdot 9_{W}6_{M}0_W] = \frac{10 * 9! 6!}{15!} = \frac{10! 6!}{15!}$

If $\displaystyle z = e^{\frac{2 \pi i}{5}}$ , then $w = 1 + z + z^2 + z^3 + 5z^4 + 4z^5 + 4z^6 + 4z^7 + 4z^8 + 5z^9 =$ ?

$0$

$4 e^{\frac{3 \pi i}{5}}$

$5 e^{\frac{4 \pi i}{5}}$

$-4 e^{\frac{2 \pi i}{5}}$

$-5 e^{\frac{3 \pi i}{5}}$

First thing to observe is that $z$ is the 5′th root of unity which carries with it two properties:

$z^5 = 1$
$\displaystyle \sum_{k = 0}^{4} z^k = 0$

Proof of property (1) is trivial, proof of property (2) is basic proof of the geometric series:

$\displaystyle \sum_{k = 0}^{n} z^k = a \implies \sum_{k = 1}^{n+1} z^k = z a \implies\displaystyle \sum_{k = 1}^{n+1}{z^k} - \sum_{k = 0}^{n}{z^k} = z a - a \implies$
$\displaystyle z^{n+1} - 1 = a (z - 1) \implies a = \frac{z^{n+1} - 1}{z - 1}$
$\displaystyle \therefore \sum_{k = 0}^{n} z^k = \frac{z^{n+1} - 1}{z - 1}$

Given these properties, we begin to evaluate $w$ by separating the terms into reducible portions:

$\displaystyle w = (1 + z + z^2 + z^3 + z^4) + 4z^4(1 + z + z^2 + z^3 + z^4) + 5(z^4 z^5) \implies$
$\displaystyle w = (1 + 4z^4)\sum_{k = 0}^{4}{z^k} + 5(z^4 1) \implies w = 5z^4 \implies w = 5 e^{\frac{8 \pi i}{5}} \implies$
$\displaystyle w = 5 e^{\pi i} e^{\frac{3 \pi i}{5}} \implies w = -5 e^{\frac{3 \pi i}{5}}$

If $\lfloor{x}\rfloor$ denotes the greatest integer not exceeding x, then $\int_{0}^{\infty}{\lfloor{x}\rfloor e^{-x} dx} =$

$\frac{e}{e^2 - 1}$

$\frac{1}{e - 1}$

$\frac{e - 1}{e}$

$1$

$+\infty$

When $x \in \mathbb{Z}$ , we have a discontinuity in the function that we are integrating because of the floor function. To resolve this issue, we need to break up the integral into a sum of integrals as follows: $\displaystyle \int_{0}^{\infty}{\lfloor{x}\rfloor e^{-x} dx} = \sum_{n = 0}^{\infty} n \int_{n}^{n+1}{e^{-x} dx} = \sum_{n = 0}^{\infty} n (-e^{-x})\rvert_{n}^{n+1} =$ $\displaystyle \sum_{n = 0}^{\infty} n (-e^{-(n+1)} + e^{-n}) = \sum_{n = 0}^{\infty} n \left( \frac{1}{e^n} - \frac{1}{e e^n} \right) = \left( 1 - \frac{1}{e} \right) \sum_{n = 0}^{\infty} n \frac{1}{e}^n$ .

Now we are left with a series that looks similar to a geometric series, so we evaluate it using a similar proof, then return to evaluating the integral.

$\displaystyle \sum_{n = 0}^{m} n x^n = a \implies x \sum_{n = 0}^{m} n x^n = x a \implies \sum_{n = 0}^{m} n x^{n+1} = x a \implies$
$\displaystyle \sum_{n = 0}^{m} n x^{n+1} - \sum_{n = 0}^{m} n x^n = (x - 1) a \implies - \sum_{n = 1}^{m}{x^n} + m x^{m+1} = a(x-1)$ .

(This last step can be produced by expanding the series.) We see that we have a geometric series (without the first term) which was proved earlier. From these facts we can conclude:

$\displaystyle a = \frac{- (\frac{x^n - 1}{x - 1} - 1) + m x^{m+1}}{x-1}$ .

As we take the limit of $m \to \infty$ of this expression (assuming $x < 1$ for convergence) we find:

$\displaystyle a = \frac{- (\frac{- 1}{x - 1} - 1)}{x-1} \implies a = \frac{\frac{1}{x - 1}+1}{x-1} \implies a = \frac{\frac{x}{x - 1}}{x-1} \implies a = \frac{x}{(x-1)^2}$
$\displaystyle \therefore \sum_{n = 0}^{\infty} n x^n = \frac{x}{(x-1)^2}$ .

Returning to evaluating the integral, we now find:

$\displaystyle \int_{0}^{\infty}{\lfloor{x}\rfloor e^{-x} dx} = \frac{e - 1}{e} \frac{\frac{1}{e}}{(\frac{1}{e}-1)^2} = \frac{e - 1}{e} \frac{e}{(1-e)^2} = \frac{e - 1}{(-1)^2(e-1)^2} = \frac{1}{e-1}$ .

The four shaded circles in Figure 1 above are congruent and each is tangent to the large circle and to two of the other shaded circles. Figure 2 is the result of replacing each of the shaded circles in Figure 1 by a figure that is geometrically similar to Figure 1. What is the ratio of the area of the shaded portion of Figure 2 to the area of the shaded portion of Figure 1?

$\frac{1}{2 \sqrt{2}}$

$\frac{1}{1 + \sqrt{2}}$

$\frac{4}{1 + \sqrt{2}}$

$\left( \frac{\sqrt{2}}{1+\sqrt{2}} \right)^2$
$\left( \frac{2}{1+\sqrt{2}} \right)^2$

The first goal is to construct one of the shaded circles in Figure 1. Once the circle has been constructed, we know that the scaling factor of Figure 1 onto Figure 2 will be the radius of the shaded circle (assuming the bounding circle is the unit circle).

Starting with the unit circle, $U$ , centered at the origin, $A$ select points $B = (1, 0)$ and $C = (0, 1)$ on $U$ . Construct the midpoint $D = (\frac{1}{2}, \frac{1}{2})$ to construct the line $\overleftarrow{A}\overrightarrow{D} = x$ . The intersection of $\overleftarrow{A}\overrightarrow{D}$ with $U$ leads to the point $E = (\frac{\sqrt{2}}{2},\frac{\sqrt{2}}{2})$ . Let $\overleftarrow{F}\overrightarrow{G} = -x + \sqrt{2}$ be perpendicular to $\overleftarrow{A}\overrightarrow{D}$ at $E$ . Points $F = (0, \sqrt{2})$ and $G (\sqrt{2}, 0)$ . We can now construct two angle bisectors: $\overleftarrow{F}\overrightarrow{J} = -(1+\sqrt{2})x + \sqrt{2}$ and $\overleftarrow{I}\overrightarrow{G} = (1 - \sqrt{2})x + 2 - \sqrt{2}$ . The intersection of these two lines is $H = ( \sqrt{2} - 1, \sqrt{2} - 1)$ . Which implies the inscribed circle $\mathcal{I} = (x - \sqrt{2} + 1)^2 + (y - \sqrt{2} + 1)^2 = (\sqrt{2} - 1)^2$ .

Thus, the area of $\mathcal{I}$ is $\pi (\sqrt{2} - 1)^2$ and the total shaded area in Figure 1 is $4 \pi (\sqrt{2} - 1)^2$ . To figure out the shaded area of Figure 2, we’ll scale Figure 1 down by a factor of $\sqrt{2} - 1$ and multiply by 4 yielding $16 \pi (\sqrt{2} - 1)^4$ . These two areas give a ratio of: $\displaystyle \frac{16 \pi (\sqrt{2} - 1)^4}{4 \pi (\sqrt{2} - 1)^2} = \frac{4 (\sqrt{2} - 1)^2}{1} \frac{(\sqrt{2}+1)^2}{(\sqrt{2}+1)^2} = \frac{(2(\sqrt{2} - 1)(\sqrt{2}+1))^2}{(\sqrt{2}+1)^2} = \left( \frac{2}{\sqrt{2}+1} \right)^2$ .

(Images in this section were generated by GeoGebra.)

Written by lewellen

2010-07-01 at 8:00 am

Posted in Mathematics

Tagged with GRE

Multi-agent ordering

leave a comment »

Given a set of agents $A = \{ a_{0}, a_{1}, \dotsc, a_{n} \}$ and a set of entities $V = \{ v_{0}, v_{1}, \dotsc, v_{m} \}$ , agents independently provide preferences $I_a \colon V^{2} \to \{ -1, 0, +1\}$ between different pairs of entities for a subset $S_{a}$ of the possible $E = V^{2}$ pairs. From these preferences, we wish to form a poset that contains the entities in order of least preferred to most preferred.

An agent may prefer an individual pair $e_{i,j}$ in the following ways: $-1$ , the agent prefers $v_{j}$ to $v_{i}$ ; $0$ , the agent is indifferent; $+1$ , the agent prefers $v_{i}$ to $v_{j}$ . Let $G_{a} = (V, S_{a})$ be the undirected graph that represents the preferences of an individual agent.

To form a consensus we are interested in $G= \sum_{a \in A} G_{a}$ . Let the adjacency matrix $M$ represent $G$ . In doing so, we may define $M$ in the following way: $M_{i,j} = \sum_{a \in A} I_a (e_{i,j})$ . It should be evident that $M_{i,i} = 0$ and $M_{i,j} = -M_{j,i}$ . We wish to simplify $M$ by turning it into a directed graph whereby all negative entries are removed. We shall define $N_{i,j} = \max{(0, M_{i,j})}$ to be this directed graph.

Let $deg_{j}^{in} = \sum_{i = 0}^{m} N_{i, j}$ be the sum of weights into vertex $v_{j}$ , $deg_{i}^{out} = \sum_{j = 0}^{m} N_{i,j}$ be the sum of weights out of vertex $v_{i}$ and $\delta_i = (deg_{i}^{in}, deg_{i}^{out})$ be the set of tuples for each vertex in $V$ . Thus, let $\delta$ along with the relation: increasing first element, then by decreasing second element be the poset $P$ of entities order from least preferred to most preferred.

Consider the following example: on a merchant’s website, customers $A = \{a,b,c,d,e,f\}$ independently provide comparisons between products that the merchant sells. The merchant wishes to know his customers’ preferences about his products $V = \{1,2,3,4,5,6\}$ .

Customers have provided the following set of preferences:

multi_agent_preferences

The consensus graph is produced from these preferences, followed by the directed graph until $\delta = \begin{pmatrix} (1,4) & (4,3) & (2,0) & (1,4) & (0,2) & (4,0) \end{pmatrix}$ is computed, which, when used to form $P$ , will yield $P = \{6,2,3,1,4,5\}$ . The merchant can infer from this analysis that customers least preferred product is 6, and most preferred product is 5.

Now, this model assumes that the agents are being truthful in their preferences. It is conceivable that a consortium of agents could force an entity to become the least or most preferred entity by providing false preferences. This premise presents two questions: how do we identity fallacious behavior, and two: how do we prevent it. The answer to the later of these two questions is simple: we add an additional term $\tau_a \in [0,1]$ when constructing $M_{i,j} = \sum_{a \in A} (1 - \tau_a) I_a (e_{i,j})$ where $0$ indicates that the agent provides factual preferences and $1$ indicates non factual preferences. The prior of these two questions is more involved and ties into calculating $\tau_a$ .

Lets consider two plausible hypotheses: a single agent tells the truth $\phi_{a}$ or he does not tell the truth about a single comparison $e$ . An agent is telling the truth about a single comparison if his preference $I_{a}(e)$ is identical to the majority preference $\sum_{a \in A} I_{a}(e)$ . An agent is not telling the truth if his preference is not identical to the majority preference. Let $\phi_{a}(e) = 1$ when $I_{a}(e)\sum_{\alpha \in A} I_{\alpha}(e) < 0$ , $0$ otherwise formalize this membership. Thus, we may define $\displaystyle \tau_{a} = \frac{\sum_{e \in S_{a}} \phi_{a}(e)}{ \lvert S_{a} \rvert }$ as the probability that a’s comparisons are majority comparisons.

This formalization assumes that there is an objective comparison and naturally, this formalization would not hold for subjective comparisons.

Consider the following example: seven robots are asked to order a collection of five blocks from smallest to largest and the result is an incorrectly ordered list. We wish to identify the defective robot(s), remove their preferences and reproduce the corrected list. The following table summarizes each robot’s preferences starting from the top left to right then down to the bottom left to right.

Going through the process above we arrive at an order of $\{0, 1, 3, 2, 4\}$ which is incorrect, thus we look at what our values are for this first iteration: $\tau_{a}^{(0)} = \begin{pmatrix}0 & 0.25 & 0 & 0 & 0 & 0.\overline{1} & 0\end{pmatrix}$ . Robots $a_{1}, a_{5}$ are more likely than their piers to deviate from the established majority. So, applying these values, we continue to perform this process iteratively until we arrive at $\{0,1,2,3,4\}$ using values of $\tau_{a}^{(2)} = \begin{pmatrix}0 & 0.41\overline{6} & 0 & 0 & 0 & 0.\overline{2} & 0 \end{pmatrix}$ . This process can be visualized the iterations of $N$ below. Edge weights represent the consensus. Thicker lines for the majority, thinner for the minority.

While not previously alluded to, it is possible that it may take several iterations to achieve a correct ordering. Which raises a couple questions: under what conditions will the correct ordering exist, how many iterations will it take to reach the correct ordering and are the number of iterations finite. The last one is left to the reader.

As a thought experiment, it is evident that there has to be at least $m - 1$ preferences for a correct ordering to be produced. Starting at any one vertex, there exists a path to every other vertex, otherwise a consensus cannot be achieved because $G$ could be decomposed into a set of disjoint subgraphs and no ordering is possible.

If there is an objective order $P^{\star}$ , it is still possible to arrive at a separate order $P^{\prime}$ that is produced after successive iterations if the majority preferences are all incorrect with respect to $P^{\star}$ . Thus, we cannot guarantee that the process will converge. If the preferences are correct with respect to $P^{\star}$ , then convergence will be achieved as the weight of dissenting agent’s preferences tends towards zero.

Written by lewellen

2010-01-01 at 8:00 am

Posted in Graph Theory

Tagged with Topological sorting

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

Generalized Interpreter Framework in C#

leave a comment »

Introduction

Like most Computer Science graduates, I had taken a Programming Languages course during college in which we built a compiler and studied the consequences of our design decisions. I had also taken a course in symbolic logic- particularly, looking at the sentential and predicate logic. Lately, I had been thinking about these two courses and thought it would be interesting to write an interpreter that would take in a sentence valid under sentential logic and output the resulting truth table.

After looking at what I had built, it became apparent that I could generalize my solution into a flexible framework where I could easily create a new interpreter for a desired (albeit simple) formal grammar and produce the desired outputs that other applications could then use for whatever purpose. The following is an overview of this framework and an example implementation for generating truth tables for sentences valid under sentential logic.

Framework

Specification

An Interpreter takes a string and produces an object hierarchy that a client can perform functions upon. The Interpreter does so by following the rules of a grammar that the string is assumed (although not garaunteed) to be valid under. The Interpreter will only be valid for context-free-grammars (CFG). A CFG consists of Terminals and Nonterminals. A Terminal is a character or sequence of characters that terminate a definition. A Nonterminal is a rule that defines the relationships between other Nonterminals and Terminals.

Design

interpreter_components

Interpreter designs differs from traditional compiler design, in that only the so-called front-end is implemented. There are three principal components involved in this process: Tokenizer, Parser and Semantic Analyzer.

A Tokenizer takes a string of characters and returns a collection of Tokens, where each Token is a tuple containing the accepted substring and meta data describing what the substring represents. Tokens are identified by looking for Terminals in the input string.

The Parser then builds (formally) a concrete syntax tree (informally a syntax tree) following the rules defined by the grammar. This is done by looking for Nonterminals in the collection of Tokens.

Finally, the Semantic Analyzer builds (formally) an abstract syntax tree (informally a semantic tree) by ignoring the syntax in the parse tree and constructing the appropriate tree for performing functions upon.

In the event of an exception, the Interpreter will write the details of the exception to a text writer (e.g., Console.Error) provided by the client of the Interpreter.

Exceptions are to be thrown under the following conditions:

The input string is null or empty, or none of the Terminals are able to accept a character in the input string.
None of the Nonterminals are able to accept a set of tokens, or the number of tokens is deficient for the Nonterminals to identify a match.

Implementation

abstract public class Tokenizer<T> {
    protected TerminalList<T> Terminals { get; private set; }

    public Tokenizer() {
        Terminals = new TerminalList<T>();
        loadTerminals();
    }

    public TokenList<T> Tokenize(string input) {
        if (string.IsNullOrEmpty(input))
            throw new UnexpectedEndOfInputException();

        TokenList<T> tokens = new TokenList<T>();
        Token<T> candiate;

        for (int n = 0; n < input.Length; n++) {
            if (char.IsWhiteSpace(input[n]))
                continue;

            if (!terminalExists(input, n, out candiate))
                throw new UnexpectedInputException(n);

            tokens.Add(candiate);
            n += candiate.Lexeme.Length - 1;
        }

        return tokens;
    }

    private bool terminalExists(string input, int n, out Token<T> candiate) {
        candiate = null;
        foreach (Terminal<T> terminal in Terminals)
            if (terminal.Exists(input, n, out candiate))
                return true;
        return false;
    }

    abstract protected void loadTerminals();
}

abstract public class Parser<T, P> {
    protected NonterminalList<T, P> Nonterminals { get; private set; }
    protected T[] FirstExpectedTokens { get; private set; }

    public Parser() {
        Nonterminals = new NonterminalList<T, P>();
        loadNonterminals();

        FirstExpectedTokens = new T[Nonterminals.Count];
        for (int n = 0; n < Nonterminals.Count; n++)
            FirstExpectedTokens[n] = Nonterminals[n].FirstExpectedToken;
    }

    public P Parse(TokenList<T> tokens) {
        return Parse(tokens, 0, tokens.Count - 1);
    }

    public P Parse(TokenList<T> tokens, int n, int m) {
        P parsedSentence;
        foreach (Nonterminal<T, P> nonterminal in Nonterminals)
            if (nonterminal.Exists(tokens, n, m, out parsedSentence))
                return parsedSentence;

        throw new UnexpectedTokenException<T>(tokens[n], FirstExpectedTokens);
    }

    abstract protected void loadNonterminals();
}

abstract public class SemanticAnalyzer<P, S> {
    abstract public S Analyze(P sentence);
}

abstract public class Interpreter<T, P, S> {
    private TextWriter externalLogger;

    abstract protected Tokenizer<T> Tokenizer { get; }
    abstract protected Parser<T, P> Parser { get; }
    abstract protected SemanticAnalyzer<P, S> SemanticAnalyzer { get; }

    public Interpreter(TextWriter externalLogger) {
        this.externalLogger = externalLogger;
    }

    public bool Interpret(string candidateSentence, out S sentence, out TokenList<T> tokens) {
        sentence = default(S);
        tokens = null;

        try {
            tokens = Tokenizer.Tokenize(candidateSentence);
            sentence = SemanticAnalyzer.Analyze(Parser.Parse(tokens));

            return true;
        } catch (UnexpectedInputException uie) {
            logger(uie.Message);
            logger("'{0}'", candidateSentence);
            logger(" {0}^", string.Empty.PadLeft(uie.AtIndex, '.'));

        } catch (UnexpectedEndOfInputException ueie) {
            logger(ueie.Message);
            logger("'{0}'", candidateSentence);
            logger(" {0}^", string.Empty.PadLeft(candidateSentence.Length, '.'));

        } catch (UnexpectedTokenException<T> pe) {
            logger(pe.Message);
            logger("'{0}'", candidateSentence);
            logger(" {0}^", string.Empty.PadLeft(pe.Unexpected.StartsAtIndex, '.'));
        }

        return false;
    }

    private void logger(string format, params object[] values) {
        externalLogger.WriteLine(format, values);
    }
}

Example

Specification

SL Backus-Naur Form (BNF)
sentence	::=	sentence_letter
	\|	~ sentence
	\|	( sentence connective sentence )

sentence_letter	::=	letter
sentence_letter	\|	letter number

connective	::=	v
	\|	^
	\|	->
	\|	<->

letter	::=	a
	\|	…
	\|	z
	\|	A
	\|	…
	\|	Z

number	::=	number digit
number	\|	digit

digit	::=	0
	\|	…
	\|	9

Sentential Logic (SL) is a simple formal system that consists of two semantic elements: primitive types called Sentence Letters that may either be true or false, and expressions built upon Sentence Letters called Sentences. A Sentence Letter is any character a-z (lower or upper) followed by an optional number. Sentences are defined in terms of other sentences, the simplest is a plain Sentence Letter. Any negated ¬ sentence is also a sentence. Any two sentences connected by a conjunction ∧, disjunction ∨, material conditional → or material biconditional ↔ and enclosed by parentheses () is also a sentence under SL. No other arrangement of symbols constitutes a sentence under SL.

Connectives and negation Truth Tables
	¬	∧	T	F	∨	T	F	→	T	F	↔	T	F
T	F		T	F		T	T		T	F		T	F
F	T		F	F		T	F		T	T		F	T

A truth table has a column for each of the sentence letter that show up in a sentence and a column for the sentence. Each row is one Cartesian coordinate in the space $\{T, F \}^{n}$ where n is the number of the number sentence letters. By convention, one lists the rows in order all sentence letters being true, to the last row where all sentence letters are false.

Implementation

TerminalSyntax and TerminalSentenceLetter are subclasses of Terminal and implement the logic for determining if a certain syntax or sentence letter exists in an input string. TerminalSyntax takes in the TokenType and string to find in the input as configuration options. If one felt inclined, one could have implemented a TerminalRegularExpression class instead of two separate classes but I felt it was unnecessary.

public class SLTokenizer : Tokenizer<TokenType> {
    protected override void loadTerminals() {
        Terminals.AddRange(new Terminal<TokenType>[] { 
            new TerminalSyntax(TokenType.Conjunction, "^"),
            new TerminalSyntax(TokenType.Disjunction, "v"), 
            new TerminalSyntax(TokenType.LeftParen, "("),
            new TerminalSyntax(TokenType.MaterialBiconditional, "<->"),
            new TerminalSyntax(TokenType.MaterialConditional, "->"),
            new TerminalSyntax(TokenType.Negation, "~"),
            new TerminalSyntax(TokenType.RightParen, ")"),
            new TerminalSentenceLetter()
        });
    }
}

The three Nonterminal types for Connectives, Negation and Letters are subclasses of Nonterminal and implement the logic for determining if a collection of tokens matches the rules associated with each Nonterminal type.

public class SLParser : Parser<TokenType, ParsedSentence> {
    protected override void loadNonterminals() {
        Nonterminals.AddRange(new Nonterminal<TokenType, ParsedSentence>[] { 
            new NonterminalSentenceConnective(this),
            new NonterminalSentenceNegation(this),
            new NonterminalSentenceLetter()
        });
    }
}

Nothing particularly exciting here, simply mapping the parsed sentences to actual sentences that will be used by the client code. Sentence implementations have the logic for evaluating themselves.

public class SLSemanticAnalyzer : SemanticAnalyzer<ParsedSentence, Sentence> {
    public override Sentence Analyze(ParsedSentence sentence) {
        ParsedSentenceLetter letter = sentence as ParsedSentenceLetter;
        if (letter != null)
            return new SentenceLetter(letter);

        ParsedSentenceUnary negation = sentence as ParsedSentenceUnary;
        if (negation != null)
            return new SentenceUnaryNegation(negation, Analyze(negation.Sentence));

        ParsedSentenceBinary connective = sentence as ParsedSentenceBinary;
        if (connective != null) {
            switch (connective.Operation.TokenType) {
                case TokenType.Conjunction:
                    return new SentenceBinaryConjunction(connective, Analyze(connective.Left), Analyze(connective.Right));
                case TokenType.Disjunction:
                    return new SentenceBinaryDisjunction(connective, Analyze(connective.Left), Analyze(connective.Right));
                case TokenType.MaterialBiconditional:
                    return new SentenceBinaryMaterialBiconditional(connective, Analyze(connective.Left), Analyze(connective.Right));
                case TokenType.MaterialConditional:
                    return new SentenceBinaryMaterialConditional(connective, Analyze(connective.Left), Analyze(connective.Right));
            }
        }

        throw new NotImplementedException();
    }
}

Nothing extra to implement in the SLInterpreter to add other than the constructor and properties.

public class SLInterpreter : Interpreter<TokenType, ParsedSentence, Sentence> {
    public SLInterpreter(TextWriter externalLogger) : base(externalLogger) { 
    
    }
    
    override protected Tokenizer<TokenType> Tokenizer {
        get { return new SLTokenizer(); }
    }

    override protected Parser<TokenType, ParsedSentence> Parser {
        get { return new SLParser(); }
    }

    override protected SemanticAnalyzer<ParsedSentence, Sentence> SemanticAnalyzer {
        get { return new SLSemanticAnalyzer(); }
    }
}

Once all of the Interpreter code has been completed, it becomes trivial to go and implement the truth table.

private void table(string input) {
    SLInterpreter interpreter = new SLInterpreter(Console.Out);
    TokenList<TokenType> tokens;
    Sentence sentence;

    if (!interpreter.Interpret(input, out sentence, out tokens))
        return;

    List<Variable> variables = tokens
       .Where((x) => x.TokenType == TokenType.Letter)
       .Select((x) => x.Lexeme)
       .Distinct()
       .Select((x) => new Variable(x))
       .ToList();

    foreach (Variable variable in variables)
        Console.Write("{0} ", variable.Identifier);
    Console.WriteLine("| {0}", input);

    for (int n = (1 << variables.Count) - 1; n >= 0; n--) {
        int N = n;
        for (int m = variables.Count - 1; m >= 0; m--) {
            variables[m].Value = (N & 1) == 1;
            N >>= 1;
        }

        IDictionary<Token<TokenType>, bool> tableRow = new Dictionary<Token<TokenType>, bool>();
        sentence.Evaluate(variables, tableRow);

        int offset = 2;
        foreach (Variable variable in variables) {
            Console.Write("{0} ", variable.Value ? "T" : "F");
            offset += 1 + variable.Identifier.Length;
        }

        Console.Write("| ");
        foreach (Token<TokenType> token in tableRow.Keys) {
            Console.CursorLeft = token.StartsAtIndex + offset;
            Console.Write(tableRow[token] ? "T" : "F");
        } Console.WriteLine();
    }
}

Written by lewellen

2009-03-01 at 12:00 am

Posted in Logic, Programming Languages

Tagged with C-Sharp, Interpreter Design, Propositional Logic

The ring of Gaussian Integers

leave a comment »

I recently started taking a foray into Complex Analysis as a means of filling in the gaps of my undergraduate mathematics knowledge. After reading about holomorphic functions, the Cauchy-Riemann equations and how to model ideal fluid flows I decided to take a break to digest it all. During this period I was reflecting on what I had done with Number Theory the previous summer and the question popped in my head: what is the complex equivalent of $w \equiv v \pmod{z}$ ? Or for that matter, are integer primes also primes in the complex domain (and vice versa)? And, are all complex $z = \displaystyle\prod_{k = 0}^{\infty} p_{k}^{e_{k}}$ factorizations unique?

After looking around the internet for a while I came across a rather well written paper [pdf] by Assistant Professor Keith Conrad of the University of Connecticut that answered all of the questions brewing in my head along with the ones that weren’t. I’m going to surmise the basic ideas behind the solutions Conrad wrote about, if you have the free time it’s worthwhile to read the original text in its entirety.

If we consider the set of complex numbers of the form $n + mi$ and restrict $n, m \in \mathbb{Z}$ we have the set of the Gaussian Integers $\mathbb{Z}[i]$ . If we take two values $z, w \in \mathbb{Z}[i]$ we can say that $z|w$ only if $w = zu$ where $u \in \mathbb{Z}[i]$ . Thus, $\displaystyle\bar{z}w = \bar{z}zu \Rightarrow u = \frac{\bar{z}w}{\bar{z}z}$ . With this last step we’ve reduced the problem to satisfying two conditions: $\displaystyle\bar{z}z|\Re{\bar{z}w} \wedge \bar{z}z|\Im{\bar{z}w}$ . For example let $z = 1 + 2i, w = -5 + 10i \Rightarrow \bar{z}z = 5, \bar{z}w = 15 + 20i$ which satisfies the conditions $5|15 \wedge 5|20 \therefore z|w$ . For the congruence $w \equiv v \pmod{z}$ with $w,v,z \in \mathbb{Z}[i]$ to hold true, it must be the case that $w = zu + v \Rightarrow z|(w - v)$ .

Under $\mathbb{Z}$ any integer $m$ that is divisible by any integer other than a unit $\pm 1$ or unit multiple $\pm m$ is said to be composite, otherwise $m$ is said to be prime. The same definition holds in $\mathbb{Z}[i]$ but our units are now $\pm 1, \pm i$ and our unit multiples are now $\pm z, \pm zi$ . As one might have guessed all primes in $\mathbb{Z}$ are primes in $\mathbb{Z}[i]$ however, the converse does not hold. For example $5$ is prime but can be written as $(1 + 2i)(1 - 2i)$ and is thus a composite number under $\mathbb{Z}[i]$ . It is also useful to mention that if $\bar{z}z \in \mathbb{P}$ then $z$ is prime. By way of the Fundamental Theorem of Arithmetic, every composite integer can be written as the product of prime integers, this is also the case under $\mathbb{Z}[i]$ and each factorization is unique.

The mechanics of finding the prime factorization of $m \in \mathbb{Z}^{+}$ can be naïvely done by trial division from $n = 2$ to $\sqrt{m}$ . Once $n|m$ reset $n = 2$ and $m_{k} = m_{k - 1} / n$ and continue this procedure until $m_{k} = 1$ . That is all well and fine for $\mathbb{Z}^{+}$ but it isn’t immediately useful for doing the same under $\mathbb{Z}[i]$ .

For example, let’s consider $z = -1395 - 12410i$ , since we are uncertain about how to factor $z$ , lets think about how to factor $\bar{z}z$ as we already know how to factor under $\mathbb{Z}$ and see where that takes us. $\bar{z}z = 155954125 = 5^{3} \cdot 61 \cdot 113 \cdot 181$ We went from $\mathbb{Z}[i] \to \mathbb{Z}$ so we’d hope that we could go the other way around by finding $p = (n + mi)(n - mi)$ for each of the prime factors of $\bar{z}z$ . In other words, we want to know $p = n^{2} + m^{2}$ which will produce four Gaussian primes $n \pm mi, m \pm ni$ . Thus, we find $5 = 1 + 2^{2}$ , $61 = 5^{2} + 6^{2}$ , $113 = 7^{2} + 8^{2}$ , $181 = 9^{2} + 10^{2}$ . Now, according to Conrad we should be able to pick any one of the four possible Guassian primes for each of the prime factors and multiply each Gaussian prime across, in doing so we get $z = (2 + i)^{2} (2 - i) (6-5i) (8 - 7i) (10 - 9i)$ . Now that’s pretty damn cool.

The general procedure for factoring a Gaussian Integer $z$ into it’s prime components is to find the integer factorization of $\bar{z}z$ and for each prime factor find $p = n^{2} + m^{2}$ . Next, explore the Cartesian space formed by each factor’s four Gaussian primes until you come across a product coordinate that equals $z$ and output said coordinate. As one can image this is a woefully poor algorithm for find the prime factorization of $z$ . If anyone has ideas or knows of more efficient methods post them in the comments.

There is clearly much, much more to be said about Gaussian Integers, but this feels like a good stopping point. If you want to find out more about how some of the traditional Number Theory constructs are defined under $\mathbb{Z}[i]$ you’ll want to read the entirety of Conrad’s paper or jump on Google and see what’s out there.

Written by lewellen

2008-08-10 at 8:05 pm

Posted in Number Theory

Tagged with Gaussian Integers

Archive for the ‘Mathematics’ Category

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

Motivation

Abstract Algebra

Group Theory

Ring Theory

Module Theory

Representation of Value Types

Enum Value Types and the Dihedral Group

Integral Value Types and the Commutative Ring with Unity over

Floating-point Value Types and the Real Vector Space

Representation of Reference Types

Strings, Computability and Monoids

Classes, Interfaces, Type Theory and Semi-rings

Delegates and Process Algebras

Examples

Conclusion

References

Introduction

Preliminaries

k-ary Cartesian Product

k-Permutations

Power Set

k-Combinations

Introduction

Framework

Specification

Design

Implementation

Example

Specification

Implementation

Pages

Categories

Tags

Archives

Meta

Follow “Antimatroid, The”

Enum Value Types and the Dihedral Group $D_8$

Integral Value Types and the Commutative Ring with Unity over $\mathbb{Z} / 2^n \mathbb{Z}$

Floating-point Value Types and the Real Vector Space $\mathbb{R}^n$