The Inner Machinations of Walker

Project Euler Problem #67

Problem #67 is the first instance of a problem that is just a more computationally intensive version of an earlier problem. The question reads:

Project Euler Problem 67: Maximum path sum II
By starting at the top of the triangle below and moving to adjacent numbers on the row below, the maximum total from top to bottom is 23.
3
7 4
2 4 6
8 5 9 3
That is, 3 + 7 + 4 + 9 = 23.
Find the maximum total from top to bottom in triangle.txt (right click and 'Save Link/Target As...'), a 15K text file containing a triangle with one-hundred rows.
NOTE: This is a much more difficult version of Problem 18. It is not possible to try every route to solve this problem, as there are 2⁹⁹ altogether! If you could check one trillion (10¹²) routes every second it would take over twenty billion years to check them all. There is an efficient algorithm to solve it. ;o)

This problem is exactly the same as Problem #18, but with a much larger triangle. Thus, the dynamic approach that was presented in my blog post for Problem #18 is mandatory for solving this problem. Here is my solution:

Solution #1: Dynamic Approach (Copy/Paste from #18)

We simply store the maximum path sum from each cell to the bottom and work our way up the triangle as we did in Problem #18. Here is an implementation of this approach in Python 2.7:

 '''
 Author: Walker Kroubalkian
 Dynamic Approach to Project Euler Problem #67
 '''
 
 import time
 
 f = open("PE67Triangle.txt","r")
 
 if(f.mode == "r"):
     contents = f.readlines()
     realContents = []
     for x in contents:
         realContents.append(map(int,x.split()))
 else:
     raise ValueError("Cannot read from file")
 
 def projectEulerProblemSixtySeven(triangle):
     rows = len(triangle)
     currentRow = triangle[rows-1]
     c = (rows-2)
     while(c>=0):
         temp = []
         for i in range(c+1):
             temp.append(triangle[c][i] + max(currentRow[i],currentRow[i+1]))
         currentRow = temp
         c-=1
     return currentRow[0]
 
 start = time.time()
 print projectEulerProblemSixtySeven(realContents)
 print ("--- %s seconds ---" % (time.time()-start))
 
 '''
 Prints
 
 7273
 --- 0.00120997428894 seconds ---
 
 for input of triangle = given triangle
 '''

And with that, we are officially 2/3 of the way to solving 100 questions on Project Euler. This is a great example of the powers of Dynamic programming. Rather than brute forcing 2^100 paths, we only need to run through the 5050 cells in the triangle once, and the program runs in just over a millisecond in Python 2.7.

Thanks for reading! See you tomorrow.

Project Euler Problem #66

Problem #66 concerns the minimal solutions to Pell Equations. The question reads:

Project Euler Problem 66: Diophantine equation
Consider quadratic Diophantine equations of the form:
x² – Dy² = 1
For example, when D=13, the minimal solution in x is 649² – 13×180² = 1.
It can be assumed that there are no solutions in positive integers when D is square.
By finding minimal solutions in x for D = {2, 3, 5, 6, 7}, we obtain the following:
3² – 2×2² = 1
2² – 3×1² = 1
9² – 5×4² = 1
5² – 6×2² = 1
8² – 7×3² = 1
Hence, by considering minimal solutions in x for D ≤ 7, the largest x is obtained when D=5.
Find the value of D ≤ 1000 in minimal solutions of x for which the largest value of x is obtained.

My solution for this problem is a bit risky because it relies on arbitrary precision. It uses the useful fact that the solutions to the Pell Equation are related to the continued fraction expansions for the square root of D. Here is my solution:

Solution #1: Arbitrary Precision Approach

We simply use the Python decimal library to obtain very precise continued fractions of the irrational square roots less than 1000. Then we keep calculating continued fractions until one provides a solution to the Pell Equation. Here is an implementation of this approach in Python 2.7:

 '''
 Author: Walker Kroubalkian
 Brute Force Approach to Project Euler Problem #66
 '''
 
 import time
 from math import sqrt
 from decimal import *
 
 getcontext().prec = 74
 
 def isSquare(n):
     a = int(sqrt(n))
     return a*a==n
 
 def gcd(a,b):
     if(a<0 or b<0):
         return gcd(abs(a),abs(b))
     if(min(a,b)==0):
         return max(a,b)
     if(a>b):
         return gcd(b,a%b)
     return gcd(a,b%a)
 
 def continuedFractionDenominators(n,x):
     a = Decimal(n).sqrt()
     denominators = []
     while(len(denominators)<x):
         v = int(a)
         denominators.append(v)
         a = Decimal(1.)/(a-Decimal(v))
     return denominators
 
 def evaluateContinuedFraction(n,x):
     a = continuedFractionDenominators(n, x)
     b = a[::-1]
     curNumerator = 1
     curDenominator = 0
     for x in b:
         n = curDenominator + curNumerator*x
         d = curNumerator
         g = gcd(n,d)
         curNumerator = n/g
         curDenominator = d/g
     return [curNumerator,curDenominator]
 
 def projectEulerProblemSixtySix(n):
     maxFound = 3
     maxIndex = 2
     for c in range(3,n+1):
         if(not isSquare(c)):
             d = 2
             while(True):
                 e = evaluateContinuedFraction(c, d)
                 f = e[0]
                 g = e[1]
                 if(f*f-c*g*g==1):
                     if(f>maxFound):
                         maxFound = f
                         maxIndex = c
                     break
                 d+=1
     return maxIndex
 
 start = time.time()
 print projectEulerProblemSixtySix(1000)
 print ("--- %s seconds ---" % (time.time()-start))
 
 '''
 Prints
 
 661
 --- 4.85595107079 seconds ---
 
 for input of n = 1000
 '''

And with that, we’re done. I needed 74 digits after the decimal place to get an accurate answer, so the program took over 4 seconds to run. I’d like to come back to this problem eventually to look for a more efficient and rigorous method.

Thanks for reading! See you tomorrow.

Project Euler Problem #65

Problem #65 concerns the infinite continued fraction for e. The question reads:

Once again, I apologize for not embedding this problem in WordPress. The problem had a ton of LaTeX this time and I thought it would be easier to just upload a screenshot. Unfortunately, the continued fractions for e cannot be reduced to solving a Pell Equation like the continued fractions for irrational square roots. Therefore, my solution uses brute force to get an answer.

Solution #1: Brute Force Approach

We simply evaluate the exact 100th expansion of the continued fraction for e. In case you cannot tell, the continued fraction for e can be written in the form e = [1,0,1,1,2,1,1,4,1,1,6,1,1,8,1,1,10, … 1,1,2k, …]. Here is a proof of this continued fraction: Link. Thus, the 100th continued fraction is composed of the first 101 terms of this array. Therefore, all we have to do is list the first 101 terms of this array and work backwards, simplifying any reducible fractions we encounter. Here is an implementation of this approach in Python 2.7:

 '''
 Author: Walker Kroubalkian
 Brute Force Approach to Project Euler Problem #65
 '''
 
 import time
 
 def gcd(a,b):
     if(min(a,b)<0):
         return gcd(abs(a),abs(b))
     if(min(a,b)==0):
         return max(a,b)
     if(a>b):
         return gcd(b,a%b)
     return gcd(a,b%a)
 
 def projectEulerProblemSixtyFive(n):
     denominatorList = [1,0]
     for c in range(n):
         if(c%3==2):
             denominatorList.append((c+1)*2/3)
         else:
             denominatorList.append(1)
     curNumerator = 1
     curDenominator = 0
     order = denominatorList[::-1]
     for x in order:
         n = curDenominator+curNumerator*x
         d = curNumerator
         g = gcd(n,d)
         if(g==0):
             g = 1
         curNumerator = n/g
         curDenominator = d/g
     total = 0
     for x in str(curNumerator):
         total+=int(x)
     return total
 
 start = time.time()
 print projectEulerProblemSixtyFive(100)
 print ("--- %s seconds ---" % (time.time() - start))
 
 '''
 Prints
 
 272
 --- 0.00275802612305 seconds ---
 
 for input of n = 100.
 '''

For context, the exact numerator of this expansion is 6963524437876961749120273824619538346438023188214475670667. Thus, while this approach may work fine in Python, it will probably require some clever implementation in languages that cannot store numbers this large. Regardless, this method runs quite quickly in Python 2.7.

Thanks for reading! See you tomorrow.

Project Euler Problem #64

Problem #64 is another Project Euler problem that concerns infinite continued fractions for irrational numbers. The question reads:

I apologize for not embedding the wording of the problem within WordPress. I was a bit lazy this time because the problem was so long. This is a good example of when a coding problem gives too many examples to build an understanding. Regardless, my solution for this problem involves a bit of arbitrary brute force with the “decimal” library in Python. Here’s my solution:

Solution #1: Arbitrary Precision Approach

We simply use the Python decimal library to get very precise values of the irrational square roots of every number less than 10,000 and evaluate the continued fraction for those precise square roots until a loop is found. This is a very risky approach because in the event that the value is not precise enough, the program could run into an infinite loop for the continued fraction or the values within the continued fraction could get too big for the program to store them accurately. As a result, I had to guess and check what level of precision was needed. Here is an implementation of this approach in Python 2.7 with the decimal library:

 '''
 Author: Walker Kroubalkian
 Brute Force Approach to Project Euler Problem #64
 '''
 
 import time
 from math import floor
 from decimal import *
 
 getcontext().prec = 214
 
 def gcd(a,b):
     if(a<0 or b<0):
         return gcd(abs(a),abs(b))
     if(min(a,b)==0):
         return max(a,b)
     if(a>b):
         return gcd(b,a%b)
     return gcd(a,b%a)
 
 def projectEulerProblemSixtyFour(n):
     total = 0
     for e in range(2,n+1):
         v = Decimal(e).sqrt()
         if(int(v)*int(v)!=e):  
             fractions = []
             root = [1,1]
             constant = [0,1]
             denominators = []
             while [root,constant] not in fractions:
                 fractions.append([root,constant])
                 a = root[0]
                 b = root[1]
                 c = constant[0]
                 d = constant[1]
                 x = int(floor(v))
                 denominators.append(x)
                 v = Decimal(1.)/(Decimal(v)-Decimal(int(v)))
                 
                 newRootN = a*b*d*d
                 newRootD = a*a*d*d*e-b*b*c*c-x*x*b*b*d*d+2*b*b*x*c*d
                 newConstantN = x*b*b*d*d-b*b*c*d
                 newConstantD = newRootD
                 
                 g1 = gcd(newRootN,newRootD)
                 g2 = gcd(newConstantN,newConstantD)
 

                 newRootN/=g1
                 newRootD/=g1
                 newConstantD/=g2
                 newConstantN/=g2
                 if(newRootD<0):
                     newRootN = -newRootN
                     newRootD = -newRootD
                 if(newConstantD<0):
                     newConstantD = -newConstantD
                     newConstantN = -newConstantN
                 root = [newRootN,newRootD]
                 constant = [newConstantN,newConstantD]
             period = len(fractions) - fractions.index([root,constant])
             if(period%2==1):
                 total+=1
             
     return total
 
 start = time.time()
 print projectEulerProblemSixtyFour(10000)
 print ("--- %s seconds ---" % (time.time()-start))
 
 '''
 Prints
 
 1322
 --- 12.4575679302 seconds ---
 
 for input of n = 10000.
 '''

And with that, we’re done. I needed 214 digits after the decimal point to precisely calculate the repeated continued fractions of all the square roots in this range. As a result, the program took over 12 seconds to run. I would definitely like to return to this problem to find an approach that does not use arbitrary precision, but it is good to know about Python libraries that allow for extra precision such as decimal.

Thanks for reading! See you tomorrow.

Project Euler Problem #63

Problem #63 concerns perfect powers and the number of digits in the perfect power. The question reads:

Project Euler Problem 63: Powerful digit counts
The 5-digit number, 16807=7⁵, is also a fifth power. Similarly, the 9-digit number, 134217728=8⁹, is a ninth power.
How many n-digit positive integers exist which are also an nth power?

My solution is heavily reliant on mathematics. Here is my solution:

Solution #1: Brute Force Approach

We know that any n-digit number must lie between 10^(n-1) and 10^n. Thus, if 9^n is not in this range, there will not be any nth powers with n digits. We can observe that because 10^(n-1) grows faster than 9^n, if 9^c < 10^(c-1) for some value of c, then 9^x < 10^(x-1) for x>c. Thus, it suffices to count the number of nth powers with this property for all powers of 10 less than 10^c. Here is an implementation of this approach in Python 2.7:

 '''
 Author: Walker Kroubalkian
 Brute Force Approach to Project Euler Problem #63
 '''
 
 import time
 from math import ceil, floor
 
 def projectEulerProblemSixtyThree(n):
     total = 0
     for c in range(1,n+1):
         first = ceil(((10.)**(float(c)-1.))**(1./float(c)))
         second = floor(((10.)**(float(c))-1.)**(1./float(c)))
         first = int(first)
         second = int(second)
         
         for x in range(first,second+1):
             if(len(str(x**c))==c):
                 total+=1
         
         if(9**c < 10**(c-1)):
             return total
     return total
 
 start = time.time()
 print projectEulerProblemSixtyThree(1000)
 print ("--- %s seconds ---" % (time.time()-start))
 
 '''
 Prints
 
 49
 --- 9.89437103271e-05 seconds ---
 
 for input of n = 1000.
 '''

And with that, we’re done. The problems on Project Euler that involve pure mathematics such as inequalities and geometry tend to have very quick solutions, so it is no surprise that this one runs in under a millisecond.

Thanks for reading! See you tomorrow.

Project Euler Problem #62

Problem #62 concerns perfect cubes whose digits are rearrangements of each other. The question reads:

Project Euler Problem 62: Cubic permutations
The cube, 41063625 (345³), can be permuted to produce two other cubes: 56623104 (384³) and 66430125 (405³). In fact, 41063625 is the smallest cube which has exactly three permutations of its digits which are also cube.
Find the smallest cube for which exactly five permutations of its digits are cube.

My solution for this problem falls into the category of brute force. It takes over 20 seconds to run, so I’d like to look for a more efficient solution to this question at some point. Here is my solution:

Solution #1: Brute Force Approach

We simply list all perfect cubes of each number of digits until we find a quintuple of cubes which are all rearrangements of the digits in the other cubes. As you might imagine, this process can take a while. To determine if 2 numbers are permutations of each other, I sorted the list of digits in each and checked if the results were equal. Here is an implementation of this method in Python 2.7:

 '''
 Author: Walker Kroubalkian
 Brute Force Approach to Project Euler Problem #62
 '''
 
 import time
 
 def isPermutation(a,b):
     c = str(a)
     d = str(b)
     if(len(c)!=len(d)):
         return False
     return sorted(c) == sorted(d)
 
 def projectEulerProblemSixtyTwo(n):
     c = 345
     found = 0
     current = -1
     cubes = []
     currCube = c**3
     while(found<n):
         c+=1
         v = c**3
         if(len(str(v))>len(str(currCube))):
             cubes = []
         currCube = v
         t = 0
         first = -1
         for a in cubes:
             if isPermutation(a,currCube):
                 if(first==-1):
                     first = a
                 t+=1
         if(t==4):
             found+=1
             current = first
         cubes.append(currCube)
     return current
 

 start = time.time()
 print projectEulerProblemSixtyTwo(1)
 print ("--- %s seconds ---" % (time.time()-start))
 
 '''
 Prints
 
 127035954683
 --- 21.7080910206 seconds ---

 for input of n = 1.
 '''

Yikes. That took over 20 seconds to run. Also, looking back on it, I never checked that the number has exactly five rearrangements. This code will only find numbers with at least five rearrangements. I would like to return to this problem to search for a more efficient solution.

Thanks for reading! See you tomorrow.

Project Euler Problem #61

Problem #61 involves the figurate numbers. The question reads:

Project Euler Problem 61: Cyclical figurate numbers
Triangle, square, pentagonal, hexagonal, heptagonal, and octagonal numbers are all figurate (polygonal) numbers and are generated by the following formulae:
Triangle
 
P_3,n=n(n+1)/2
 
1, 3, 6, 10, 15, ...
Square
 
P_4,n=n²
 
1, 4, 9, 16, 25, ...
Pentagonal
 
P_5,n=n(3n−1)/2
 
1, 5, 12, 22, 35, ...
Hexagonal
 
P_6,n=n(2n−1)
 
1, 6, 15, 28, 45, ...
Heptagonal
 
P_7,n=n(5n−3)/2
 
1, 7, 18, 34, 55, ...
Octagonal
 
P_8,n=n(3n−2)
 
1, 8, 21, 40, 65, ...
The ordered set of three 4-digit numbers: 8128, 2882, 8281, has three interesting properties.
The set is cyclic, in that the last two digits of each number is the first two digits of the next number (including the last number with the first).
Each polygonal type: triangle (P_3,127=8128), square (P_4,91=8281), and pentagonal (P_5,44=2882), is represented by a different number in the set.
This is the only set of 4-digit numbers with this property.
Find the sum of the only ordered set of six cyclic 4-digit numbers for which each polygonal type: triangle, square, pentagonal, hexagonal, heptagonal, and octagonal, is represented by a different number in the set.

This is one of the problems that horrified me for years of doing Project Euler problems because I really didn’t want to implement a solution to it. My solution definitely falls into the category of brute force, but it does run quite fast so I believe it is acceptable. Here is my solution:

Solution #1: Brute Force Approach

We simply list all of the triangular, square, pentagonal, hexagonal, heptagonal, and octagonal numbers less than 10000. Next we iterate through all of the octagonal starting numbers and find all of the figurate numbers which start with the last two digits of the corresponding octagonal. For each 1st figurate number, we list all 2nd figurate numbers with the same digits. We do this five times until we have obtained a 6-tuple of figurate numbers. Finally, we check all 6-tuples obtained in this process and store the maximum total of the 6-tuples that satisfy the property. Like I said, this solution was meant solely to solve the problem no matter the cost, so it is far from elegant. Here is an implementation of it in Python 2.7:

 '''
 Author: Walker Kroubalkian
 Brute Force Approach to Project Euler Problem #61
 '''
 
 import time
 
 def projectEulerProblemSixtyOne():
     triangles = []
     c = 2
     v = 1
     while(v<10000):
         if(v>=1000):
             triangles.append(v)
         v+=c
         c+=1
     squares = []
     c = 32
     v = c*c
     while(v<10000):
         squares.append(v)
         c+=1
         v += (2*c-1)
     pentagons = []
     c = 2
     v = 1
     while(v<10000):
         if(v>=1000):
             pentagons.append(v)
         v+=(3*c-2)
         c+=1
     hexagons = []
     c = 2
     v = 1
     while(v<10000):
         if(v>=1000):
             hexagons.append(v)
         v+=(4*c-3)
         c+=1
     heptagons = []
     c = 2
     v = 1
     while(v<10000):
         if(v>=1000):
             heptagons.append(v)
         v+=(5*c-4)
         c+=1
     octagons = []
     c = 2
     v = 1
     while(v<10000):
         if(v>=1000):
             octagons.append(v)
         v+=(6*c-5)
         c+=1
     total = triangles[:]
     for x in squares:
         if x not in total:
             total.append(x)
     for x in pentagons:
         if x not in total:
             total.append(x)
     for x in hexagons:
         if x not in total:
             total.append(x)
     for x in heptagons:
         if x not in total:
             total.append(x)
     for x in octagons:
         if x not in total:
             total.append(x)
     possibleTuples = []
     for o in octagons:
         a = o%100
         p1 = []
         for x in total:
             if(x/100==a):
                 p1.append(x)
         for o1 in p1:
             b = o1%100
             p2 = []
             for x in total:
                 if(x/100==b):
                     p2.append(x)
             for o2 in p2:
                 c = o2%100
                 p3 = []
                 for x in total:
                     if(x/100==c):
                         p3.append(x)
                 for o3 in p3:
                     d = o3%100
                     p4 = []
                     for x in total:
                         if(x/100==d):
                             p4.append(x)
                     for o4 in p4:
                         e = o4%100
                         p5 = []
                         for x in total:
                             if(x/100==e and x%100 == o/100):
                                 p5.append(x)
                         for o5 in p5:
                             possibleTuples.append([o,o1,o2,o3,o4,o5])
     realPossible = []
     for a in possibleTuples:
         whole = []
         hpRep = False
         heRep = False
         peRep = False
         sqRep = False
         trRep = False
         repeat = False
         for x in a:
             if x not in whole:
                 whole.append(x)
             else:
                 repeat = True
                 break
             if(x in hexagons):
                 heRep = True
             elif(x in triangles):
                 trRep = True
             elif(x in heptagons):
                 hpRep = True
             elif(x in pentagons):
                 peRep = True
             elif(x in squares):
                 sqRep = True
         if(hpRep and heRep and peRep and sqRep and trRep and not repeat):
             realPossible.append(a)
     maxFound = -1
     for x in realPossible:
         t = 0
         for y in x:
             t+=y
         if(t>maxFound):
             maxFound = t
     return maxFound
     
 start = time.time()
 print projectEulerProblemSixtyOne()
 print ("--- %s seconds ---" % (time.time()-start))
 
 '''
 Prints
 
 28684
 --- 0.0513989925385 seconds ---
 
 for given problem.
 '''

And with that, we’re done. 150 lines of code for a program that runs in a 20th of a second. Like I said in the intro, I dreaded solving this problem because it seemed like a pain to implement. I hope you can understand the messy solution I made.

Thanks for reading! See you tomorrow.

Project Euler Problem #60

Problem #60 is one of the many Project Euler problems that involves manipulating digits in prime numbers. The question reads:

Project Euler Problem 60: Prime pair sets
The primes 3, 7, 109, and 673, are quite remarkable. By taking any two primes and concatenating them in any order the result will always be prime. For example, taking 7 and 109, both 7109 and 1097 are prime. The sum of these four primes, 792, represents the lowest sum for a set of four primes with this property.

Find the lowest sum for a set of five primes for which any two primes concatenate to produce another prime.

My solution for this question barely doesn’t run in under 30 seconds, so I definitely need to retry this problem at some point. Here is my current solution to this problem:

Solution #1: Brute Force Approach

We simply iterate through the possible values of the largest prime in a quintuple, checking all pairs of primes for the concatenation property until 5 primes are found for which all pairs within the quintuple have the concatenation property. This can be done by manually checking each new maximum prime with all smaller primes and storing a list of all concatenation primes respective to each new prime. (Sorry if that’s confusing, I’m not sure how to word it) For all primes which concatenate with at least 4 smaller primes, we can list all 4-element subsets of the smaller primes which concatenate with it and then check all pairs within each of those subsets. Once a subset satisfies the condition for all pairs, we’re done. Here is an implementation of this approach in Python 2.7:

 '''
 Author: Walker Kroubalkian
 Brute Force Approach to Project Euler Problem #60
 '''
 
 import time
 from math import sqrt
 
 def isPrime(n):
     if(n%2==0):
         return False
     c = 3
     while(c<=sqrt(n)):
         if(n%c==0):
             return False
         c+=2
     return True
 
 def concatenatePrime(a,b):
     p1 = int(str(a)+str(b))
     if not isPrime(p1):
         return False
     return isPrime(int(str(b)+str(a)))
 
 def subsets(myList,n):
     a = myList[0]
     l = len(myList)
     if(n>l):
         return []
     if(n==l):
         return [myList]
     if(n==0):
         return [[]]
     first = subsets(myList[1:],n-1)
     total = subsets(myList[1:],n)
     for x in first:
         c = x[:]
         c.append(a)
         total.append(sorted(c))
     return total
 
 def projectEulerProblemSixty(n):
     consideredPrimes = []
     primeConcats = []
     c = 3
     found = 0
     indexCheck = subsets([0,1,2,3],2)
     current = [-1]
     minSum = -1
     while(found<n):
         if(isPrime(c)):
             if c not in consideredPrimes:
                 concatPossible = []
                 l = len(consideredPrimes)
                 for a in range(l):
                     if(concatenatePrime(consideredPrimes[a], c)):
                         concatPossible.append(a)
                         primeConcats[a].append(l)
                 if(len(concatPossible)>=4):       
                     toCheck = subsets(concatPossible,4)
                     for a in toCheck:
                         haveFound = False
                         for b in indexCheck:
                             x = b[0]
                             y = b[1]
                             if a[x] not in primeConcats[a[y]]:
                                 haveFound = True
                                 break
                         if not haveFound:
                             found+=1
                             e = [c]
                             for z in a:
                                 e.append(consideredPrimes[z])
                             subTotal = 0
                             for t in e:
                                 subTotal+=t
                             if(current!=[-1]):
                                 if(subTotal<minSum):
                                     current = e
                                     minSum = subTotal
                             else:
                                 current = e
                                 minSum = subTotal
                                 
                 consideredPrimes.append(c)
                 primeConcats.append(concatPossible)
         c+=2
     return minSum
 
 start = time.time()
 print projectEulerProblemSixty(1)
 print ("--- %s seconds ---" % (time.time()-start))
 
 '''
 Prints
 
 26033
 --- 30.2959289551 seconds ---
 
 for input of n = 1.
 '''

And with that we’re done. Looking back on this, I really need to add some comments and make my code more clean. It’s pretty hideous and it took me a while to figure out my own code when writing this post (I wrote this code several months ago) In addition, I need to look for a more efficient algorithm because this took over 30 seconds to run.

Thanks for reading! See you tomorrow.

Project Euler Problem #59

Problem #59 concerns a type of cipher in cryptography known as XOR encryption. The question reads:

Project Euler Problem 59: XOR decryption
Each character on a computer is assigned a unique code and the preferred standard is ASCII (American Standard Code for Information Interchange). For example, uppercase A = 65, asterisk (*) = 42, and lowercase k = 107.

A modern encryption method is to take a text file, convert the bytes to ASCII, then XOR each byte with a given value, taken from a secret key. The advantage with the XOR function is that using the same encryption key on the cipher text, restores the plain text; for example, 65 XOR 42 = 107, then 107 XOR 42 = 65.

For unbreakable encryption, the key is the same length as the plain text message, and the key is made up of random bytes. The user would keep the encrypted message and the encryption key in different locations, and without both "halves", it is impossible to decrypt the message.

Unfortunately, this method is impractical for most users, so the modified method is to use a password as a key. If the password is shorter than the message, which is likely, the key is repeated cyclically throughout the message. The balance for this method is using a sufficiently long password key for security, but short enough to be memorable.

Your task has been made easy, as the encryption key consists of three lower case characters. Using p059_cipher.txt (right click and 'Save Link/Target As...'), a file containing the encrypted ASCII codes, and the knowledge that the plain text must contain common English words, decrypt the message and find the sum of the ASCII values in the original text.

As it turns out, XOR Decryption is exactly the same as XOR encryption. Setting c = a XOR b, we can find that c XOR b = a for all a and b. This fact will be used to greatly simplify our brute force.

Solution #1: Brute Force Approach

It suffices to check all 26*26*26 = 17,576 possible 3-letter keys and store the ones that decrypt the cipher text into a reasonable deciphered text. It takes a bit of guessing and checking to determine what characters are acceptable in a reasonable deciphered text, but after going through the most common characters in typical English text, we can find that only one key provides a reasonable answer. Here is an implementation of this method in Python 2.7:

 '''
 Author: Walker Kroubalkian
 Basic Cryptography Approach to Project Euler Problem #59
 '''
 
 import time
 
 f = open("PE59CipherText.txt", "r")
 
 if(f.mode == "r"):
     contents = f.readlines()
     realContents = []
     for x in contents:
         realContents.append(x.split(","))
 else:
     raise ValueError("Cannot read from file")
 
 finalContents = []
 for x in realContents[0]:
     finalContents.append(int(x))
 
 def isAcceptable(myChar):
     return 'a'<=myChar<='z' or 'A'<=myChar<='Z' or myChar == ' ' or myChar == '\'' or myChar == ',' or myChar == '"' or myChar == '[' or myChar == ']' or myChar == ':' or '0'<=myChar<='9' or myChar == '/' or myChar == '.' or myChar == '(' or myChar == ')' or myChar == ';' or myChar=='+'
     
 def projectEulerProblemFiftyNine(cipherText):
     tl = len(cipherText)
     alphabet = "abcdefghijklmnopqrstuvwxyz"
     possible = []
     for a in alphabet:
         for b in alphabet:
             for c in alphabet:
                 myGuess = a+b+c
                 asciiIndices = [ord(a),ord(b),ord(c)]
                 s = ""
                 bad = False
                 for x in range(tl):
                     myChar = chr(cipherText[x] ^ asciiIndices[x%3])
                     if not isAcceptable(myChar):
                         bad = True
                         break
                     s+=myChar
                 if not bad:
                     possible.append(s)
     total = -1
     for x in possible:
         subTotal = 0
         for y in x:
             subTotal+=ord(y)
         if(total==-1):
             total = subTotal
     return total
 
 start = time.time()
 print projectEulerProblemFiftyNine(finalContents)
 print ("--- %s seconds ---" % (time.time()-start))
 
 '''
 Prints
 
 129448
 --- 0.169327974319 seconds ---
 
 for input of given cipher text.
 '''

And with that, we’re done. In case you’re wondering, the key was “exp” and the decrypted text was “An extract taken from the introduction of one of Euler’s most celebrated papers, “De summis serierum reciprocarum” [On the sums of series of reciprocals]: I have recently found, quite unexpectedly, an elegant expression for the entire sum of this series 1 + 1/4 + 1/9 + 1/16 + etc., which depends on the quadrature of the circle, so that if the true sum of this series is obtained, from it at once the quadrature of the circle follows. Namely, I have found that the sum of this series is a sixth part of the square of the perimeter of the circle whose diameter is 1; or by putting the sum of this series equal to s, it has the ratio sqrt(6) multiplied by s to 1 of the perimeter to the diameter. I will soon show that the sum of this series to be approximately 1.644934066842264364; and from multiplying this number by six, and then taking the square root, the number 3.141592653589793238 is indeed produced, which expresses the perimeter of a circle whose diameter is 1. Following again the same steps by which I had arrived at this sum, I have discovered that the sum of the series 1 + 1/16 + 1/81 + 1/256 + 1/625 + etc. also depends on the quadrature of the circle. Namely, the sum of this multiplied by 90 gives the biquadrate (fourth power) of the circumference of the perimeter of a circle whose diameter is 1. And by similar reasoning I have likewise been able to determine the sums of the subsequent series in which the exponents are even numbers.”

Thanks for reading! See you tomorrow.

Project Euler Problem #58

Problem #58 is another example of a Project Euler problem where the process of Engineer’s Induction can be very useful. The question reads:

Project Euler Problem 58: Spiral primes
Starting with 1 and spiralling anticlockwise in the following way, a square spiral with side length 7 is formed.

37 36 35 34 33 32 31
38 17 16 15 14 13 30
39 18  5  4  3 12 29
40 19  6  1  2 11 28
41 20  7  8  9 10 27
42 21 22 23 24 25 26
43 44 45 46 47 48 49

It is interesting to note that the odd squares lie along the bottom right diagonal, but what is more interesting is that 8 out of the 13 numbers lying along both diagonals are prime; that is, a ratio of 8/13 ≈ 62%.

If one complete new layer is wrapped around the spiral above, a square spiral with side length 9 will be formed. If this process is continued, what is the side length of the square spiral for which the ratio of primes along both diagonals first falls below 10%?

My solution for this problem definitely falls into the category of brute force. Here is my solution:

Solution #1: Engineer’s Induction and Brute Force Approach

We can notice by Engineer’s Induction that the four numbers on the diagonals in the (c-1)/2-th layer of the spiral have values c*c+c+1, c*c+2*c+2, c*c+3*c+3, and c*c+4*c+4. By checking whether each of these numbers is prime through the conventional method, the proportion of primes on the diagonals of the spiral can be determined. Here is an implementation of this method in Python 2.7:

'''
 Author: Walker Kroubalkian
 Engineer's Induction Approach to Project Euler Problem #58
 '''

 import time
 from math import sqrt

 def isPrime(n):
     if(n<=1):
         return False
     if(n<4):
         return True
     if(n%2==0):
         return False
     c = 3
     while(c<=sqrt(n)):
         if(n%c==0):
             return False
         c+=2
     return True

 def projectEulerProblemFiftyEight(p):
     total = 5
     found = 3
     c = 3
     while(found100>=totalp):
         t = c*c+c+1
         for a in range(3):
             if isPrime(t):
                 found+=1
             total+=1
             t+=(c+1)
         total+=1
         c+=2
     return c

 start = time.time()
 print projectEulerProblemFiftyEight(10)
 print ("--- %s seconds ---" % (time.time()-start))

 '''
 Prints
 26241
 --- 3.82181882858 seconds ---
 for input of p = 10
 '''

And with that, we’re done. This method took 3.8 seconds to run in Python 2.7. Looking back on it, it could be improved by 25% by noticing that c*c+4*c+4 is a perfect square for all integer values of c. There are probably many other ways to make this program more efficient, but unfortunately I am too lazy to look for them at this moment. I may come back to this problem eventually.

Thanks for reading! See you tomorrow.