Problem 143 of Project Euler is an interesting geometry problem, and one notorious for having one of the lowest solvers to time released ratios of all the problems on the site: under 600 solvers in over 3 years as time of writing.

So what is the problem?

The problem is quite simple. We have a triangle , with all angles smaller than . Point is inside the triangle such that is at its minimum ( can always be determined uniquely).

We define a *Torricelli triangle* as one where , , , , , are all integers.

Generate all Torricelli triangles with .

### The Fermat point

In triangle , with being the point such that is minimized, is called the **Fermat point** of triangle .

We can use Ptolemy’s theorem to locate the Fermat point in the triangle.

Let be a point such that is equilateral, then construct the circumcircle of :

We now show that must be at the intersection of and the circumcircle of .

Assume that does *not* lie on the circumcircle of . Then, according to Ptolemy’s theorem:

Equality only occurs when is cyclic. Since is equilateral, we have so we can divide it out to get

Adding to both sides:

Now if *did* lie on the circumcircle of , we would have , so must lie on the circumcircle of in order to be optimal.

So we know is on the circle. But now we assume that is not on . As we had before, .

Next if is not on , then , or by substitution,

Of course if were actually on , then , and the sum would be optimal.

This proves the optimal place for to be on the intersection of and the circumcircle of .

If we repeat this for the other three sides, we notice that is the intersection of the three circumcircles, and also the intersection of , , and :

Further, we see that quadrilaterals , , and are all cyclic. As , , and similarly and .

### Designing an algorithm

We have enough information now to start working on an algorithm. Let us come back to the previous diagram:

So knowing that the central angles are all , we can apply the cosine law ():

A similar formula applies to sides and . We call a *pair* if is a perfect square.

We have found a Torricelli triangle if for three integers , , , all of , , are all pairs.

This leaves us with an outline of an algorithm:

- Generate all pairs under the limit (with and being a square)
- Sort the list of pairs and index them (to be explained in step 3)
- For each pair in the list, search through the list to check if there exists some where is a pair and is a pair. We index the list to drastically improve searching time.
- If an triple is found and , then mark as found. This is easily implemented as a bit array of size 120000 with all bits initially 0 and flipped to 1 when a sum is encountered.
- Add up the sums and output

Step 1 takes up most of the time in this algorithm. To do this by ‘brute force’, we need an outer loop from 1 up to 120000, then another inner loop again from 1 to 120000 (actually a bit less), which is essentially operations, which is in the billions. An implementation of this algorithm would take about 3 minutes.

### Parameterizing the a^2 + ab + b^2 = c^2 equation

It is well known that all primitive triples of the form can be generated by , , where and are coprime and are of opposite parity.

Can a similar parameterization be found for the equation equation?

Letting , and and dividing by , we get

which is the equation of an ellipse.

Originally we wanted to find *integer *solutions to the equation . This is equivalent to finding all *rational *points on the ellipse :

It is easy to see why: if and are rational points such that then and where , , are positive integers, and we arrive back at the original form of the equation.

Also in order for , , to be positive, we will only consider points on the ellipse in the first quadrant, with both and being positive. We do this by choosing a point on the ellipse, and from there combine the equations of the ellipse and the line.

Let us choose the point to be the first point of the line. Then the equation of the line is , where is valid only if it is positive and greater than 1.

Substituting into the equation of the ellipse:

Simplifying:

Now evaluating for :

Notice now that for and to be rational, we just need the slope to be rational. So we can write as where and are positive integers and .

Simplifying the expressions for and :

Since we defined and , we have a valid triple if , , and .

If and are not coprime, then neither will , , and be coprime as they will share a common divisor. So in order for to be primitive we will need and to be coprime.

Still this isn’t quite sufficient. Consider what would happen if . We have which can be written as which is divisible by 3. Next which is divisible by 3. Finally, which can be written as , again divisible by 3! So if , then the resulting triple is not primitive!

But it turns out that we can ignore all cases where . Given a triple , with , it is possible to find an identical parameterization giving .

As , we have and also . Thus we can have and where and are positive integers:

Multiplying by 3, , and . Combining the two and substituting we get

If we substitute and for and in the triple , we get:

So this proves that when , the case is nonprimitive and has already been covered by a smaller (possibly primitive) case. We are done.

### Implementation in C++

Now using the new parameterization scheme, we can generate all pairs in about operations instead of operations. So we can loop up to , which cuts down the time down to about 2 seconds:

#include <cstdio>
#include <vector>
#include <algorithm>
#include <utility>
using namespace std;
// Limit and square root of limit
#define L 120000
#define LSQ 347
typedef long long int int64;
int64 gcd(int64 a, int64 b)
{
int64 c = a % b;
while(c != 0)
{
a = b;
b = c;
c = a % b;
}
return b;
}
int main()
{
// Store pairs in here
vector< pair<int64,int64> > pairs;
// Use the parameterization
for(int64 u=1; u<LSQ; u++){
for(int64 v=1; v<u; v++){
if(gcd(u,v) != 1) continue;
if((u-v) % 3 == 0) continue;
int64 a = 2*u*v + v*v;
int64 b = u*u - v*v;
if(a+b > L) break;
// From coprime pairs make composite pairs
for(int k=1; k*(a+b)<L; k++){
pairs.push_back(make_pair(k*a,k*b));
pairs.push_back(make_pair(k*b,k*a));
}
}
}
// Sort pairs list
sort(pairs.begin(),pairs.end());
// Create index
int index[L];
for(int i=0; i<L; i++) index[i] = -1;
for(int i=0; i<pairs.size(); i++){
if(index[pairs[i].first] == -1)
index[pairs[i].first] = i;
}
// Which sums have been reached?
bool sums[L];
for(int i=0; i<L; i++) sums[i] = false;
// Iterate through all pairs
for(int i=0; i<pairs.size(); i++){
int64 a = pairs[i].first;
int64 b = pairs[i].second;
// Construct vectors for indices
vector<int64> va;
vector<int64> vb;
// Fetch indices
int ia = index[a];
int ib = index[b];
while(ia<pairs.size()){
pair<int64,int64> next = pairs[ia];
if(next.first != a) break;
va.push_back(next.second);
ia++;
}
while(ib<pairs.size()){
pair<int64,int64> next = pairs[ib];
if(next.first != b) break;
vb.push_back(next.second);
ib++;
}
// Find common objects between va and vb
for(int ja=0; ja<va.size(); ja++){
for(int jb=0; jb<vb.size(); jb++){
if(va[ja]==vb[jb]){
// Potential c found
int64 c = va[ja];
if(a+b+c<L) sums[a+b+c] = true;
}
}
}
}
// Tally up sums
int64 s = 0;
for(int i=0; i<L; i++)
if(sums[i]) s+=i;
printf("%lld\n",s);
return 0;
}