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:

  1. Generate all pairs (x,y) under the limit (with x+y < 120000 and x^2 + xy + y^2 being a square)
  2. Sort the list of pairs and index them (to be explained in step 3)
  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.
  4. 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.
  5. 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;
}

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

  1. Anonymous says:

    implemented in python (w/o pythony optimization) ~ 18.5 s

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 63 other followers

%d bloggers like this: