Construct Pythagorean triangles with given constraints

This post will present some ideas related to the generation of all Pythagorean triangles satisfying a certain criterion, from an algorithmic point of view. Of course, there are infinitely many such triangles (integer sided with a right angle). Denoting by and the edges adjacent to the right angle and with the hypothenuse, we have the classical relation

When are all integers, it is possible to give a more precise result regarding the possible values of . Indeed, for any such triangle, there exist integers with , not both odd and coprime such that

These formulas are attributed to Euclid. The case corresponds to coprime and such a triangle is caled “primitive”. What is important to note is that this formula generates ALL pythagorean triangles exactly once, and this allows us to solve the following two questions:

1. Generate all pythagorean triangles with edges .

2. Generate all pythagorean triangles with or .

3. Generate all pythagorean triangles with .

Now let’s answer these questions one at a time:

1. In order to generate pythagorean triangles with edges we need to loop over all , coprime, not both odd with such that and for each primitive triangle, add to the list all its “multiples” by looking at the appropriate . One possible algorithm is

• set an empty list
• loop for to
  • loop for to , , not both odd
  • loop for to and add the triangle to the list
• in the end will contain the desired triangles

2. In this question one of the two legs or should be equal to . The only difficulty here is that we need to loop on the divisors of . Indeed, we have two cases:

• : for all divisors of , we should find the possible factorizations with not both odd and coprime, and then add the corresponding triangle to the list.
• : find all factorizations and check again that obtained are coprime and not both odd.

3. In this case . Therefore, one should loop on all divisors of and in each case solve

where , are coprime, not both odd. This can be done again with a loop.

These tools might come in handy when working on Project Euler problems, since often when dealing with integer sided quantities in a triangle, things can be reduced to pythagorean triangles. When you reach this step, it is enough to loop on these triangles and perform the requested operations.

Integer sided triangle and IA, IB and IC integers

Let be a triangle whose side lengths are positive integers. Denote by the incenter of the triangle and suppose also that the segments have integer lengths. Prove that the inradius of the triangle is an integer.

Solution: Denote by the projections of the incenter on , respectively. Use the classical notation for the lengths of the sides of the triangle. Moreover, use the notations , , . Using Pythagora’s theorem in triangles determined by we obtain

Moreover, it can be proved that if then . Using the hypothesis, it follows that are integers. Using this and the Pythagora’s relations above we find that should be an integer. However, this is not enough to conclude that would also be an integer.

Looking at the triangle , denoting and , and applying the sine rule we get

Now note that and , which gives

This complicated relation allows us to deduce that is rational. This means that with and coprime integers. Moreover, we saw that is an integer, which means that or . In the end we find that should always be an integer.

Now, is it possible that is only a half integer, i.e. is odd? The Pythagora’s relations above imply that

If is an odd integer then is also odd and of the form . Moreover, is also an integer, which by the above relation should also be odd, which means that is an integer of the form . In the end we arrive at

which is a contradiction. Therefore must be an integer!

Sum of floors of multiples of the Golden Ratio

Propose an algorithm for computing

for , where is the golden ratio. The sum is the sum of the first terms in the so called lower Wythoff sequence.

Solution: Computing a floor function and a multiplication is not complicated, therefore proposing a algorithm for computing is trivial. However, such an algorithm is not viable for .

The path to a much more efficient algorithm goes through the notion of Beatty sequence. For a positive irrational number the associated Beatty sequence is . This notion becomes interesting when noting the following fact. If is defined by then the Beatty sequences and cover all the positive integers. The latter sequence is called the associated complementary Beatty sequence. Two proofs of this fact can be found in the Wikipedia link above.

It is obvious now that is just the partial sum of the Beatty sequence associated to . Now let’s see what is the complementary sequence. A brief computation shows that the associated is given by , which shows that the terms which are not of the form are not far from it. In order to see how this can give an efficient algorithm, let’s follow the instructions below:

• denote by , the largest term in .
• looking at all the integers up to , in view of the result regarding the Beatty sequences, they are either of the form or of the form (in the associated complementary sequence).
• denote by the largest integer such that , which is given by . Then, it is not difficult to see that is the difference between the sum of all positive integers up to and .
• In the end we obtain

  

• by convention, state that or .

The last item above shows that a recursive algorithm can be implemented in order to compute . Moreover, since the computation of is reduced to the computation of where , the algorithm will converge very fast, since the sequence of upper bounds for converges exponentially to zero. Therefore, the algorithm obtained will have complexity and will work for extremely large values of , provided the language you are using can handle the expected precision.

The algorithm for computing shown above can be used, for example, in computing sums related to elements in the Wythoff table, which can be expressed using Fibonacci numbers and terms in .

Other values of , like for example (with ) lead to other types of sums for which the same type of algorithm can be applied. It is likely that even cases where is not explicitly related to through an integer may be handled through a similar procedure (to be confirmed).

Sum of the Euler Totient function

Given a positive integer , the Euler totient function is defined as the number of positive integers less than which are co-prime with (i.e. they have no common factors with ). There are formulas for computing starting from the factorization of . One such formula is

where the product is made over all primes dividing .

If you have to compute for all numbers less than a threshold then another property could be useful: is arithmetic, that is, whenever . Therefore you could store all values computed until and for computing the value there are two possibilities: is a prime power and then or is composite and with . Then use the stored values to compute .

I now come to the main point of this post: computing the sum of all values of the totient function up to a certain :

One approach is to compute each and sum them. I will call this the brute-force approach. For all numerical purposes I will use Pari-GP in this post. On my computer it takes less than a second to compute and about seconds to compute . This is super linear in time, since the algorithm computes the factorization for each and then sums the values. Using the sieve approach could improve the timing a bit, but the algorithm is still super linear.

In some Project Euler problems it is not uncommon to have to compute something like or even larger. Therefore, there must be more efficient ways to compute out there, so let’s study some of the properties of . In another post I dealt with the acceleration of the computation of the sum of the divisor function.

We have which is the number of pairs with such that . It is not difficult to see that the total number of such pairs is . Moreover, the possible values of are . Now, if for we search instead for pairs satisfying then we have with and we get

There fore the number of pairs with gcd equal to is . Now we arrive at an interesting recursive formula:

At a first sight this looks more complicated, but there is a trick to keep in mind whenever you see a summation over of terms of the form : these quantities are constant on large intervals. Indeed,

Therefore we can change the index of summation from to . The range of for which the interval contains more than one integer is of order . Indeed, . Therefore for we should have at least one integer in the interval . The part where is larger than corresponds to smaller than . Therefore, we can split into two sums, each of order . and get that

where in the last sum we must make sure that in order to avoid duplicating terms in the sum.

Therefore we replaced a sum until to two sums with upper bound . The complexity is not , but something like since we have a recursive computation. Nevertheless, with this new formula and using memoization, to keep track of the values of already computed, we can compute very fast:

is computed instantly (vs second with brute force)

takes second (vs seconds with brute force)

takes seconds (vs over minutes with brute force)

takes seconds

takes about minutes

etc. Recall that these computations are done in Pari GP, which is not too fast. If you use C++ you can compute in seconds, in second and in seconds and in under a minute, if you manage to get past overflow errors.

A hint for Project Euler Pb 613

The text for Problem 613 can be found here. The hint is the following picture 🙂

Project Euler – Problem 264

Today I managed to solve problem 264 from Project Euler. This is my highest rating problem until now: 85%. You can click the link for the full text of the problem. The main idea is to find all triangles ABC with vertices having integer coordinates such that

• the circumcenter O of each of the triangles is the origin
• the orthocenter H (the intersection of the heights) is the point of coordinates (0,5)
• the perimeter is lower than a certain bound

I will not give detailed advice or codes. You can already find a program online for this problem (I won’t tell you where) and it can serve to verify the final code, before going for the final result. Anyway, following the hints below may help you get to a solution.

The initial idea has to do with a geometric relation linking the points A, B, C, O and H. Anyone who did some problems with vectors and triangles should have come across the needed relation at some time. If not, just search for important lines in triangles, especially the line passing through O and H (and another important point).

Once you find this vectorial relation, it is possible to translate it in coordinates. The fact that points A, B, C are on a circle centered in O shows that their coordinates satisfy an equation of the form , where is a positive integer, not necessarily a square… It is possible to enumerate all solutions to the following equation for fixed , simply by looping over and . This helps you find all lattice points on the circle of radius .

Once these lattice points are found one needs to check the orthocenter condition. The relations are pretty simple and in the end we have two conditions to check for the sum of the x and y coordinates. The testing procedure is a triple loop. We initially have a list of points on a circle, from the previous step. We loop over them such that we dont count triangles twice: i from 1 to m, j from i+1 to m, k from j+1 to m, etc. Once a suitable solution is found, we compute the perimeter using the classical distance formula between two points given in coordinates. Once the perimeter is computed we add it to the total.

Since the triple loop has cubic complexity, one could turn it in a double loop. Loop over pairs and construct the third point using the orthocenter condition. Then just check if the point is also on the circle. I didn’t manage to make this double loop without overcounting things, so I use it as a test: use double loops to check every family of points on a given circle. If you find something then use a triple loop to count it properly. It turns out that cases where the triple loop is needed are quite rare.

So now you have the ingredients to check if on a circle of given radius there are triangles with the desired properties. Now we just iterate over the square of the radius. The problem is to find the proper upper bound for this radius in order to get all the triangles with perimeter below the bound. It turns out that a simple observation can get you close to a near optimal bound. Since in the end the radii get really large and the size of the triangles gets really large, the segment OH becomes small, being of fixed length 5. When OH is very small, the triangle is almost equilateral. Just use the upper bound for the radius for an equilateral triangle of perimeter equal to the upper bound of 100000 given in the problem.

Using these ideas you can build a bruteforce algorithm. Plotting the values of the radii which give valid triangles will help you find that you only need to loop over a small part of the radii values. Factoring these values will help you reduce even more the search space. I managed to  solve the problem in about 5 hours in Pari GP. This means things could be improved. However, having an algorithm which can give the result in “reasonable” time is fine by me.

I hope this will help you get towards the result.

Project Euler 607

If you like solving Project Euler problems you should try Problem number 607. It’s not very hard, as it can be reduced to a small optimization problem. The idea is to find a path which minimizes time, knowing that certain regions correspond to different speeds. A precise statement of the result can be found on the official page. Here’s an image of the path which realizes the shortest time:

 

Project Euler tips

A few years back I started working on Project Euler problems mainly because it was fun from a mathematical point of view, but also to improve my programming skills. After solving about 120 problems I seem to have hit a wall, because the numbers involved in some of the problems were just too big for my simple brute-force algorithms.

Recently, I decided to try and see if I can do some more of these problems. I cannot say that I’ve acquired some new techniques between 2012-2016 concerning the mathematics involved in these problems. My research topics are usually quite different and my day to day programming routines are more about constructing new stuff which works fast enough than optimizing actual code. Nevertheless, I have some experience coding in Matlab, and I realized that nested loops are to be avoided. Vectorizing the code can speed up things 100 fold.

So the point of Project Euler tasks is making things go well for large numbers. Normally all problems are tested and should run within a minute on a regular machine. This brings us to choosing the right algorithms, the right simplifications and finding the complexity of the algorithms involved.

Project Euler Problem 285

Another quite nice problem from Project Euler is number 285. The result of the problem depends on the computation of a certain probability, which in turn is related to the computation of a certain area. Below is an illustration of the computation for k equal to 10.

To save you some time, here’s a picture of the case k=1 which I ignored and spent quite a lot of time debugging… Plus, it only affects the last three digits or so after the decimal point…

Here’s a Matlab code which can construct the pictures above and can compute the result for low cases. To solve the problem, you should compute explicitly all these areas.

function problem285(k)

N = 100000;

a = rand(1,N);
b = rand(1,N);

ind = find(abs(sqrt((k*a+1).^2+(k*b+1).^2)-k)<0.5);

plot(a(ind),b(ind),'.');
axis equal

M = k;
pl = 1;

for k=1:M
if mod(k,100)==0
k
end
r1 = (k+0.5)/k;
r2 = (k-0.5)/k;

f1 = @(x) (x<=(-1/k+r1)).*(x>=(-1/k-r1)).*(sqrt(r1^2-(x+1/k).^2)-1/k).*(x>=0).*(x<=1); f1 = @(x) f1(x).*(f1(x)>=0);
f2 = @(x) (x<=(-1/k+r2)).*(x>=(-1/k-r2)).*(sqrt(r2^2-(x+1/k).^2)-1/k).*(x>=0).*(x<=1); f2 = @(x) f2(x).*(f2(x)>=0);

if k == pl
thetas = linspace(0,pi/2,200);
hold on
plot(-1/k+r1*cos(thetas),-1/k+r1*sin(thetas),'r','LineWidth',2);
plot(-1/k+r2*cos(thetas),-1/k+r2*sin(thetas),'r','LineWidth',2);
plot([0 1 1 0 0],[0 0 1 1 0],'k','LineWidth',3);
hold off
axis off
end

A(k) = integral(@(x) f1(x)-f2(x),0,1);

end

xs = xlim;
ys = ylim;

w = 0.01;
axis([xs(1)-w xs(2)+w ys(1)-w ys(2)+w]);

sum((1:k).*A)

Fireworks – Project Euler 317

You can read the text of the problem here. The idea is that we have an explosion at a given height in a uniform gravitational field (no friction/wind). Supposing that all particles go outwards from the explosion point with a constant initial speed, what is the shape of the body formed by these particles? A hint is given in the picture below with a nice Matlab simulation.

Now what happens if we add a little bit of wind? Here’s the perturbation obtained when adding a 10m/s uniform wind speed vs the initial configuration. Something like this would be a lot more challenging 🙂

Project Euler Problem 144 – Laser light escaping an ellipse

Project Euler – problem 144 – visualization of the solution in Matlab

Linear programming #1 – Project Euler 185

I was recently faced with a very nice challenge from Project Euler – Problem 185. The idea is to find a number with a fixed number of digits by making guesses. Each time you are told how many digits you got right. By got right, I mean guessing the right digit on the right position. If you find a digit which is in the number, but on a different position, you don’t know…

There is a test case with 5 digits:

90342 ;2 correct
70794 ;0 correct
39458 ;2 correct
34109 ;1 correct
51545 ;2 correct
12531 ;1 correct

with unique solution 39542 and the difficult case with 16 digits

5616185650518293 ;2 correct
3847439647293047 ;1 correct
5855462940810587 ;3 correct
9742855507068353 ;3 correct
4296849643607543 ;3 correct
3174248439465858 ;1 correct
4513559094146117 ;2 correct
7890971548908067 ;3 correct
8157356344118483 ;1 correct
2615250744386899 ;2 correct
8690095851526254 ;3 correct
6375711915077050 ;1 correct
6913859173121360 ;1 correct
6442889055042768 ;2 correct
2321386104303845 ;0 correct
2326509471271448 ;2 correct
5251583379644322 ;2 correct
1748270476758276 ;3 correct
4895722652190306 ;1 correct
3041631117224635 ;3 correct
1841236454324589 ;3 correct
2659862637316867 ;2 correct

While the small case could be tackled using a pure brute force approach (only ) cases to check, the second case becomes intractable. Looping through all cases would take ages. One could try genetic algorithms or other algorithms based on randomness. Recursive approaches are also possible.

There is however a way to express this problem which makes it solvable right away using a linear programming solver: that is an algorithm which finds valid solutions given some constraints. How do we find these constraints? There are essentially two types of constraints in our case:

• each of the positions of the result contains one digit from 0 to 9
• in each of the guesses given above we have a fixed number of correct digits