Screen Shot 2019-05-27 at 9.06.57 AM

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  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 :

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:

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:
 if sum(list)==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:
 Result1 = sum(list)
 for m in range(1,200):
  list2 = []
  for j in range (1,m+1):
   if m%j ==0:
  Result2 = sum(list2)

  if Result2/m ==Result1/n:
   if 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:
  n =n/2
 elif n%2 !=0:
  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
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)

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]

 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]

 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():

 z = int(0.5*a*10**k)

 for x in range(1,z):
 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):