# Calculate the Fermat Point of a Triangle

This is somewhat similar to The centers of a triangle, but with a different point. The Fermat Point is the point P in triangle ABC such that the value of AP + BP + CP is minimized. There are two cases:

If there is an angle greater than 120 degrees, that vertex is the fermat point. Otherwise, draw equilateral triangles on each of the sides of ABC. Connect the far vertex of each equilateral triangle to the opposite vertex of triangle ABC. Doing this for each of the three equilateral triangles results in a single common point of intersection for all three lines, which is the Fermat Point.

It should run within 5 seconds on a reasonable machine.

Input: A set of 3 points, not necessarily integers. This can be taken as a nested array, string, list of tuples, etc. (whatever suits your language).

Output: The coordinates of the Fermat point, again, however your language best handles points. Floating point inaccuracies will not be counted against you.

Test Cases:

[[1, 1], [2, 2], [1, 2]] --> [1.2113248654051871, 1.788675134594813]
[[-1, -1], [-2, -1], [0, 0]] --> [-1, -1]
[[-1, -1], [1, -1], [0, 1]] --> [0, -0.42264973081037427]
[[0, 0], [0.5, 0.8660254037844386], [-5, 0]] --> [0, 0]
[[0, 0], [0, -5], [-0.8660254037844386, 0.5]] --> [0, 0]


This is code golf so shortest code wins!

• Is it OK to try all points in increments of floating point precision and select the one that minimizes the total distance?
– xnor
May 10 '16 at 3:26
• @xnor If you can do it within 5 seconds. May 10 '16 at 21:34
• Up to how many significant figures must the output be accurate to? Also, is it okay if -0.0 is output in place of some 0.0s? May 15 '16 at 2:56
• @R. Kap I'd say about 5 or 6 significant figures. There isn't so much much that rounding errors should be a problem. As for the second question, that seems fine. May 16 '16 at 20:22

infixl 5£
z=zipWith
(?)=z(-)
t[a,b]=[-b,a]
a¤b=sum$z(*)a b a%b=t a¤b r a b c=[c%b/a%b,c%a/a%b] x£y=2*x¤y<= -sqrt(x¤x*y¤y) f[a,b,c]|a?b£c?b=b|a?c£b?c=c|b?a£c?a=a|[n,m,p,o]<-c?k a b c++a?k b c a=r[m,o][n,p][c%[n,m],a%[p,o]] k a b c=map(/2)$z(+)a b?map(signum((b?a)%(c?a))*sqrt 3*)(t$b?a)  The same code with some explanations infixl 5£ z=zipWith -- operator ? : difference of two vectors (?)=z(-) -- function t : rotate a vector by +90 degrees t[a,b]=[-b,a] -- operator ¤ : scalar product of two vectors ( a¤b = a0 * b0 + a1 * b1 ) a¤b=sum$z(*)a b

-- operator % : "cross product" of two vectors ( a%b = a0 * b1 - a1 * b0 )
--      this returns actually the z coordinate of the 3d cross vector
--      other coordinates are nul since a and b are in the xy plan
a%b=t a¤b

-- function r : solves the system of two linear equations with two variables x0,x1 :
--      a0*x0 - b0*x1 = c0
--      a1*x0 - b1*x1 = c1
r a b c=[c%b/a%b,c%a/a%b]

-- operator £ : returns true if the angle between two vectors is >= 120 degrees
--      x¤y = ||x|| * ||y|| * cos(xyAngle)
--      so xyAngle>=120° is equivalent to : x¤y / (||x|| * ||y||) <= -0.5
x£y=2*x¤y<= -sqrt(x¤x*y¤y)

-- function k : takes 3 points A B C of a triangle and constructs the point C'
--              of the equilateral triangle ABC' which is opposite to C:
--              C' = (A+B)/2 - ((+/-) sqrt(3)/2 * t(AB))
--
--      the sign +/- is given by the sign of the cross vector of AB an AC ((b?a)%(c?a))
--      which is >0 if the angle between AB and AC is positive
--      and <0 otherwise.
k a b c=map(/2)$z(+)a b?map(signum((b?a)%(c?a))*sqrt 3*)(t$b?a)

-- function f : returns the fermat point of a triangle
f[a,b,c]
|a?b£c?b=b  -- return B if angle ABC >= 120°
|a?c£b?c=c  -- return C if angle BCA >= 120°
|b?a£c?a=a  -- return A if angle CAB >= 120°
|[n,m,p,o]<-c?k a b c++a?k b c a= -- calculate the two segments C'C and A'A
r[m,o][n,p][c%[n,m],a%[p,o]]  -- return their intersection


Tests:

main = do
print $f [[1, 1], [2, 2], [1, 2]] print$ f [[-1, -1], [-2, -1], [0, 0]]
print $f [[-1, -1], [1, -1], [0, 1]] print$ f [[0, 0], [0.5, 0.8660254037844386], [-5, 0]]
print \$ f [[0, 0], [0, -5], [-0.8660254037844386, 0.5]]


Output:

[1.2113248654051871,1.7886751345948126]
[-1.0,-1.0]
[0.0,-0.42264973081037427]
[0.0,0.0]
[0.0,0.0]

• How are you counting bytes? £ and ¤ are 2 bytes each in UTF-8, and I don’t know of a Haskell compiler that accepts ISO-8859-1. (You have plenty of free 1-byte ASCII operators, though, so this is easy to fix.) May 17 '16 at 22:15
• I'm counting it with my editor which actually counts characters. I didn't know that these were 2 bytes, but anyway, as you said, I could replace it by other 1 byte operators. This code compiles with GHC 7.8.3 May 19 '16 at 6:04
• GHC compiles your code when encoded as UTF-8 with £ and ¤ as 2 byte operators, but not when encoded as ISO-8859-1 with £ and ¤ as 1 byte operators. The available 1 byte operators in UTF-8 are !, #, %, &, ?. You should replace the 2 byte operators or adjust your byte count. May 19 '16 at 19:18

# Mathematica, 39 bytes

Sum[Norm[p-{x,y}],{p,#}]~NArgMin~{x,y}&


Constructs an equation based on the distances between the vertices and a point {x,y}. Then uses the NArgMin function to find a global minimum for that equation, which will be the Fermat Point by definition.

• 39 bytes, when the next shortest answer is 285... May 15 '16 at 7:34

# Python, 475448 440 bytes

Any help to golf further is appreciated.

from math import *
d=lambda x,y:((x[0]-y[0])**2+(x[1]-y[1])**2)**0.5
s=lambda A,B,C:(d(B,C), d(C,A), d(A,B))
j=lambda a,b,c:acos((b*b+c*c-a*a)/(2*b*c))
t=lambda a,b,c:1/cos(j(a,b,c)-pi/6)
b=lambda A,B,C,p,q,r:[(p*A[i]+q*B[i]+r*C[i])/(p+q+r) for i in [0,1]]
f=lambda A,B,C:A if j(*s(A,B,C)) >= 2*pi/3 else B if j(*s(B,C,A)) >= 2*pi/3 else C if j(*s(C,A,B)) >= 2*pi/3 else b(A,B,C,d(B,C)*t(*s(A,B,C)),d(C,A)*t(*s(B,C,A)),d(A,B)*t(*s(C,A,B)))


Ungolfed:

from math import *
#distance between two points
d = lambda x,y: ((x[0]-y[0])**2+(x[1]-y[1])**2)**0.5

#given the points, returns the sides
s = lambda A,B,C : (d(B,C), d(C,A), d(A,B))

#given the sides, returns the angle
j = lambda a,b,c : acos((b*b+c*c-a*a)/(2*b*c))

#given the sides, returns secant of that angle
t = lambda a,b,c: 1/cos(j(a,b,c)-pi/6)

#given the sides and the Trilinear co-ordinates, returns the Cartesian co-ordinates
b = lambda A,B,C,p,q,r: [(p*A[i]+q*B[i]+r*C[i])/(p+q+r) for i in [0,1]]

#this one checks if any of the angle is >= 2π/3 returns that point else computes the point
f = lambda A,B,C: A if j(*s(A,B,C)) >= 2*pi/3 else B if j(*s(B,C,A)) >= 2*pi/3 else C if j(*s(C,A,B)) >= 2*pi/3 else b(A,B,C,d(B,C)*t(*s(A,B,C)),d(C,A)*t(*s(B,C,A)),d(A,B)*t(*s(C,A,B)))


Input:

print('{}'.format(f([1, 1], [2, 2], [1, 2])))
print('{}'.format(f([-1, -1], [-2, -1], [0, 0])))
print('{}'.format(f([-1, -1], [1, -1], [0, 1])))
print('{}'.format(f([0, 0], [0.5, 0.8660254037844386], [-5, 0])))
print('{}'.format(f([0, 0], [0, -5], [-0.8660254037844386, 0.5])))


Output:

[1.2113248652983113, 1.7886751347016887]
[-1, -1]
[0.0, -0.42264973086764884]
[0, 0]
[0, 0]

• from math import* is a pretty common golf. This will also let you use pi instead of hard coding it (same length for 2*pi/3). You can also drop a lot of whitespace like: d=lambda x,y:(.... May 12 '16 at 13:26

# Python 3.5, 10191016998982969 953 bytes:

from math import*
def H(z,a,b):c=complex;T=lambda A,B:abs(c(*A)-c(*B));d=T(z,a);e=T(z,b);f=T(a,b);g=[d,e,f];h=max(g);g.remove(h);i=acos((sum(i*i for i in g)-(h*h))/(2*g[0]*g[-1]));_=[[z,a],[z,b],[a,b]];j,s,t=cos,sin,atan;N=[[b,a]for a,b in zip([b,a,z],[max(i,key=i.get)if i!=''else''for i in[{(g[0]+(h*j(t(l))),g[1]+(h*s(t(l)))):T(k,(g[0]+(h*j(t(l))),g[1]+(h*s(t(l))))),(g[0]-(h*j(t(l))),g[1]-(h*s(t(l)))):T(k,(g[0]-(h*j(t(l))),g[1]-(h*s(t(l)))))}if l else{(g[0]+h,g[1]):T(k,(g[0]+h,g[1])),(g[0]-h,g[1]):T(k,(g[0]-h,g[1]))}if l==0else''for g,h,l,k in zip([((a[0]+b[0])/2,(a[1]+b[1])/2)for a,b in _],[(3**0.5)*(i/2)for i in[d,e,f]],[-1/p if p else''if p==0else 0for p in[((a[1]-b[1])/(a[0]-b[0]))if a[0]-b[0]else''for a,b in _]],[b,a,z])]])if b!=''];I=N[0][0][1];J=N[0][0][0];K=N[1][0][1];G=N[1][0][0];A=(N[0][1][1]-I)/(N[0][1][0]-J);B=I-(A*J);C=(K-N[1][1][1])/(G-N[1][1][0]);D=K-(C*G);X=(D-B)/(A-C);Y=(A*X)+B;return[[X,Y],[[a,b][h==d],z][h==f]][i>2.0943]


Incredibly long compared to other answers, but hey, at least it works! I could not be happier with the result I got as this has got to be one of the hardest challenges I have ever done. I am just so happy that it actually works! :D Now, onto the more technical notes:

• This function takes each ordered pair in as a list or a tuple. For instance, H((1,1),(2,2),(1,2)) will work, but so will H([1,1],[2,2],[1,2]).
• Outputs the coordinates of the points in either a list of integers or floating points depending on whether or not one angle more than or equal to 120º exists.
• This may output -0.0 in place of 0.0 for some inputs. For instance, the output for the input [-1, -1], [1, -1], [0, 1] is [-0.0, -0.4226497308103744]. I hope this is okay, although if it isn't, I will change it, though it will cost me a few more bytes. This is okay, as confirmed by OP.
• Should be accurate up to at least 13 to 14 significant figures.

I will try and golf this more over time. An explanation, possibly very long, coming soon.

Try It Online! (Ideone)