Vrftsjtryk

I am trying to solve Project Euler Problem 224:

Let us call an integer sided triangle with sides a ≤ b ≤ c barely obtuse if the sides satisfy a² + b² = c² - 1.

How many barely obtuse triangles are there with perimeter ≤ 75,000,000?

For large v_max and v_max_n, the generator below yields very slow when I have 750000000 and 500000000, respectively. Having tried for days, I am still looking for an algorithm that can do this in under a minute.

from itertools import product

import math



n = 15 * 10 ** 8

v_max = math.floor(n / 2 - 1 / (2 * n))

v_max_n = n // 3



(i for i in product(range(2, v_max_n + 1, 2), 

                    range(2, v_max + 1, 2), 

                    range(3, v_max + 1)) 

 if i[0] <= i[1] and i[1] <= i[2] and sum(i) <= n and

 i[0] ** 2 + i[1] ** 2 == i[2] ** 2 - 1)

I have tried using it as map(lambda x: x, gen), nested for loops, and have generally been racking my brain to no avail.

I tried re-writing b, c as expressions. This tends to be quick but misses some of the triples.

n = 12

v_max = math.floor(n / 2 - 1 / (2 * n))

v_max_n = n // 3

j = 0  # tested this with low n and compared: 12, 150 to see 

for a in range(1, v_max_n + 1, 1):

    b = (n ** 2 - 2 * n * a - 1) / (2 * (n - a))

    c = (a ** 2 + b ** 2 + 1) ** 0.5

    if c % 1 == 0 and b % 1 == 0 and a + int(b) + int(c) <= n:

        # print(a, int(b), int(c))

        j += 1

print(j)

It seems like you want to find a solution for a²+b²=c²-1 with some added constraints, like a<=b<=c, a+b+c<=n, etc. by trying all possible combinations of a, b, and c and testing whether they satisfy all your conditions at once. This is pretty wasteful.

Instead of itertools.product, you can use a generator expression with three regular for loops. This has several advantages:

• 
  

  you could move some of the if clauses up after the first or second loop, i.e. you don't have to run the other loops if those conditions already fail at this point, e.g.

  
  
  

   for x in ... if test(x) for y in ... if test(y) ...
  
  

  
  



• you can use the values of your previous variables in the lower and upper bounds for the later variables (in your case this replaces most of your if conditions)



• you can basically calculate the value for c from a and b without testing all possible values

This is how I'd do it:

gen = ((a, b, int(c))

       for a in range(2, v_max_n + 1, 2)

       for b in range(a, min(v_max + 1, n - a - a), 2) # use a for lower and upper bounds

       for c in [(a ** 2 + b ** 2 + 1)**0.5]           # just a single candidate for c

       if  c % 1 == 0)                                 # whole-numbered c found?

Note that the calculation of c using very_large_number**0.5 might be imprecise with float; using decimal might work, though. However, even with those optimizations, testing much fewer values than your original loop (on the order of O(n²) instead of O(n³)), it might not be feasible to find solution for large values of a, b and c.

Also note that I did not thoroughly test this for off-by-one errors, since performance was the main concern.

Another thing that you might try: Invert the loops for a and for b, i.e. instead of testing all b for each a, no matter how large the difference, test all a up to b first. This does not decrease the number of combinations to test, but it seems to yield much more valid combinations in the "smaller" number much faster. Like this (upper bound for b can probably be reduced):

       for b in range(2, v_max + 1, 2)

       for a in range(2, b+1, 2)

v_max_n = n // 3

The name of this variable is a mystery to me, but you seem to be using it as an upper bound on $a$. You can get a notably better upper bound on $a$ very easily: since $a le b$ and $c^2 = a^2 + b^2 + 1$, $a + b + c ge 2a + sqrt{2a^2 + 1} > 2a + sqrt{2a^2} = (2 + sqrt 2)a$. Then your upper bound on $a$ can be reduced by 12.13%.

Given a value of $a$, the constraints on $b$ are $a le b$ and $a + b + sqrt{a^2+b^2+1} le n$. Rearrange and square: $a^2 + b^2 + 1 le n^2 + a^2 + b^2 - 2an +2ab -2bn$ or $b le frac{n^2 - 2an - 1}{2n - 2a}$. I'm not sure why you appear to be using this bound as a forced value in the code you added to the question.

$$sum_{a=1}^{frac{n}{2 + sqrt 2}} left(frac{n^2 - 2an - 1}{2n - 2a} - aright)
= left(sum_{a=1}^{frac{n}{2 + sqrt 2}} frac{n^2 - 2an - 1}{2n - 2a}right) - frac{n(n + 2 + sqrt 2)}{2(2 + sqrt 2)^2} \
= frac{n^2 + 1}2 left(psileft(frac{n}{2 + sqrt 2}-n+1right) - psi(1-n)right) + frac{n^2}{2 + sqrt 2} - frac{n(n + 2 + sqrt 2)}{2(2 + sqrt 2)^2} \
= frac{n^2 + 1}2 left(psileft(1-frac{1+sqrt2}{2 + sqrt 2}nright) - psi(1-n)right) + frac{(3 + 2sqrt 2)n^2 + (2 + sqrt 2)n}{2(2 + sqrt 2)^2} \
> 0.9n^2
$$
The step from the sum to an expression in terms of the digamma function is courtesy of Wolfram Alpha.

As a rule of thumb, you don't want to be repeating a loop body $10^{12}$ times, so looping over all values of $a$ and all values of $b$ is a non-starter.

Hint: try rearranging $a^2 + b^2 = c^2 - 1$.

One improvement over @tobias_k s answer is to add the extra upper bound on b of b<1 + a*a//2. This is valid because after that, the distance between b^2 and (b+1)^2 will be bigger than a^2, so c^2 does not exist. This provides a massive speedup at least a first.