\documentclass[11pt]{article}
\usepackage{fullpage,amsmath,amsfonts,mathpazo,microtype,nicefrac,algorithm2e,graphicx}

% Set-up for hypertext references
\usepackage{hyperref,color,textcomp}
\definecolor{webgreen}{rgb}{0,.35,0}
\definecolor{webbrown}{rgb}{.6,0,0}
\definecolor{RoyalBlue}{rgb}{0,0,0.9}
\hypersetup{
   colorlinks=true, linktocpage=true, pdfstartpage=3, pdfstartview=FitV,
   breaklinks=true, pdfpagemode=UseNone, pageanchor=true, pdfpagemode=UseOutlines,
   plainpages=false, bookmarksnumbered, bookmarksopen=true, bookmarksopenlevel=1,
   hypertexnames=true, pdfhighlight=/O,
   urlcolor=webbrown, linkcolor=RoyalBlue, citecolor=webgreen,
   pdfauthor={Chris H. Rycroft},
   pdfsubject={Harvard AM225 (Spring 2021)},
   pdfkeywords={},
   pdfcreator={pdfLaTeX},
   pdfproducer={LaTeX with hyperref}
}
\hypersetup{pdftitle={AM225: Assignment 1}}

\SetKwProg{Fn}{}{}{}\SetKwFunction{FRecurs}{void SudokuSolve}
\SetKwProg{Fm}{}{}{}\SetKwFunction{FRecursb}{void PuzzleSolve}

\graphicspath{{am225_hw1_files/}}

% Macro definitions
\newcommand{\N}{\mathbb{N}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\Q}{\mathbb{Q}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\p}{\partial}
\renewcommand{\vec}[1]{\mathbf{#1}}
\newcommand{\vx}{\vec{x}}
\newcommand{\vp}{\vec{p}}
\newcommand{\Atol}{\textit{Atol}}
\newcommand{\Rtol}{\textit{Rtol}}
\newcommand{\Trans}{\mathsf{T}}

\begin{document}
\section*{AM225: Assignment 1 (due 5~PM, February 11)}
\noindent \textit{This homework is designed to get you up to speed with
multithreaded programming in OpenMP, and highlight interesting challenges that
can occur. Complete at least \textbf{four out of six} problems. If you submit
answers for more, your grade will be calculated using the four best scores.
The questions are roughly ordered from easiest to hardest.}

\begin{enumerate}
  \item \textbf{A randomized Casino game.} Consider a casino game that costs
    \$250 to play. The casino draws successive numbers from a uniform
    distribution $U(0,1)$ until the sum exceeds 1. The player then receives
    \$100 for each number that was drawn.
    \begin{enumerate}
      \item Using OpenMP, write a multithreaded code that simulate $10^9$
	instances of this game, and calculate the expected amount of winnings,
	reporting your answer to at least five significant figures. Measure the
	wall clock time of your code using 1, 2, and 4 threads. Based on the
	expected winnings is this a game worth playing?

	[Note: A difficulty with this problem is that the
	\href{http://www.cplusplus.com/reference/cstdlib/rand/}{C++ \texttt{rand}
	function} for generating random numbers cannot be accessed by multiple
	threads simultaneously, since it contains an internal state that becomes
	corrupted due to race conditions. In the program files, you will find
	a program \texttt{custom\_rng.cc} that provides a custom self-contained
	random number generator. You will need to make multiple instances of
	this random number generator, each of which is used by a separate
	thread.]
      \item \textbf{Optional.} Analytically derive the expected winnings.
    \end{enumerate}
  \item \textbf{Cellular automaton mazes.} Consider a cellular automaton on a
    regular periodic $m \times n$ grid, with cells indexed as $(i,j)$ where
    $i\in \{0,1,\ldots,m-1\}$ and $j\in \{0,1,\ldots n-1\}$. Each cell can
    either be alive or empty. Define $N_{i,j}$ to be the number alive neighbors
    of $(i,j)$ in the adjacent orthogonal and diagonal cells. Hence $0\le
    N_{i,j} \le 8$.

    In one generation, all the cells are simultaneously updated according to the
    following rules:
    \begin{itemize}
      \item an alive cell at $(i,j)$ survives if $N_{i,j} \in S_\textit{survive}$,
      \item an empty cell at $(i,j)$ becomes alive if $N_{i,j} \in S_\textit{birth}$.
    \end{itemize}

    \begin{enumerate}
      \item Implement this cellular automaton, using
	$S_\textit{survive}=\{1,2,3,4,5\}$ and $S_\textit{birth}=\{3\}$. Use a
	grid with $(m,n)=(80,40)$ that is initially empty except for the square
	\smash{$\frac{m}2-6\le i < \frac{m}2+6$} and \smash{$\frac{n}2-6 \le j
	< \frac{n}2+6$} where each cell is set alive with
	probability \smash{$\frac{3}{4}$}. Print a snapshot of the initial
	condition.\footnote{A simple way to do this is to use
	\href{https://en.wikipedia.org/wiki/ASCII_art}{ASCII art} and the
	\texttt{printf}, \texttt{puts}, \textit{etc.}\ functions
	presented in the lectures.} Perform 150 generations, and print a
	snapshot of the system after every 25 generations.
      \item Extend your program to use OpenMP to process the grid in parallel.
	Consider simulating on a square $n\times n$ grid, using the same
	initial condition as in part (a). Consider grid sizes of $n=16, 32, 64,
	128, 256, 512, 1024$ using $T=1, 2, 3, 4$ threads. For each combination
	of $n$ and $T$ measure the wall clock time $w(n,T)$ to compute a
	generation.\footnote{To get an accurate measure, you may need to average
	over the wall clock time to compute a number of generations. For smaller
	grids, more generations may be required.}
      \item For $T=2,3,4$, make a plot of the parallel efficiency
	$p(n,T)=w(n,1)/(T w(n,T))$ as function of $\log_2 n$.
      \item \textbf{Optional.} On a very large grid, what proportion of cells
	are alive once the system reaches a steady state?
      \item \textbf{Optional.} Investigate the alternative growth model
	$S_\textit{survive}=\{5,6,7,8\}$ and $S_\textit{birth}=\{3,5,6,7,8\}$
	and make a movie showing the typical patterns that emerge.
      \item \textbf{Optional.} Is it possible to write the code so that the
	memory $M$ required scales like $M \sim mn\text{~byte}$? Is possible to
	write the code so that $M \sim mn\text{~bit}$?
    \end{enumerate}
  \item \textbf{A Mersenne unprime.} The
    \href{https://en.wikipedia.org/wiki/Mersenne_prime}{Mersenne primes} are a
    special sequence of prime numbers of the form $2^n-1$ for an integer $n$.
    In December 2018, a new Mersenne prime of $M=2^{82589933}-1$ was
    discovered, which is currently the largest known prime number. In decimal
    notation it has 24,862,048 digits.

    The number $M$ is far too large to be represented using a standard computer
    integer. However, it can be represented in a long expansion as
    \begin{equation}
      \label{eq:long_expansion}
      M = \sum_{k=0}^K d_k B^k
    \end{equation}
    where $B$ is the base and $d_k\in \{0,1,\ldots,B-1\}$ are the
    digits.\footnote{Hence, if $B=10$ this is equivalent to the standard
    decimal notation for a number.} For this question, it is convenient to set
    $B$ as a large power of two that fits within a standard computer integer.
    \begin{enumerate}
      \item Write a program to find\footnote{This could be done using the
	\href{https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes}{Sieve of
	Eratosthenes}, for example. The Sieve is well-suited to being
	multithreaded.} all prime numbers less than $2 \times 10^5$.
      \item Write a program that performs division with remainder on a
	number represented in the format of Eq.~\ref{eq:long_expansion}.
	You can assume that the divisor is smaller than $B$.
      \item Use your program from part (b) to test how many primes less
	than $2\times 10^5$ are factors of $M$. Extend your program to use OpenMP
	so that several primes can be simultaneously tested using multiple
	threads.
      \item The number $N=2^{82589932}-1$ is presumably not prime. Using your
	program from (c), determine the number of primes less than $10^6$ that
	are factors of $N$.
      \item \textbf{Optional.} Compute the decimal expansion of $M$. Count the
	frequencies of the numbers $0,1,\ldots,9$ and show that they pass a
	$\chi^2$ test for randomness.
    \end{enumerate}
  \item \textbf{A multithreaded Sudoku solver.}
    \href{https://en.wikipedia.org/wiki/Sudoku}{Sudoku} is a popular puzzle
    that involves filling numbers from 1 to 9 in a $9\times 9$ grid according
    to simple rules. In magazines and newspapers, you will often find partially
    completed grids of numbers from which you can fill in the remaining numbers
    to uniquely complete the puzzle.

    Suppose that the grid squares are numbered from $0$ to $80$ in some order.
    Then a simple recursive algorithm to solve a Sudoku is shown in
    Algorithm~\ref{alg:sudoku}.
    \begin{enumerate}
      \item Implement Algorithm~\ref{alg:sudoku} and use it to find the unique
	solution to the partial grid shown in Fig.~\ref{fig:sudoku}(a). Find
	the wall clock time required to solve the puzzle. To gain accurate
	statistics, you may need to time $N$ repeats, and then take the average.
      \item Use your program to count the total number of solutions to the
	partial grid shown in Fig.~\ref{fig:sudoku}(b). Show that there are
	exactly 283576 solutions.
      \item Extend your program from the previous section to use OpenMP.
	Calculate wall clock times using 1, 2, and 4 threads.
      \item Find the total number of solutions to the partial grid shown in
	Fig.~\ref{fig:sudoku}(c).
      \item \textbf{Optional.} It turns out that many Sudoku puzzles found
	online and in magazines have redundancies, meaning that some numbers
	can be removed while still achieving a unique solution. Find an example
	of real puzzle with a redundancy and use your code to demonstrate this.
      \item \textbf{Optional.} Use your program to collect some interesting
	statistics, such as how the number of solutions scales with the number
	of initially filled-in squares.
      \item \textbf{Optional.} Implement some algorithmic improvements to the
	solver.
    \end{enumerate}
  \item \textbf{Testing thread performance on the diffusion equation.}
    Consider the diffusion equation
    \begin{equation}
      \label{eq:diff_untrans}
      \frac{\p u}{\p t} = \nabla^2 u
    \end{equation}
    in the two-dimensional plane with coordinate system $\vec{x}=(x,y)$.
    Suppose that instead of solving Eq.~\ref{eq:diff_untrans} directly, the
    equation is solved in a transformed coordinate system $\vec{X}=(X,Y)$
    where $\vec{x}=T\vec{X}$ for some constant $2\times 2$ matrix $T$.
    Show that Eq.~\ref{eq:diff_untrans} can be rewritten as
    \begin{equation}
      \label{eq:diff_trans}
      \frac{\p u}{\p t} = \alpha \frac{\p^2 u}{\p X^2} + \beta \frac{\p^2 u}{\p
      X \p Y} + \gamma \frac{\p^2 u}{\p Y^2}
    \end{equation}
    and determine $\alpha$, $\beta$, and $\gamma$ in terms of the components of
    $T$. Consider discretizing Eq.~\ref{eq:diff_trans} using an explicit finite
    difference scheme with equal spacing $h$ in the $X$ and $Y$ directions, and
    timestep $\Delta t$ in time, so \smash{$u^{n}_{j,k} = u(jh,kh,n\Delta t)$}.
    Use an explicit Euler step for $\p u/\p t$, the standard centered finite
    difference formula for $\p^2 u/\p X^2$ and $\p^2 u/\p Y^2$, and the stencil
    \begin{equation}
      \left[\frac{\p^2 u}{\p X \p Y} \right]^n_{j,k} = \frac{u^n_{j+1,k+1} - u^n_{j-1,k+1}-u^n_{j+1,k-1}+u^n_{j-1,k-1}}{4h^2}
    \end{equation}
    for the cross-derivative term.
    \begin{enumerate}
      \item Suppose that $\alpha=\gamma$. Perform a Fourier stability analysis\footnote{Notes from \href{https://courses.seas.harvard.edu/courses/am205/material.html}{AM205 Unit 3} may be useful.} by
	substituting in the ansatz \smash{$u^n_{j,k}=[\lambda(l,m)]^n e^{i(jl+km)h}$}
	into the discretized formula, where $\lambda(l,m)$ is the amplification
	factor. You may find it useful to express $\lambda$ in terms of
	\begin{equation}
	  q= \cos \frac{(l+m)h}{2}, \qquad r = \cos \frac{(l-m)h}{2}.
	\end{equation}
	Show that there exists a $\beta_*$ such that if $|\beta|>\beta_*$ the
	method is unconditionally unstable, but that if $|\beta| \le \beta_*$
	then the method is conditionally stable. In the latter case, determine
	the maximum allowable timestep.
      \item Write a multithreaded code using OpenMP that integrates this
	scheme on a $N \times N$ grid by dividing the gridpoint updates among
	$P$ threads. Use the periodic domain $[0,1]^2$, the matrix
	\begin{equation}
	  T = \left(
	  \begin{array}{cc}
	    1 & \nicefrac{1}{2} \\
	    0 & 1
	  \end{array}
	  \right)
	  \label{eq:transform}
	\end{equation}
	and the initial condition
	\begin{equation}
	  u(X,Y,0)=e^{-\cos 2\pi X -\cos 2\pi Y}.
	\end{equation}
	Choose a timestep that satisfies the restriction in part (a). Find the
	wall clock time to run 1000 timesteps for $N=64,128,256,512,1024,2048$
	and $P=1,2,3,4$.
      \item Define
	\begin{equation}
	  \langle u^k \rangle(t) = \int_0^1 dX \int_0^1 dY \,[u(X,Y,t)]^k
	\end{equation}
	and
	\begin{equation}
	  S(t) = \sqrt{\langle u^2\rangle - [\langle u\rangle]^2}.
	\end{equation}
	Using $N=400$ and your choice of threads, plot $S(t)$ over the range
	$0\le t \le \nicefrac1{10}$ for the transform $T$ and initial
	conditions used in part (b). On the same graph, plot $S(t)$ when the
	transform $T$ is just the identity. Which of the two $S(t)$ curves
	decreases more rapidly, and why?
      \item \textbf{Optional.} Extend your analysis from part (a) to the case
	when $\alpha \ne \gamma$.
    \end{enumerate}
  \item \textbf{A wooden puzzle.} Chris has a large collection of wooden
    puzzles, which he mainly received as gifts from his father. Recently, Chris
    received the
    \href{https://www.professorpuzzle.com/products/einsteins-letter-puzzle/}{\textit{Letter
    Block Puzzle}} in a series from Professor
    Puzzle, a UK-based company. It consists of fourteen wooden pieces that
    spell out the letters of ``ALBERT EINSTEIN'' (Fig.~\ref{fig:puz_pieces}).
    The puzzle contains several different challenges, and the hardest one is to
    arrange the pieces into $4\times 4 \times 5$ cuboid. This can be done with
    a recursive algorithm.
    \begin{enumerate}
      \item Consider the nine distinct pieces in Fig.~\ref{fig:puz_pieces}. For
	each piece, write down the number of different orientations in which it could
	be placed into the cuboid. You should find that the pieces can be
	divided into four sets with different orientation groups.
      \item Let the individual blocks in the $4 \times 4\times 5$ cuboid be
	labeled from $i=0, 1, \ldots, 79$ in some order.\footnote{For example,
	if the blocks are indexed as $(r,s,t)$ with $r\in \{0,3\}$, $s \in
	\{0,3\}$, and $t\in \{0,4\}$, then the blocks could be labeled using
	$i=r+4s+16t$.} Let the nine distinct pieces be labeled $j=0,1,\ldots,8$
	and let $N_j$ be the number of available pieces of type $j$. Using Algorithm
	\ref{alg:woodpuz} write a
	program to find a valid arrangement of the pieces in the $4 \times
	4\times 5$ cuboid.
      \item Extend your program to use OpenMP, and find the total number of
	ways that the pieces can be assembled to make the cuboid.
      \item \textbf{Optional.} Extend your program to find the following solutions:
	\begin{enumerate}
	  \item Using two complete sets of pieces, assemble them into a
	    $2\times 8 \times 10$ cuboid.
	  \item Using three complete sets of pieces, assemble them into a
	    $5 \times 6 \times 8$ cuboid.
	  \item Using four complete sets of pieces, assemble them into the
	    shape in Fig.~\ref{fig:xtreme_challenge}.
	\end{enumerate}
	For these cases, the search tree of possible solutions becomes
	considerably larger, and you may need to implement code optimizations
	in order to find.
    \end{enumerate}
\end{enumerate}

\RestyleAlgo{boxed}
\begin{algorithm}
  \SetAlgoLined
  \Fn{\FRecurs($i$)}{
    \eIf{$i=81$}{
      The grid is complete\;
    }{
      \eIf{square $i$ is empty}{
	\For{$j\in\{1,2,\ldots,9\}$}{
	  \If{$j$ allowed at square $i$}{
	    Add $j$ to square $i$\;
	    Call \texttt{SudokuSolve}($i+1$)\;
	    Remove $j$ from square $i$\;
	  }
	}
      }{
	Call \texttt{SudokuSolve}($i+1$)
      }
    }
  }
  \caption{A recursive Sudoku solver. Here the squares in the $9\times 9$
  Sudoku grid are labeled from $i=0$ to $i=80$.\label{alg:sudoku}}
\end{algorithm}

\begin{figure}
  \begin{center}
    \includegraphics[width=\textwidth]{am225_hw1_files/sudoku}
  \end{center}
  \vspace{-1em}
  \caption{Three partial Sudoku grids.\label{fig:sudoku}}
\end{figure}

\newlength{\pheight}
\setlength{\pheight}{3.6cm}
\begin{figure}
  \begin{center}
    \begin{tabular}{cccc}
      \includegraphics[height=\pheight]{puz_pieces/charA} &
      \includegraphics[height=\pheight]{puz_pieces/charT} &
      \includegraphics[height=\pheight]{puz_pieces/charR} &
      \includegraphics[height=\pheight]{puz_pieces/charN} \\
      A & T ($\times 2$) & R & N ($\times 2$)
    \end{tabular}
    \begin{tabular}{ccccc}
      \includegraphics[height=\pheight]{puz_pieces/charL} &
      \includegraphics[height=\pheight]{puz_pieces/charB} &
      \includegraphics[height=\pheight]{puz_pieces/charE} &
      \includegraphics[height=\pheight]{puz_pieces/charI} &
      \includegraphics[height=\pheight]{puz_pieces/charS} \\
      L & B & E ($\times 3$) & I ($\times 2$) & S
    \end{tabular}
  \end{center}
  \vspace{-1em}
  \caption{The nine pieces that make up the wooden puzzle in question 6. The brackets show how many times that piece is repeated.\label{fig:puz_pieces}}
\end{figure}

\RestyleAlgo{boxed}
\begin{algorithm}
  \SetAlgoLined
  \Fm{\FRecursb($i$)}{
    \eIf{$i=80$}{
      The grid is complete\;
    }{
      \eIf{block $i$ is empty}{
	\For{$j\in\{0,1,\ldots,8\}$}{
	  \If{$N_j>0$}{
	    \For{all valid positions $P$ of piece $j$ at $i$} {
	      Insert piece $j$ into grid with position $P$\;
	      Decrement $N_j$\;
	      Call \texttt{PuzzleSolve}($i+1$)\;
	      Increment $N_j$\;
	      Remove piece $j$ from grid with position $P$\;
	    }
	  }
	}
      }{
	Call \texttt{PuzzleSolve}($i+1$)
      }
    }
  }
  \caption{A recursive algorithm to solve the wooden puzzle. Here, the blocks in
  the $4\times 4\times 5$ cuboid are labeled from $i=0$ to $i=79$. $N_j$ is the
  number of available pieces of type $j$.\label{alg:woodpuz}}
\end{algorithm}

\begin{figure}
  \begin{center}
    \includegraphics[width=0.6\textwidth]{puz_pieces/triangle}
  \end{center}
  \vspace{-2em}
  \caption{A target shape to make using four complete sets of the
  \textit{Letter Block Puzzle} with pieces shown in
  Fig.~\ref{fig:puz_pieces}.\label{fig:xtreme_challenge}}
\end{figure}

\end{document}
