Fermat points and parameterizing the 120 degree integer triangle (Project Euler 143)

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 $ABC$ , with all angles smaller than $120^\circ$ . Point $T$ is inside the triangle such that $p+q+r$ is at its minimum ( $T$ can always be determined uniquely).

We define a Torricelli triangle as one where $a$ , $b$ , $c$ , $p$ , $q$ , $r$ are all integers.

Generate all Torricelli triangles with $p + q + r \leq 120000$ .

The Fermat point

In triangle $ABC$ , with $P$ being the point such that $PA+PB+PC$ is minimized, $P$ is called the Fermat point of triangle $ABC$ .

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

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

We now show that $P$ must be at the intersection of $AD$ and the circumcircle of $\triangle BCD$ .

Assume that $P$ does not lie on the circumcircle of $\triangle BCD$ . Then, according to Ptolemy’s theorem:

$BP \times CD + PC \times BD > PD \times BC$

Equality only occurs when $PCDB$ is cyclic. Since $BCD$ is equilateral, we have $BC = CD = BD$ so we can divide it out to get

$PB + PC > PD$

Adding $PA$ to both sides:

$PA + PB + PC > PA + PD$

Now if $P$ did lie on the circumcircle of $\triangle BCD$ , we would have $PA + PB + PC = PA + PD$ , so $P$ must lie on the circumcircle of $\triangle BCD$ in order to be optimal.

So we know $P$ is on the circle. But now we assume that $P$ is not on $AD$ . As we had before, $PA + PB + PC = PA + PD$ .

Next if $P$ is not on $AD$ , then $AP + PD > AD$ , or by substitution,

$PA + PB + PC > AD$

Of course if $P$ were actually on $AD$ , then $PA + PB + PC = AD$ , and the sum $PA+PB+PC$ would be optimal.

This proves the optimal place for $P$ to be on the intersection of $AD$ and the circumcircle of $\triangle BCD$ .

If we repeat this for the other three sides, we notice that $P$ is the intersection of the three circumcircles, and also the intersection of $AD$ , $BE$ , and $CF$ :

Further, we see that quadrilaterals $BPCD$ , $APCE$ , and $FAPB$ are all cyclic. As $\angle BDC = 60^\circ$ , $\angle BPC = 120^\circ$ , and similarly $\angle APC = 120^\circ$ and $\angle APB = 120^\circ$ .

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 $120^\circ$ , we can apply the cosine law ( $\cos 120^\circ = -\frac{1}{2}$ ):

$a^2 = q^2 + r^2 - 2qr \cos 120$

$a^2 = q^2 + r^2 +qr$

A similar formula applies to sides $b$ and $c$ . We call $(x,y)$ a pair if $x^2 + xy + y^2$ is a perfect square.

We have found a Torricelli triangle if for three integers $p$ , $q$ , $r$ , all of $(p,q)$ , $(p,r)$ , $(q,r)$ are all pairs.

This leaves us with an outline of an algorithm:

Generate all pairs $(x,y)$ under the limit (with $x+y < 120000$ and $x^2 + xy + y^2$ being a square)
Sort the list of pairs and index them (to be explained in step 3)
For each pair $(a,b)$ in the list, search through the list to check if there exists some $c$ where $(a,c)$ is a pair and $(b,c)$ is a pair. We index the list to drastically improve searching time.
If an $(a,b,c)$ triple is found and $a+b+c \leq 120000$ , then mark $a+b+c$ 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 $120000^2$ 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 $a^2 + b^2 = c^2$ can be generated by $a = m^2 - n^2$ , $b = 2mn$ , $c = m^2 + n^2$ where $m$ and $n$ are coprime and are of opposite parity.

Can a similar parameterization be found for the equation $a^2 + ab + b^2$ equation?

Letting $x = \frac{a}{c}$ , and $y = \frac{b}{c}$ and dividing by $c^2$ , we get

$x^2 + xy + y^2 = 1$

which is the equation of an ellipse.

Originally we wanted to find integer solutions to the equation $a^2 + ab + b^2 = c^2$ . This is equivalent to finding all rational points on the ellipse $x^2 + xy + y^2 = 1$ :

It is easy to see why: if $x$ and $y$ are rational points such that $x^2 + xy + y^2 = 1$ then $x = \frac{a}{c}$ and $y = \frac{b}{c}$ where $a$ , $b$ , $c$ are positive integers, and we arrive back at the original form of the equation.

Also in order for $a$ , $b$ , $c$ to be positive, we will only consider points on the ellipse in the first quadrant, with both $x$ and $y$ 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 $(0,-1)$ to be the first point of the line. Then the equation of the line is $y = tx-1$ , where $t$ is valid only if it is positive and greater than 1.

Substituting into the equation of the ellipse:

$x^2 +tx^2 - x + t^2 x^2 - 2tx + 1 = 1$

Simplifying:

$x (t^2+t+1) = 2t + 1$

$x = \frac{2t+1}{t^2 + t + 1}$

Now evaluating for $y$ :

$y = t (\frac{2t+1}{t^2 + t + 1})-1 = \frac{t^2-1}{t^2+t+1}$

Notice now that for $x$ and $y$ to be rational, we just need the slope $t$ to be rational. So we can write $t$ as $\frac{m}{n}$ where $m$ and $n$ are positive integers and $m>n$ .

Simplifying the expressions for $x$ and $y$ :

$x = \frac{2 \frac{m}{n} + 1}{(\frac{m}{n})^2 + \frac{m}{n} + 1} = \frac{2mn+n^2}{m^2 + mn + n^2}$

$y = \frac{(\frac{m}{n})^2 - 1}{(\frac{m}{n})^2 + \frac{m}{n} + 1} = \frac{m^2 - n^2}{m^2 + mn + n^2}$

Since we defined $x = \frac{a}{c}$ and $y = \frac{b}{c}$ , we have a valid triple if $a = 2mn+n^2$ , $b = m^2 - n^2$ , and $c = m^2 + mn + n^2$ .

If $m$ and $n$ are not coprime, then neither will $a$ , $b$ , and $c$ be coprime as they will share a common divisor. So in order for $(a,b,c)$ to be primitive we will need $m$ and $n$ to be coprime.

Still this isn’t quite sufficient. Consider what would happen if $m \equiv n (\mod 3)$ . We have $a = 2mn+n^2$ which can be written as $n(n-m+3m)$ which is divisible by 3. Next $b = (m+n) (m-n)$ which is divisible by 3. Finally, $c = m^2 + mn + n^2$ which can be written as $(m-n)^2 + 3mn$ , again divisible by 3! So if $m \equiv n \mod 3$ , then the resulting triple is not primitive!

But it turns out that we can ignore all cases where $m \equiv n \mod 3$ . Given a triple $(2mn+n^2, m^2-n^2, m^2+mn+n^2)$ , with $m \equiv n \mod 3$ , it is possible to find an identical parameterization giving $(\frac{2mn+n^2}{3}, \frac{m^2-n^2}{3}, \frac{m^2 + mn + n^2}{3})$ .

As $m \equiv n mod 3$ , we have $m+2n \equiv 0 \mod 3$ and also $m-n \equiv 0 \mod 3$ . Thus we can have $u$ and $v$ where $u$ and $v$ are positive integers:

$u = \frac{m+2n}{3}$

$v = \frac{m-n}{3}$

Multiplying by 3, $3u = m+2n$ , and $3v = m-n$ . Combining the two and substituting we get

$m = u+2v$

$n = u-v$

If we substitute $u$ and $v$ for $m$ and $n$ in the triple $(2mn+n^2, m^2-n^2, m^2+mn+n^2)$ , we get:

$2mn + n^2 = 3 (u^2-v^2)$

$m^2 - n^2 = 3(2uv + v^2)$

$m^2 + mn + n^2 = 3 (u^2 + uv + v^2)$

So this proves that when $m \equiv n \mod 3$ , 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 $n$ operations instead of $n^2$ operations. So we can loop up to $\sqrt{120000}$ , 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;
}

Lucky's Notes

Notes on math, coding, and other stuff