You are currently browsing the category archive for the ‘computing’ category.
Generating e through probability and hypercubes
This is a really beautiful solution to an interesting probability problem posed by fellow IB teacher Daniel Hwang, for which I’ve outlined a method for solving suggested by Ferenc Beleznay. The problem is as follows:
On average, how many random real numbers from 0 to 1 (inclusive) are required for the sum to exceed 1?
1 number
Clearly if we choose only 1 number then we can’t exceed 1.
2 numbers
Here we imagine the 2 numbers we pick as x and y and therefore we can represent them as a coordinate pair. The smallest pair (0,0) and the largest pair (1,1). This means that the possible coordinates fit inside the unit square shown above. We want to know for what coordinate pairs we have the inequality x + y > 1. This can be rearrange to give y > 1-x. The line y = 1-x is plotted and we can see that any coordinate points in the triangle BCD satisfy this inequality. Therefore the probability of a random coordinate pair being in this triangle is 1/2.
3 numbers
This time we want to find the probability that we exceed 1 with our third number. We can consider the numbers as x, y, z and therefore as 3D coordinates (x,y,z). From the fact that we are choosing a third number we must already have x +y <1. We draw the line x+y = 1, which in 3D gives us a plane. The volume in which our coordinate point must lie is the prism ABDEFG.
We now also add the constraint x+y+z >1. This creates the plane as shown. If our coordinate lies inside the pyramid ABDE then our coordinates will add to less than 1, outside this they will add to more than 1.
The volume of the pyramid ABDE = 1/3 (base area)(perpendicular height).
The volume of the prism ABDEFG = (base area)(perpendicular height).
Given that they share the same perpendicular height and base area then precisely 1/3 of the available volume would give a coordinate point that adds to less than 1, and 2/3 of the available volume would give a coordinate point that adds to more than 1.
Therefore we have the following tree diagram:
Exceeds 1 with 2 numbers = 1/2
Does not exceed 1 with 2 numbers, exceeds 1 with 3 numbers = 1/2 x 2/3 = 1/3.
Does not exceed 1 with 2 numbers, does not exceed 1 with 3 numbers = 1/2 x 1/3 = 1/6.
4 numbers
If you been following so far this is where things get interesting! We can now imagine a 4 dimensional unit cube (image above from Wikipedia) and a 4D coordinate point (x,y,z,a).
Luckily all we care about is the ratio of the 4-D pyramid and the 4-D prim formed by our constraints x+y+z <1 and x+y+z+a >1.
We have the following formula to help:
The n-D volume of a n-D pyramid = 1/n (base)(perpendicular height).
Therefore:
The 4-D volume of a 4-D pyramid = 1/4 (base 3D volume)(perpendicular height).
The 4-D volume of the prism ABDEFG = (base 3D volume)(perpendicular height).
Given that the 2 shapes share the same base and perpendicular height, the hyper-pyramid occupies exactly 1/4 of the 4-D space of the hyper-prism. So the probability of being in this space is 1/4 and 3/4 of being outside this space.
We can now extend our tree diagram:
Does not exceed 1 with 2 numbers, does not exceed 1 with 3 numbers, exceeds with 4 numbers = 1/2 x 1/3 x 3/4 = 1/8
Does not exceed 1 with 2 numbers, does not exceed 1 with 3 numbers, does not exceed with 4 numbers = 1/2 x 1/3 x 1/4 = 1/24.
In general a hyper-pyramid in n dimensional space occupies exactly 1/n of the space of the hyper-prism – so we can now continue this tree diagram.
Expected value
We can make a table of probabilities to find how many numbers we expect to use in order to exceed one.
Which gives us the following expected value calculation:
Which we can rewrite as:
But we have:
Therefore this gives:
So on average we would need to pick e numbers for the sum to exceed one! This is quite a remarkable result – e, one of the fundamental mathematical constants has appeared as if by magic on a probability question utilizing hyper-dimensional shapes.
Demonstrating this with Python
Running the Python code shown above will simulate doing this experiment. The computer generates a “random” number, then another and carries on until the sum is greater than 1. It then records how many numbers were required. It then does this again 1 million times and finds the average from all the trials.
1 million simulations gives 2.7177797177797176. When we compare this with the real answer for e, 2.7182818284590452353602874713527, we can see it has taken 1 million simulations to only be correct to 4sf.
Even 5 million simulations only gives 2.7182589436517888, so whilst we can clearly see that we will eventually get e, it’s converging very slowly. This may be because we are reliant on a random number generator which is not truly random (and only chooses numbers to a maximum number of decimal places rather than choosing from all values between 0 and 1).
I think this is a beautiful example of the unexpected nature of mathematics – we started out with a probability problem and ended up with e, via a detour into higher dimensional space! We can also see the power of computers in doing these kinds of brute force calculations.
Have you got a Super Brain?
Adapting and exploring maths challenge problems is an excellent way of finding ideas for IB maths explorations and extended essays. This problem is taken from the book: The first 25 years of the Superbrain challenges. I’m going to see how many different ways I can solve it.
The problem is to find all the integer solutions to the equation above. Finding only integer solutions is a fundamental part of number theory – a branch of mathematics that only deals with integers.
Method number 1: Brute force
This is a problem that computers can make short work of. Above I wrote a very simple Python program which checked all values of x and y between -99 and 99. This returned the only solution pairs as:
Clearly we have not proved these are the only solutions – but even by modifying the code to check more numbers, no more pairs were found.
Method number 2: Solving a linear equation
We can notice that the equation is linear in terms of y, and so rearrange to make y the subject.
We can then use either polynomial long division or the method of partial fractions to rewrite this. I’ll use partial fractions. The general form for this fraction can be written as follows:
Next I multiply by the denominator and the compare coefficients of terms.
This therefore gives:
I can now see that there will only be an integer solution for y when the denominator of the fraction is a factor of 6. This then gives (ignoring non integer solutions):
I can then substitute these back to find my y values, which give me the same 4 coordinate pairs as before:
Method number 3: Solving a quadratic equation
I start by making a quadratic in x:
I can then use the quadratic formula to find solutions:
Which I can simplify to give:
Next I can note that x will only be an integer solution if the expression inside the square root is a square number. Therefore I have:
Next I can solve a new quadratic as follows:
As before I notice that the expression inside my square root must be a square number. Now I can see that I need to find m and n such that I have 2 square numbers with a difference of 24. I can look at the first 13 square numbers to see that from the 12th and 13th square numbers onwards there will also be a difference of more than 24. Checking this list I can find that m = 1 and m = 5 will satisfy this equation.
This then gives:
which when I solve for integer solutions and then sub back into find x gives the same four solutions:
Method number 4: Graphical understanding
Without rearranging I could imagine this as a 3D problem by plotting the 2 equations:
This gives the following graph:
We can see that the plane intersects the curve in infinite places. I’ve marked A, B on the graph to illustrate 2 of the coordinate pairs which we have found. This is a nice visualization but doesn’t help find our coordinates, so lets switch to 2D.
In 2D we can use our rearranged equation:
This gives the following graph:
Here I have marked on the solution pairs that we found. The oblique asymptote (red) is y = 2x-1 because as x gets large the fraction gets very small and so the graph gets closer and closer to y = 2x -1.
All points on this curve are solutions to the equation – but we can see that the only integer solution pairs will be when x is small. When x is a large integer then the curve will be close to the asymptote and hence will return a number slightly bigger than an integer.
So, using this approach we would check all possible integer solutions when x is small, and again should be able to arrive at our coordinate pairs.
So, 4 different approaches that would be able to solve this problem. Can you find any others?
Sierpinski Triangle: A picture of infinity
This pattern of a Sierpinski triangle pictured above was generated by a simple iterative program. I made it by modifying the code previously used to plot the Barnsley Fern. You can run the code I used on repl.it. What we are seeing is the result of 30,000 iterations of a simple algorithm. The algorithm is as follows:
Transformation 1:
xi+1 = 0.5xi
yi+1= 0.5yi
Transformation 2:
xi+1 = 0.5xi + 0.5
yi+1= 0.5yi+0.5
Transformation 3:
xi+1 = 0.5xi +1
yi+1= 0.5yi
So, I start with (0,0) and then use a random number generator to decide which transformation to use. I can run a generator from 1-3 and assign 1 for transformation 1, 2 for transformation 2, and 3 for transformation 3. Say I generate the number 2 – therefore I will apply transformation 2.
xi+1 = 0.5(0) + 0.5
yi+1= 0.5(0)+0.5
and my new coordinate is (0.5,0.5). I mark this on my graph.
I then repeat this process – say this time I generate the number 3. This tells me to do transformation 3. So:
xi+1 = 0.5(0.5) +1
yi+1= 0.5(0.5)
and my new coordinate is (1.25, 0.25). I mark this on my graph and carry on again. The graph above was generated with 30,000 iterations.
Altering the algorithm
We can alter the algorithm so that we replace all the 0.5 coefficients of x and y with another number, a.
a = 0.3 has disconnected triangles:
When a = 0.7 we still have a triangle:
By a = 0.9 the triangle is starting to degenerate
By a = 0.99 we start to see the emergence of a line “tail”
By a = 0.999 we see the line dominate.
And when a = 1 we then get a straight line:
When a is greater than 1 the coordinates quickly become extremely large and so the scale required to plot points means the disconnected points are not visible.
If I alternatively alter transformations 2 and 3 so that I add b for transformation 2 and 2b for transformation 3 (rather than 0.5 and 1 respectively) then we can see we simply change the scale of the triangle.
When b = 10 we can see the triangle width is now 40 (we changed b from 0.5 to 10 and so made the triangle 20 times bigger in length):
Fractal mathematics
This triangle is an example of a self-similar pattern – i.e one which will look the same at different scales. You could zoom into a detailed picture and see the same patterns repeating. Amazingly fractal patterns don’t fit into our usual understanding of 1 dimensional, 2 dimensional, 3 dimensional space. Fractals can instead be thought of as having fractional dimensions.
The Hausdorff dimension is a measure of the “roughness” or “crinkley-ness” of a fractal. It’s given by the formula:
D = log(N)/log(S)
For the Sierpinski triangle, doubling the size (i.e S = 2), creates 3 copies of itself (i.e N =3)
This gives:
D = log(3)/log(2)
Which gives a fractal dimension of about 1.59. This means it has a higher dimension than a line, but a lower dimension than a 2 dimensional shape.
Sphere packing problem: Pyramid design
Sphere packing problems are a maths problems which have been considered over many centuries – they concern the optimal way of packing spheres so that the wasted space is minimised. You can achieve an average packing density of around 74% when you stack many spheres together, but today I want to explore the packing density of 4 spheres (pictured above) enclosed in a pyramid.
Considering 2 dimensions
First I’m going to consider the 2D cross section of the base 3 spheres. Each sphere will have a radius of 1. I will choose A so that it is at the origin. Using some basic Pythagoras this will give the following coordinates:
Finding the centre
Next I will stack my single sphere on top of these 3, with the centre of this sphere directly in the middle. Therefore I need to find the coordinate of D. I can use the fact that ABC is an equilateral triangle and so:
3D coordinates
Next I can convert my 2D coordinates into 3D coordinates. I define the centre of the 3 base circles to have 0 height, therefore I can add z coordinates of 0. E will be the coordinate point with the same x and y coordinates as D, but with a height, a, which I don’t yet know:
In order to find a I do a quick sketch, seen below:
Here I can see that I can find the length AD using trig, and then the height DE (which is my a value) using Pythagoras:
Drawing spheres
The general equation for spheres with centre coordinate (a,b,c) and radius 1 is:
Therefore the equation of my spheres are:
Plotting these on Geogebra gives:
Drawing a pyramid
Next I want to try to draw a pyramid such that it encloses the spheres. This is quite difficult to do algebraically – so I’ll use some technology and a bit of trial and error.
First I look at creating a base for my pyramid. I’ll try and construct an equilateral triangle which is a tangent to the spheres:
This gives me an equilateral triangle with lengths 5.54. I can then find the coordinate points of F,G,H and plot them in 3D. I’ll choose point E so that it remains in the middle of the shape, and also has a height of 5.54 from the base. This gives the following:
As we can see, this pyramid does not enclose the spheres fully. So, let’s try again, this time making the base a little bit larger than the 3 spheres:
This gives me an equilateral triangle with lengths 6.6. Taking the height of the pyramid to also be 6.6 gives the following shape:
This time we can see that it fully encloses the spheres. So, let’s find the density of this packing. We have:
Therefore this gives:
and we also have:
Therefore the density of our packaging is:
Given our diagram this looks about right – we are only filling less than half of the available volume with our spheres.
Comparison with real data
[Source: Minimizing the object dimensions in circle and sphere packing problems]
We can see that this task has been attempted before using computational power – the table above shows the average density for a variety of 2D and 3D shapes. The pyramid here was found to have a density of 46% – so our result of 44% looks pretty close to what we should be able to achieve. We could tweak our measurements to see if we could improve this density.
So, a nice mixture of geometry, graphical software, and trial and error gives us a nice result. You could explore the densities for other 2D and 3D shapes and see how close you get to the results in the table.
Square Triangular Numbers
Square triangular numbers are numbers which are both square numbers and also triangular numbers – i.e they can be arranged in a square or a triangle. The picture above (source: wikipedia) shows that 36 is both a square number and also a triangular number. The question is how many other square triangular numbers we can find?
The equation we are trying to solve is:
a2 = 0.5(b2+b)
for some a, b as positive integers. The LHS is the formula to generate square numbers and the RHS is the formula to generate the triangular numbers.
We can start with some simple Python code (which you can run here):
for c in range(1,10001):
for d in range(1,10001):
if c**2 == (d**2+d)/2:
print(c**2, c,d)
This checks the first 10000 square numbers and the first 10000 triangular numbers and returns the following:
1 1 1
36 6 8
1225 35 49
41616 204 288
1413721 1189 1681
48024900 6930 9800
i.e 1225 is the next square triangular number after 36, and can be formed as 352 or as 0.5(492+49). We can see that there are very few square triangular numbers to be found in the first 50 million numbers. The largest we found was 48,024,900 which is made by 69302 or as 0.5(98002+9800).
We can notice that the ratio between each consecutive pair of square triangular numbers looks like it converges as it gives:
36/1 = 36
1225/36 = 34.027778
41616/1225 = 33.972245
1413721/41616 = 33.970612
48024900/1413721 = 33.970564
So, let’s use this to predict that the next square triangular number will be around
48024900 x 33.9706 = 1,631,434,668.
If we square root this answer we get approximately 40391
If we solve 0.5(b2+b) = 1,631,434,668 using Wolfram we get approximately 57120.
Therefore let’s amend our code to look in this region:
for c in range(40380,40400):
for d in range(57100,57130):
if c**2 == (d**2+d)/2:
print(c**2, c,d)
This very quickly finds the next solution as:
1631432881 40391 57121
This is indeed 403912 – so our approximation was very accurate. We can see that this also gives a ratio of 1631432881/48024900 = 33.97056279 which we can then use to predict that the next term will be 33.970563 x 1631432881 = 55,420,693,460. Square rooting this gives a prediction that we will use the 235,416 square number. 235,4162 gives 55,420,693,056 (using Wolfram Alpha) and this is indeed the next square triangular number.
So, using a mixture of computer code and some pattern exploration we have found a method for finding the next square triangular numbers. Clearly we will quickly get some very large numbers – but as long as we have the computational power, this method should continue to work.
Using number theory
The ever industrious Euler actually found a formula for square triangular numbers in 1778 – a very long time before computers and calculators, so let’s have a look at his method:
We start with the initial problem, and our initial goal is to rearrange it into the following form:
Next we make a substitution:
Here, when we get to the equation 1 = x2 – 2y2 we have arrived at a Pell Equation (hence the rearrangement to get to this point). This particular Pell Equation has the solution quoted above where we can define Pk as
Therefore we have
Therefore for any given k we can find the kth square triangular number. The a value will give us the square number required and the b value will give us the triangular number required. For example with k = 3:
This tells us the 3rd square triangular number is the 35th square number or the 49th triangular number. Both these give us an answer of 1225 – which checking back from our table is the correct answer.
So, we have arrived at 2 possible methods for finding the square triangular numbers – one using modern computational power, and one using the skills of 18th century number theory.
When do 2 squares equal 2 cubes?
Following on from the hollow square investigation this time I will investigate what numbers can be written as both the sum of 2 squares, 2 cubes and 2 powers of 4. i.e a2+b2 = c3+d3 = e4+f4.
Geometrically we can think of this as trying to find an array of balls such that we can arrange them into 2 squares, or we can rearrange them and stack them to form 2 cubes, or indeed we can arrange them into 2 4-dimensional cubes. I’ll add the constraints that all of a,b,c,d,e,f should be greater than 1 and that the pair of squares or cubes (etc) must be distinct. Therefore we can’t for example have 2 squares the same size.
Infinite solutions
Let’s look at why we can easily find infinite solutions if the squares or cubes (etc) can be the same size.
We want to find solutions to:
a2+b2 = c3+d3 = e4+f4.
so we look at the powers 2,3,4 which have LCM of 12. Therefore if we choose powers with the same base we can find a solution. For example we chose to work with base 2. Therefore we choose
a = 26, b = 26, which gives 212+212
c = 24, d = 24, which gives 212+212
e = 23, f = 23, which gives 212+212
Clearly these will be the same. So we can choose any base we wish, and make the powers into the same multiples of 12 to find infinite solutions.
Writing some code
Here is some code that will find some other solutions:
list1=[]
for a in range(2, 200):
for b in range(2,200):
list1.append(a**2+b**2)
list2=[]
for j in list1:
for c in range(2,200):
for d in range(2,200):
if c**3+d**3 == j:
list2.append(c**3+d**3)
print(list2)
for k in list2:
for e in range(2,200):
for f in range(2,200):
if k == e**4+f**4:
print(k,e,f)
This returns the following solutions: 8192, 18737, 76832. Of these we reject the first as this is the solution 212+212 which we found earlier and which uses repeated values for the squares, cubes and powers of 4. The 3rd solution we also reject as this is formed by 14 4 + 14 4. Therefore the only solution up to 79202 (we checked every value up to and including 1992 + 1992) is:
18737 = 642+1212 = 173+243 = 114+84.
Therefore if we had 18,737 balls we could arrange them into 2 squares, a 64×64 square and a 121×121 square. Alternatively we could rearrange them into 2 cubes, one 17x17x17 and one 24x24x24. Or we could enter a higher dimensional space and create 2 tesseracts one with sides 11x11x11x11 and the other with 14x14x14x14.
With only 1 solution for around the first 80,000 numbers it looks like these numbers are quite rare – could you find another one? And could you find one that also satisfies g5+h5?
Hollow Cubes investigation
Hollow cubes like the picture above [reference] are an extension of the hollow squares investigation done previously. This time we can imagine a 3 dimensional stack of soldiers, and so try to work out which numbers of soldiers can be arranged into hollow cubes.
Therefore what we need to find is what numbers can be formed from a3-b3
Python code
We can write some Python3 code to find this out (this can be run here):
for k in range(1,200):
for a in range(0, 100):
for b in range(0,100):
if a**3-b**3 == k :
print(k,a,b)
This gives the following: (the first number is the number of soldiers and the 2 subsequent numbers are the 2 cubes).
1 1 0
7 2 1
8 2 0
19 3 2
26 3 1
27 3 0
37 4 3
56 4 2
61 5 4
63 4 1
64 4 0
91 6 5
98 5 3
117 5 2
124 5 1
125 5 0
127 7 6
152 6 4
169 8 7
189 6 3
We could perhaps investigate any patterns in these numbers, or explore how we can predict when a hollow cube has more than one solution. I’ll investigate which numbers can be written as both a hollow square and also a hollow cube.
Hollow squares and hollow cubes
list1=[]
for a in range(2, 50):
for b in range(2,50):
if a**2-b**2 !=0:
if a**2-b**2 > 0:
list1.append(a**2-b**2)
list2=[]
for j in list1:
for c in range(2,50):
for d in range(2,50):
if c**3-d**3 == j:
list2.append(c**3-d**3)
print(list2)
This returns the following numbers which can all be written as both hollow squares and hollow cubes.
[56, 91, 19, 117, 189, 56, 208, 189, 217, 37, 279, 152, 117, 448, 513, 504, 448, 504, 387, 665, 504, 208, 875, 819, 936, 817, 61, 999, 988, 448, 728, 513, 189, 1216, 936, 784, 335, 469, 1323, 819, 1512, 1352, 1197, 992, 296, 152, 1519, 1512, 1197, 657, 1664, 1323, 1647, 1736, 1701, 1664, 936, 504, 2107, 1387, 1216, 1027, 91, 2015, 279, 2232]
Hollow squares, cubes and hypercubes
Taking this further, can we find any number which can be written as a hollow square, hollow cube and hollow hypercube (4 dimensional cube)? This would require our soldiers to be able to be stretch out into a 4th dimensional space – but let’s see if it’s theoretically possible.
Here’s the extra code to type:
list1=[]
for a in range(2, 200):
for b in range(2,200):
if a**2-b**2 !=0:
if a**2-b**2 > 0:
list1.append(a**2-b**2)
list2=[]
for j in list1:
for c in range(2,200):
for d in range(2,200):
if c**3-d**3 == j:
list2.append(c**3-d**3)
print(list2)
for k in list2:
for e in range(2,200):
for f in range(2,200):
if k == e**4-f**4:
print(k)
Very pleasingly this does indeed find some solutions:
9919: Which can be formed as either 1002-92 or 223-93 or 104-34.
14625: Which can be formed as either 1212-42 or 253-103 or 114-24.
Given that these took some time to find, I think it’ll require a lot of computer power (or a better designed code) to find any number which is a hollow square, hollow cube, hollow hypercube and hollow 5-dimensional cube, but I would expect that there is a number out there that satisfies all criteria. Maybe you can find it?
The Van Eck Sequence
This is a nice sequence as discussed in the Numberphile video above. There are only 2 rules:
- If you have not seen the number in the sequence before, add a 0 to the sequence.
- If you have seen the number in the sequence before, count how long since you last saw it.
You start with a 0.
0
You have never seen a 0 before, so the next number is 0.
00
You have seen a 0 before, and it was 1 step ago, so the next number is 1.
001
You have never seen a 1 before, so the next number is 0.
0010
You have seen a 0 before, it was 2 steps ago, so the next number is 2.
00102.
etc.
I can run a quick Python program (adapted from the entry in the Online Encyclopedia of Integer Sequences here) to find the first 100 terms.
A181391 = [0, 0]
for n in range(1, 10**2):
for m in range(n-1, -1, -1):
if A181391[m] == A181391[n]:
A181391.append(n-m)
break
else:
A181391.append(0)
print(A181391)
This returns:
[0, 0, 1, 0, 2, 0, 2, 2, 1, 6, 0, 5, 0, 2, 6, 5, 4, 0, 5, 3, 0, 3, 2, 9, 0, 4, 9, 3, 6, 14, 0, 6, 3, 5, 15, 0, 5, 3, 5, 2, 17, 0, 6, 11, 0, 3, 8, 0, 3, 3, 1, 42, 0, 5, 15, 20, 0, 4, 32, 0, 3, 11, 18, 0, 4, 7, 0, 3, 7, 3, 2, 31, 0, 6, 31, 3, 6, 3, 2, 8, 33, 0, 9, 56, 0, 3, 8, 7, 19, 0, 5, 37, 0, 3, 8, 8, 1, 46, 0, 6, 23]
I then assigned each term an x coordinate value, i.e.:
0 , 0
1 , 0
2 , 1
3 , 0
4 , 2
5 , 0
6 , 2
7 , 2
8 , 1
9 , 6
10 , 0
11 , 5
12 , 0
13 , 2
14 , 6
15 , 5
16 , 4
17 , 0
18 , 5
19 , 3
20 , 0
etc.
This means that you can then plot the sequence as a line graph, with the y values corresponding to the sequence terms. As you can see, every time we hit a new peak the following value is 0, leading to the peaks and troughs seen below:
Let’s extend the sequence to the first 1000 terms:
We can see that the line y = x provides a reasonably good upper bound for this data:
But it is not known if every number would actually appear in the sequence somewhere – so this bound may not hold for larger values.
Length of steps before new numbers appear.
We can also investigate how long we have to wait to see each number for the first time by running the following Python code:
A181391 = [0, 0]
for n in range(1, 10**3):
for m in range(n-1, -1, -1):
if A181391[m] == A181391[n]:
A181391.append(n-m)
break
else:
A181391.append(0)
for m in range(1,50):
if A181391[n]==m:
print(m, ",", n+1)
break
This returns the following data:
1 , 3
2 , 5
6 , 10
5 , 12
4 , 17
3 , 20
9 , 24
14 , 30
15 , 35
17 , 41
11 , 44
8 , 47
42 , 52
20 , 56
32 , 59
18 , 63
7 , 66
31 , 72
33 , 81
19 , 89
etc.
The first coordinate tells us the number we are interested in, and the second number tells us how long we have to wait in the sequence until it appears. So (1 , 3) means that we have to wait until 3 terms in the sequence to see the number 1 for the first time.
Plotting this for numbers 1-50 on a graph returns the following:
So, we can see (for example that we wait 66 terms to first see a 7, and 173 terms to first see a 12. There seems to be a general trend that as the numbers get larger we have to wait longer to see them. Testing this with a linear regression we can see a weak to moderate correlation:
Checking for the numbers up to 300 we get the following:
For example this shows that we have to wait 9700 terms until we see the number 254 for the first time. Testing this with a linear correlation we have a weaker positive correlation than previously.
So, a nice and quick investigation using a combination of sequences, coding, graphing and regression, with lots of areas this could be developed further.
Computers can brute force a lot of simple mathematical problems, so I thought I’d try and write some code to solve some of them. In nearly all these cases there’s probably a more elegant way of coding the problem – but these all do the job! You can run all of these with a Python editor such as Repl.it. Just copy and paste the below code and see what happens.
1) Happy Numbers.
Happy numbers are defined by the rule that you start with any positive integer, square each of the digits then add them together. Now do the same with the new number. Happy numbers will eventually spiral down to a number of 1. Numbers that don’t eventually reach 1 are called unhappy numbers.
As an example, say we start with the number 23. Next we do 2²+3² = 13. Now, 1²+3² = 10. Now 1²+0² = 1. 23 is therefore a happy number.
k= int(input("type a 2 digit number "))
a = int(k%10)
c = int(k//100)
b = int(k//10 -10*c)
print (a**2+b**2+c**2)
for k in range (1,20):
k = a**2+b**2 + c**2
a = int(k%10)
c = int(k//100)
b = int(k//10 -10*c)
print (a**2+b**2+c**2)
2) Sum of 3 cubes
Most (though not all) numbers can be written as the sum of 3 cubes. For example:
13 + 23 + 23 = 17. Therefore 17 can be written as the sum of 3 cubes.
This program allows you to see all the combinations possible when using the integers -10 to 10 and trying to make all the numbers up to 29.
for k in range(1,30):
for a in range(-10, 10):
for b in range(-10,10):
for c in range (-10, 10):
if a**3+b**3+c**3 == k :
print(k,a,b,c)
3) Narcissistic Numbers
A 3 digit narcissistic number is defined as one which the sum of the cubes of its digits equal the original number. This program allows you to see all 3 digit narcissistic numbers.
for a in range (0,10):
for b in range(0, 10):
for c in range(0,10):
if a**3 + b**3 + c**3 ==100*a + 10*b + c:
print(int(100*a + 10*b + c))
4) Pythagorean triples
Pythagorean triples are integer solutions to Pythagoras’ Theorem. For example:
32 + 42 = 52 is an integer solution to Pythagoras’ Theorem.
This code allows you to find all integer solutions to Pythagoras’ Theorem for the numbers in the range you specify.
k = 100
for a in range(1, k):
for b in range(1,k):
for c in range (1, 2*k):
if a**2+b**2==c**2:
print(a,b,c)
5) Perfect Numbers
Perfect numbers are numbers whose proper factors (factors excluding the number itself) add to the number. This is easier to see with an example.
6 is a perfect number because its proper factors are 1,2,3 and 1+2+3 = 6
8 is not a perfect number because its proper factors are 1,2,4 and 1+2+4 = 7
Perfect numbers have been known about for about 2000 years – however they are exceptionally rare. The first 4 perfect numbers are 6, 28, 496, 8128. These were all known to the Greeks. The next perfect number wasn’t discovered until around 1500 years later – and not surprisingly as it’s 33,550,336.
The code below will find all the perfect numbers less than 10,000.
for n in range(1,10000):
list = []
for i in range (1,n):
if n%i ==0:
list.append(i)
if sum(list)==n:
print(n)
Friendly Numbers
Friendly numbers are numbers which share a relationship with other numbers. They require the use of σ(a) which is called the divisor function and means the addition of all the factors of a. For example σ(7) = 1 + 7 = 8 and σ(10) = 1 +2 +5 + 10 = 18.
Friendly numbers therefore satisfy:
σ(a)/a = σ(b)/b
As an example:
σ(6) / 6 = (1+2+3+6) / 6 = 2,
σ(28) / 28 = (1+2+4+7+14+28) / 28 = 2
σ(496)/496 = (1+2+4+8+16+31+62+124+248+496)/496 = 2
Therefore 28 and 6 are friendly numbers because they share a common relationship.
This code will help find some Friendly numbers (though these are very difficult to find, as we need to check against every other integer until we find a relationship).
The code below will find some Friendly numbers less than 200, and their friendly pair less than 5000:
for n in range(1,5000):
list = []
for i in range (1,n+1):
if n%i ==0:
list.append(i)
Result1 = sum(list)
for m in range(1,200):
list2 = []
for j in range (1,m+1):
if m%j ==0:
list2.append(j)
Result2 = sum(list2)
if Result2/m ==Result1/n:
if n != m:
print(n,m)
Hailstone numbers
Hailstone numbers are created by the following rules:
if n is even: divide by 2
if n is odd: times by 3 and add 1
We can then generate a sequence from any starting number. For example, starting with 10:
10, 5, 16, 8, 4, 2, 1, 4, 2, 1…
we can see that this sequence loops into an infinitely repeating 4,2,1 sequence. Trying another number, say 58:
58, 29, 88, 44, 22, 11, 34, 17, 52, 26, 13, 40, 20, 10, 5, 16, 8, 4, 2, 1, 4, 2, 1…
and we see the same loop of 4,2,1.
The question is, does every number end in this loop? Well, we don’t know. Every number mathematicians have checked do indeed lead to this loop, but that is not a proof. Perhaps there is a counter-example, we just haven’t found it yet.
Run the code below, and by changing the value of n you can see how quickly the number enters the 4,2,1 loop.
n = 300
for k in range(1,40):
if n%2 ==0:
print(n/2)
n =n/2
elif n%2 !=0:
print(3*n+1)
n =3*n+1
Generating the Golden ratio
The Golden ratio can be approximated by dividing any 2 successive terms of the Fibonacci sequence. As we divide ever larger successive terms we get a better approximation for the Golden ratio. This code returns successive terms of the Fibonacci sequence and the corresponding approximation for the Golden ratio.
a = 0
b = 1
print(a)
print(b)
for k in range(1,30):
a = a+b
b = a+b
print(a,b, b/a)
Partial sums
We can use programs to see if sums to infinity converge. For example with the sequence 1/n, if I add the terms together I get: 1/1 + 1/2 + 1/3 + 1/4…In this case the series (surprisingly) diverges. The code below shows that the sum of the sequence 1/n2 converges to a number (pi2/6).
list = []
for n in range(1,100):
n = 1/(n**2)
list.append(n)
print(sum(list))
Returning to 6174
This is a nice number trick. You take any 4 digit number, then rearrange the digits so that you make the largest number possible and also the smallest number possible. You then take away the smallest number from the largest number, and then start again. For example with the number 6785, the largest number we can make is 8765 and the smallest is 5678. So we do 8765 – 5678 = 3087. We then carry on with the same method. Eventually we will arrive at the number 6174!
k= int(input("type a 4 digit number "))
a = int(k%10)
d = int(k//1000)
c = int(k//100 - 10*d)
b = int(k//10 -10*c-100*d)
for n in range(1,10):
list = []
list = [a,b,c,d]
list.sort()
a = list[0]
d = list[3]
c = list[2]
b = list[1]
print(1000*d+100*c+10*b+a -1000*a-100*b-10*c-d)
k = int(1000*d+100*c+10*b+a -1000*a-100*b-10*c-d)
a = int(k%10)
d = int(k//1000)
c = int(k//100 - 10*d)
b = int(k//10 -10*c-100*d)
list = []
list = [a,b,c,d]
list.sort()
a = list[0]
d = list[3]
c = list[2]
b = list[1]
print(1000*d+100*c+10*b+a -1000*a-100*b-10*c-d)
Maximising the volume of a cuboid
If we take a cuboid of length n, and cut squares of size x from the corner, what value of x will give the maximum volume? This code will look at initial squares of size 10×10 up to 90×90 and find the value of x for each which give the maximum volume.
def compute():
list1=[]
k=6
z = int(0.5*a*10**k)
for x in range(1,z):
list1.append((10*a-2*x/10**(k-1))*(10*a-2*x/10**(k-1))*(x/10**(k-1)))
print("length of original side is, ", 10*a)
y= max(list1)
print("maximum volume is, ", max(list1))
q = list1.index(y)
print("length of square removed from corner is, ", (q+1)/10**(k-1))
for a in range(1,10):
print(compute())
Stacking cannonballs – solving maths with code
Numberphile have recently done a video looking at the maths behind stacking cannonballs – so in this post I’ll look at the code needed to solve this problem.
Triangular based pyramid.
A triangular based pyramid would have:
1 ball on the top layer
1 + 3 balls on the second layer
1 + 3 + 6 balls on the third layer
1 + 3 + 6 + 10 balls on the fourth layer.
Therefore a triangular based pyramid is based on the sum of the first n triangular numbers.
The formula for the triangular numbers is:
and the formula for the sum of the first n triangular numbers is:
We can simplify this by using the identity for the sum of the first n square numbers and also the identity for the sum of the first n natural numbers:
Therefore:
and the question we want to find out is whether there is triangular based pyramid with a certain number of cannonballs which can be rearranged into a triangular number i.e.:
here n and m can be any natural number. For example if we choose n = 3 and m = 4 we see that we have the following:
Therefore we can have a triangular pyramid of height 3, which has 10 cannonballs. There 10 cannonballs can then be rearranged into a triangular number.
Square based pyramids and above.
For a square based pyramid we would have:
1 ball on the top layer
1 + 4 balls on the second layer
1 + 4 + 9 balls on the third layer
1 + 4 + 9 + 16 balls on the fourth layer.
This is the sum of the first n square numbers. So the formula for the square numbers is:
and the sum of the first n square numbers is:
For a pentagonal based pyramid we have:
1 ball on the top layer
1 + 5 balls on the second layer
1 + 5 + 12 balls on the third layer
1 + 5 + 12 + 22 balls on the fourth layer.
This is the sum of the first n pentagonal numbers. So the formula for the pentagonal numbers is:
and the formula for the first n pentagonal numbers is:
For a hexagonal based pyramid we have:
The formula for the first n hexagonal numbers:
and the formula for the sum of the first n hexagonal numbers:
For a k-agon based pyramid we have
and the formula for the sum of the first n k-agon numbers:
Therefore the general case is to ask if a k-agonal pyramid can be rearranged into a k-agon number i.e:
Computers to the rescue
We can then use some coding to brute force some solutions by running through large numbers of integers and seeing if any values give a solution. Here is the Python code. Type it (taking care with the spacing) into a Python editor and you can run it yourself.
You can then change the k range to check larger k-agons and also change the range for a and b. Running this we can find the following. (The first number is the value of k, the second the height of a k-agonal pyramid, the third number the k-agon number and the last number the number of cannonballs used).
Solutions:
3 , 3 , 4 , 10
3 , 8 , 15 , 120
3 , 20 , 55 , 1540
3 , 34 , 119 , 7140
4 , 24 , 70 , 4900
6 , 11 , 22 , 946
8 , 10 , 19 , 1045
8 , 18 , 45 , 5985
10 , 5 , 7 , 175
11 , 25 , 73 , 23725
14 , 6 , 9 , 441
14 , 46 , 181 , 195661
17 , 73 , 361 , 975061
20 , 106 , 631 , 3578401
23 , 145 , 1009 , 10680265
26 , 190 , 1513 , 27453385
29 , 241 , 2161 , 63016921
30 , 17 , 41 , 23001
32 , 298 , 2971 , 132361021
35 , 361 , 3961 , 258815701
38 , 430 , 5149 , 477132085
41 , 204 , 1683 , 55202400
41 , 505 , 6553 , 837244045
43 , 33 , 110 , 245905
44 , 586 , 8191 , 1408778281
50 , 34 , 115 , 314755
88 , 15 , 34 , 48280
145, 162, 1191, 101337426
276, 26, 77, 801801)
322, 28, 86, 1169686
823, 113, 694, 197427385
2378, 103, 604, 432684460
31265, 259, 2407, 90525801730
For example we can see a graphical representation of this. When k is 6, we have a hexagonal pyramid with height 11 or the 22nd hexagonal number – both of which give a solution of 946. These are all the solutions I can find – can you find any others? Leave a comment below if you do find any others and I’ll add them to the list!
