How do I solve this billion iteration sum? Project Euler problem 94, Python

Question

Question -- Problem 94, Project Euler -- Python -- Almost equilateral triangles

*It is easily proved that no equilateral triangle exists with integral length sides and integral area. However, the almost equilateral triangle 5-5-6 has an area of 12 square units.

We shall define an almost equilateral triangle to be a triangle for which two sides are equal and the third differs by no more than one unit.

Find the sum of the perimeters of all almost equilateral triangles with integral side lengths and area and whose perimeters do not exceed one billion (1,000,000,000).*

Output: 518408346

I have a very inefficient solution and it is as follows:

1. I tried to bruteforce(it took a considerable amount of time) (using a for loop)
2. Then I added in the condition to check whether it was a triangle (sum of two sides is greater than the third)
3. Next, I defined the values in a perimeter variable and also used simple geometry for the area variable
4. Next, to check whether they were valid, for the perimeter I used a type method to check whether it was int and since I couldn't do the same for area(since even if it was an integer, it would be part of the float class), I used modulo 1 == 0 to check if it was an integer
5. If it satisfied the condition, I added it to my final 'sum of perimeter' variable.

The code is as follows:

import time
start = time.time()
plst = 0
o = 100000
for i in range(1,o):
    j = i+1
    if 2*i > j and j > 0:
        p = 2*i + j
        a = 0.5 * j* ((i**2 - (j/2)**2)**0.5)
        # if p > 100:
        #     break
        # else:
        #     pass
        if type(p) == int and a % 1 == 0:
            print(i)
            print('\t', p)
            plst += p
        else:
            pass
    else:
        pass
for k in range(1,o):
    b = k-1
    if 2*k > b and b > 0:
        p = 2*k + b
        a = 0.5 * b* ((k**2 - (b/2)**2)**0.5)
        # if p > 100:
        #     break
        # else:
        #     pass
        if type(p) == int and a % 1 == 0:
            print(k)
            print('\t', p)
            plst += p
        else:
            pass
    else:
        pass
print(plst)
print(f'{time.time() - start} seconds')

Note:

1. Yes, I know my o value is 100000.

2. Yes, i commented out the condition to check whether the perimeter variable crossed a billion, as I was facing problems before reaching one billion.

3. Don't mind the commented code. I have two problems for which I require help In order of precedence

4. When I run the code, uptill 10000, it works fine and fast, but as soon as I jump up to 100000, i realise that there are some unnecesary values that are popping up in the output.(like 93686,93686,93687)

   A. I call it unnecesary because when I searched this up, one website showed a table of values that will finally result in the answer and 93686 was not a part of it.(This is the website, https://blog.dreamshire.com/project-euler-94-solution/)

   B. But I could not find a mistake in my program and I tried out the area of these 'unnecesary' numbers and surprisingly the area turned out to be an integer and that surprised me and this resulted in me getting a wrong answer

So could someone explain what the mistake is?

2. Could someone give me an efficient solution with a proper explanation in Python?

I use python 3.8

Answer 1

You can obtain the result without using square roots at all. In fact it can be obtained using only additions and multiplications:

If we define an almost equilateral triangle as (a,a,c) then c is either a+1 or a-1. The height of the (a,a,c) triangle can be computed using a right triangle formed of half the base (c/2) with a as the hypothenuse.

                                               Pythagorean
                        /|\                |\  Triangle  
                       / | \               | \   
 c = a +/- 1        a /  |  \ a            |  \ a 
                     /  h|   \           h |   \  
                    /    |    \            |    \
                   /_____|_____\           |_____\
                         c                   c/2
                     c/2   c/2 

For this to produce an integer surface, h and c/2 must be integers and (h,c/2,a) must be a pythagorean triple.

Based on this we can go through possible (even) values of c while advancing a corresponding h sequentially. As c increases, so will h which will allow us to simply match h^2 with a^2 - (c/2)^2 to detect pythagorean triples.

def almostEqui(maxPerimeter):
    h = 1
    for c in range(2,maxPerimeter//3+2,2):
        cc4 = c*c//4
        for a in (c-1,c+1):
            hh = a*a - cc4
            while h*h < hh:
                h += 1
            if h*h == hh and 2*a+c <= maxPerimeter:
                yield a,a,c
        h -= 1 

def euler94():
    return sum(a+b+c for a,b,c in almostEqui(1000000000))

The euler94() function returns 518408346 in 2.21 minutes (133 seconds) on my laptop. This is roughly 2.5x faster than computing square roots.

[EDIT] An even faster method would be to generate pythagorean primitive triples and check that they form an almost equilateral triangle when doubled up:

This is a generator function based on Euclid's coprime formula for Primitive Pythagorean triples :

from math import gcd
def genTriples(k): # primitive pythagorean triangles a,b,c where a<b<c<k
    n,m = 1,2
    while m*m+1<k:                  # while c<k (for largest m producing c)
        if n>=m: n,m = m%2,m+1      # n reached m, advance m, reset n
        c = m*m+n*n                 # compute c 
        if c >= k: n=m;continue     # skip remaining n when c >= k
        if gcd(n,m) == 1:           # trigger on coprimes
            yield m*m-n*n,2*m*n,c   # return a,b,c triple
        n += 2                      # advance n, odd with evens

output:

from time import time
start=time()

N = 1000000000
result = 0
for a,b,c in genTriples(N//3+2):
    if abs(2*a-c)==1:               # (c,c,2*a) almost equilateral triangle
        result += 2*a+2*c
    if abs(2*b-c)==1:               # (c,c,2*b) almost equilateral triangle
        result += 2*b+2*c

print(time()-start)  # 51 seconds
print(result)        # 518408346

This is more than 2x faster than the previous approach.

Answer 2

You could decide to use functions:

def y(x):
  s = x**2 - 1
  if s % 3 == 0:
    h = (s / 3) ** 0.5
    if int(h) == h and h > 0:
      return int(h)

def hypoteneus(x):
  s = 2 * x
  uppr, lowr = s + 1, s - 1
  if uppr % 3 == 0:
    return uppr // 3
  if lowr % 3 == 0:
    return lowr // 3

def triangle(x):
  hyp = hypoteneus(x)
  height = y(x)
  if height is not None and hyp is not None:
    base = (hyp ** 2 - height **2) ** 0.5
    if base == int(base) and base > 0:
      return hyp, hyp, int(base) * 2

def  almost_eq_triangles(x):
  per = 0
  p1 = 0
  for i in range(1, int(1e9)):
    b = triangle(i)
    if b:
      p1 = int(sum(b))
      per +=p1
      if p1 > x:
        print('-'*(68 + len(str(p1))))
        print(f'{i:<10}: {str(b):<40}\t Perimeter: {p1}\tSum:{per}')
        print(f'\n\t\t{"Total sum of Perimeters (exclude last)"}: {per - p1}')
        print(f'\t\t{"Total sum of Perimeters (include last)"}: {per}')
        break
      print(f'{i:<10}: {str(b):<40}\t Perimeter: {p1}')

For example if I want all almost equilateral triangles whose perimeter is less than 10000000 we could do:

     almost_eq_triangles(10000000)
7         : (5, 5, 6)                                    Perimeter: 16
26        : (17, 17, 16)                                 Perimeter: 50
97        : (65, 65, 66)                                 Perimeter: 196
362       : (241, 241, 240)                              Perimeter: 722
1351      : (901, 901, 902)                              Perimeter: 2704
5042      : (3361, 3361, 3360)                           Perimeter: 10082
18817     : (12545, 12545, 12546)                        Perimeter: 37636
70226     : (46817, 46817, 46816)                        Perimeter: 140450
262087    : (174725, 174725, 174726)                     Perimeter: 524176
978122    : (652081, 652081, 652080)                     Perimeter: 1956242
3650401   : (2433601, 2433601, 2433602)                  Perimeter: 7300804
----------------------------------------------------------------------------
13623482  : (9082321, 9082321, 9082320)                  Perimeter: 27246962

                Total sum of Perimeters (exclude last): 9973078
                Total sum of Perimeters (include last): 37220040

Answer 3

The 'unnecessary' numbers are due to floating point precision and rounding error. For example, for a triangle of sides (93686, 93686, 93687), computing the parts of the formula separately we get area = 93687 * 81134.16730176011 / 2. This is clearly not an integer, yet area % 1 == 0 returns True.

You could try using Decimal to increase precision. Or instead, you could use the theorem that an integer is either a perfect square or its square root is an irrational number. For the area to be integral the term inside the square root must be a perfect square. b must also be an even number, otherwise a**2 - (b/2)**2 cannot be an integer.

Using a function is_square to check if the term is a square, your code can be corrected and slightly optimized like this:

from math import isqrt

def is_square(i):
    return i == isqrt(i) ** 2

for a in range(3, o, 2):  
    if is_square(a ** 2 - ((a+1)//2) ** 2):
        p = 3 * a + 1
        print(a, p)
        plst += p
    if is_square(a ** 2 - ((a-1)//2) ** 2):
        p = 3 * a - 1
        print(a, p)
        plst += p

(The condition type(p) == int is redundant because p will always be an integer; 2*a > b and b > 0 are also redundant if we start the loop at a = 2)

This is probably about as fast as you can get implementing a brute force algorithm and I leave it to others to demonstrate much more efficient algorithms which must exist using maths.