MTHS 3600 Project 2

Three Puzzles

More Computational and Experimental Mathematics

February 28, 2014

Team Members: _____________________________

_________________________________

_____________________________

Triangle Puzzle

Consider an isoceles triangle with 4 slots on each side and 5 slots on the base. We can place the integers 0 through 9 on the slots of the triangle as follows.

1 9

2 8

3 4 5 6 7

In this case the left side adds up to 6, the right side adds up to 24, and the base adds up to 25. The goal of the triangle puzzle is to find arrangements of the integers 0 through 9 so that all three sides add up to the same number.

The goal of the programming exercise is to write a program which will generate ALL the distinct solutions to the triangle puzzle in the most efficient way possible --- faster programs get higher scores. My personal program is the bar to shoot for and it is almost three orders of magnitude faster than Triangle_Puzzle01.

The following block of code builds on the triple for loop solution to the sampling without replacement exercise in HW09. In this case it enumerates all $10!$ permutions of the integers 0 through 9 and then checks the sums for the three sides. If the sums are equal, then the count of the number of solutions is incremented. The first 5 solutions are printed out and the total number of solutions is returned.

def Triangle_Puzzle01(): available = [True for k in range(10)] count = 0 for t0 in range(10): available[t0] = False for t1 in range(10): if available[t1]: available[t1] = False for t2 in range(10): if available[t2]: available[t2] = False for t3 in range(10): if available[t3]: available[t3] = False for t4 in range(10): if available[t4]: available[t4] = False for t5 in range(10): if available[t5]: available[t5] = False for t6 in range(10): if available[t6]: available[t6] = False for t7 in range(10): if available[t7]: available[t7] = False for t8 in range(10): if available[t8]: available[t8] = False for t9 in range(10): if available[t9]: left_sum = t0+t1+t2+t3 bottom_sum = t3+t4+t5+t6+t7 right_sum = t7+t8+t9+t0 if left_sum==right_sum and right_sum==bottom_sum: count += 1 if count < 6: print ' ',t0 print ' ',t1,t9 print ' ',t2,' ',t8 print t3,t4,t5,t6,t7 print available[t8] = True available[t7] = True available[t6] = True available[t5] = True available[t4] = True available[t3] = True available[t2] = True available[t1] = True available[t0] = True return count

time Triangle_Puzzle01()

There are two major strategies for improving the speed of a program like this.

The first strategy is called Pruning. Think of a tree with a single root node that represents the empty solution. From this node we have 10 possible children. For each of these children, there are 9 possible grandchildren. And so forth. At the tenth level, the bottom of this tree, there are $10!$ leaves. Pruning is the strategy of cutting off branches of the tree as soon as we know that that branch can not lead to a solution. This strategy includes several implementation tactics:

Do the equality tests earlier
Rearrange the order in which you choose entries for each slot so that tests can be done earlier
Exploit any mathematical properties of the puzzle that might constrain the values that can be placed in certain slots.

The second strategy exploits Symmetries. This puzzle is symmetric about the vertical axis. Consequently we could insist that the lower left corner always be smaller (or larger) than the lower right corner. Once we have enumerated all the solutions that satisfy this constraint, then we only have to multiply by two to find the total number of solutions. In this puzzle there are other symmetries that can also be exploited.

In the following blocks enter a copy of a working code that implements one of these ideas. Run it and verify that the timing has improved. Use a separate block for each improved implementation, and be sure to include comments or preliminary text descriptions of your rationale for each improvement.

def Triangle_Puzzle02()

time Triangle_Puzzle02()

def Triangle_Puzzle03()

time Triangle_Puzzle03()

def Triangle_Puzzle04()

time Triangle_Puzzle04()

def Triangle_Puzzle05()

time Triangle_Puzzle05()

Exploring Mersenne Primes

Many early writers felt that the numbers of the form $2^p-1$ were prime for all primes $p$, but in 1536 Hudalricus Regius showed that 2¹¹-1 = 2047 was not prime (it is 23^.89). By 1603 Pietro Cataldi had correctly verified that 2¹⁷-1 and 2¹⁹-1 were both prime, but then incorrectly stated $2^p-1$ was also prime for 23, 29, 31 and 37. In 1640 Fermat showed Cataldi was wrong about 23 and 37; then Euler in 1738 showed Cataldi was also wrong about 29. Sometime later Euler showed Cataldi's assertion about 31 was correct. The French monk Marin Mersenne (1588-1648) stated in the preface to his Cogitata Physica-Mathematica (1644) that the numbers $2^p-1$ were prime for

p = 2, 3, 5, 7, 13, 17, 19, 31, 67, 127 and 257

and were composite for all other positive integers $p<257$. Mersenne's (incorrect) conjecture fared only slightly better than earlier conjectures, but it still got his name attached to these numbers.

Definition: When $2^p-1$ is prime it is said to be a Mersenne prime.

To underscore the enormous advantage that we have over these earlier mathematicians because of computers, the next two blocks use the Sage built-in function, is_prime(), to determine all the Mersenne Primes less than 1,000. Specifically, for all primes $p$<N, Mersenne1(N) determines whether $2^p-1$ is a prime and if it is, then it prints $p$, $2^p-1$, and the words Mersenne Prime.

def Mersenne1(N): small_primes = [p for p in range(N) if is_prime(p)] for p in small_primes: Mp = 2^p - 1 if is_prime(Mp): print p,'\t',Mp,'\tMersenne Prime'

time Mersenne1(1000)

Complete the definition of the function Mersenne2, which should use the Sage built-in functions, is_prime() and factor(), to determine for all primes $p<260$ whether $2^p-1$ is a Mersenne prime and if not what is its factorization. This function should print the results in a table that displays $p$, $2^p-1$ and either the words Mersenne Prime or the factorization.

def Mersenne2(N): # # # # # # #

time Mersenne2(260)

List the primes $p<= 257$ which contradict Mersenne's conjecture.

You may also wish to comment on the apparent difference in the computational complexity of the functions, is_prime() and factor().

The function, Lehmer, is defined recursively as

Lehmer(1) = 4
Lehmer(n) = (Lehmer(n-1))^2 - 2

If we start the calculation with $n = p-1$, where $p$ is an odd prime, then $2^p-1$ is a Mersenne prime if and only if Lehmer(n) mod $(2^p-1)$ is zero, that is, if and only if Lehmer($p-1)$ is divisible by $2^p-1$.

Start by completing the function, Lehmer1, using the recursive definition. Then define and run the function Mersenne3 with an argument of $N = 18$.

def Lehmer1(n): # # # # return s

def Mersenne3(N): odd_primes = [p for p in range(3,N) if is_prime(p)] for p in odd_primes: Mp = 2^p-1 s = Lehmer1(p-1) if s%Mp == 0: print p,'\t',Mp,'\tMersenne Prime','\t',s else: print p,'\t',Mp,'\t',s

Mersenne3(18)

The results from Mersenne3(18) clearly indicate that the value of $s$ is growing very rapidly and increasing $p$ will cause some problems. (Feel free to verify this for yourself.) Since we only want to determine that $s$ mod $(2^p-1)$ is 0 we can move the mod function inside the definition of Lehmer. To do this, we will need to have two arguments n and Mp. On the initial call n will have the value $p-1$ and Mp will have the value $2^p-1$. The argument n will be decremented on each recursive call but the argument Mp will stay the same. The values of s will never grow larger than the value of Mp and we can easily calculate Mersenne4(150).

def Lehmer2(n,Mp): # # # # return s

def Mersenne4(N): odd_primes = [p for p in range(3,N,2) if is_prime(p)] for p in odd_primes: Mp = 2^p-1 s = Lehmer2(p-1,Mp) if s%Mp == 0: print p,'\t',Mp,'\tMersenne Prime','\t',s else: print p,'\t',Mp,'\t',s

Mersenne4(150)

This seems like an adequate solution, but try Mersenne4(1000). You will encounter a different memory problem, one which is actually due to the memory limitation on recursive functions. A recursive function has to have enough memory to maintain copies of the state for every instantiation of the function. In this case 1,000 copies is too much. However, each instance of the Lehmer2 function just invokes itself and does a simple calculation with the result. This is known as tail recursion. Since we know that the base step returns 4, we can easily convert the recursion into a simple loop. Complete the definition of Lehmer3 using this simple loop structure, and see how well it runs on all the primes less than 1,000.

def Lehmer3(p): # # # # return s def Mersenne5(N): odd_primes = [p for p in range(3,N,2) if is_prime(p)] for p in odd_primes: Mp = 2^p-1 if Lehmer3(p-1)==0: print p,'\t',Mp,'\tMersenne Prime'

Mersenne5(1000)

What is the largest known Mersenne Prime?

When was it discovered and by whom?

The Graph Class

class Graph(object): '''Represents a graph''' def __init__(self, vertices, edges): '''A Graph is defined by its set of vertices and its set of edges.''' self.V = set(vertices) self.E = set([]) self.Adj = {} self.add_edges(edges) print '(Initializing a graph with %d vertices and %d edges)' % (len(self.V),len(self.E)) def add_vertices(self,vertex_list): ''' This method will add the vertices in the vertex_list to the set of vertices for this graph. Since V is a set, duplicate vertices will not be added to V. ''' for v in vertex_list: self.V.add(v) self.build_Adj() def add_edges(self,edge_list): ''' This method will add a list of edges to the graph It will insure that the vertices of each edge are included in the set of vertices (and not duplicated). It will also insure that edges are added to the list of edges and not duplicated. ''' for s,t in edge_list: if (s,t) not in self.E and (t,s) not in self.E: self.V.add(s) self.V.add(t) self.E.add((s,t)) self.build_Adj() def build_Adj(self): self.Adj = {} for v in self.V: self.Adj[v] = [] for e in self.E: s,t = e self.Adj[s].append(t) self.Adj[t].append(s) def Degree(self,v): return len(self.Adj[v]) def plot(self,Vc): # Vc is a dictionary with the coordinates for each vertex g = Graphics() # Create the list of coordinates Vcoords = [Vc[v] for v in self.V] # Include the vertices as red dots g += point2d(Vcoords,rgbcolor=(1,0,0),size=25) # Include the names of the vertices for v in self.V: g += text(v,(1.1*Vc[v][0], 1.1*Vc[v][1])) # Include the edges for edge in self.E: u,v = edge g += line([Vc[u],Vc[v]]) return g def build_Components(self): Visited = {} for v in self.V: Visited[v] = False components = [] Q = Queue() for v0 in self.V: if not Visited[v0]: comp = [] Visited[v0] = True Q.enqueue(v0) while 0 < len(Q): u = Q.dequeue() comp.append(v) for v in self.Adj[u]: if not Visited[v]: Visited[v] = True Q.enqueue(v) components.append(comp) return components def is_connected(self): Visited = {} for v in self.V: Visited[v] = False Q = Queue() for v0 in self.V: if not Visited[v0]: Visited[v0] = True Q.enqueue(v0) while 0 < len(Q): u = Q.dequeue() for v in self.Adj[u]: if not Visited[v]: Visited[v] = True Q.enqueue(v) break for v in self.V: if not Visited[v]: # The Graph is not connected return False # The Graph is connected return True

The Prime Maze

The following blocks define and display the Prime Maze.

# The dictionary of edges and weights for the Prime Maze WPM = {('A', 'B'): 2, ('A', 'C'): 8, ('A', 'D'): 10, ('A', 'E'): 5,\ ('B', 'D'): 14, ('B', 'J'): 17, ('C', 'D'): 12, ('C', 'G'): 3,\ ('C', 'H'): 10, ('D', 'H'): 4, ('E', 'F'): 2, ('E', 'N'): 6,\ ('F', 'G'): 6, ('F', 'O'): 6, ('G', 'O'): 4, ('H', 'T'): 8,\ ('I', 'J'): 8, ('I', 'K'): 4, ('I', 'L'): 6, ('J', 'M'): 10,\ ('J', 'S'): 6, ('K', 'P'): 12, ('K', 'Q'): 12, ('L', 'M'): 12,\ ('L', 'Q'): 8, ('M', 'R'): 8, ('N', 'O'): 2, ('N', 'V'): 4,\ ('O', 'V'): 12, ('P', 'T'): 6, ('P', 'U'): 4, ('Q', 'R'): 10,\ ('R', 'S'): 12, ('S', 'W'): 13, ('T', 'U'): 14, ('U', 'W'): 10,\ ('V', 'W'): 2 } # print 'WPM',WPM WW = copy(WPM) for edge in WPM.keys(): u,v = edge WW[(v,u)] = WPM[edge] # print 'WW',WW

# The coordinates for displaying the graph PMV = ['A','B','C','D','E','F','G','H','I','J','K','L','M','N','P','Q','R','S','T','U','V','W'] PMcoord = {'B':(-3,3),'A':(0,3),'E':(3,3)} PMcoord['J']=(-3,2); PMcoord['D']=(-1,2); PMcoord['C']=(1,2); PMcoord['F']=(2,2) PMcoord['I']=(-1,1); PMcoord['H']=(0,1); PMcoord['G']=(1,1) PMcoord['M']=(-2,0); PMcoord['L']=(-1,0); PMcoord['K']=(0,0); PMcoord['T']=(1,0); PMcoord['O']=(2,0); PMcoord['N']=(3,0) PMcoord['R']=(-2,-1); PMcoord['Q']=(-1,-1); PMcoord['P']=(0,-1); PMcoord['U']=(2,-1) PMcoord['S']=(-3,-2); PMcoord['W']=(0,-2); PMcoord['V']=(3,-2) # print 'PMcoord', PMcoord

# Construct the Graph of the Prime Maze and display it PMG = Graph([],WPM.keys()) print 'Graph of the Prime Maze' Gpm = PMG.plot(PMcoord) for edge in PMG.E: s,t = edge xs,ys = PMcoord[s] xt,yt = PMcoord[t] w = WPM[edge] Gpm += text(w,((xs+xt)/2, 0.2+(ys+yt)/2)) show(Gpm,axes=False)

Write a program that finds all the cycles through the Prime Maze that satisfy the following property.

The distance from the root node to every intermediate node is a prime number.

#def prime_maze(new_vertex, distance, count): # Is the distance a prime number? # If not return count # If it is, append this new_vertex to the path # In the early testing you probably should print the path # If new_vertex is the same as the first vertex in the path # then we have found a complete cycle # print the result, add 1 to the count and return count # Otherwise press on, check out - recursively - every edge connected to new_vertex # This will be very similar to the loop inside the Prime Maze Search Code. return count

# Prime Maze Search Code for start_node in PMG.V: path = [start_node] count = 0 for v in PMG.Adj[start_node]: edge = (start_node,v) count = prime_maze(v, WW[edge], count) print start_node, ': ', count

MATH 3600 Project 2

3382 days ago by warner