Prolog: Find all equilateral triangles which are built from DB points

You have a couple of fundamental issues here.

From memory, the two-dimensional coordinates of an equilateral triangle cannot all be integers. You must have at least one irrational number. Computers are not good at representing irrational numbers.

So, you're also confusing real maths with computer maths. You can't simply get the sqrt and check for equality.

To illustrate I produced these three points:

point(j1,0,0).
point(j2,1,0).
point(j3,0.5,0.866025403784439). % 1/2 & sqrt(3)/2

That's an equilateral triangle.

When I run that with your code, the lengths that get produced are like this:

[1.0,1.0000000000000004,1.0000000000000004]

They are all effectively 1.0, but of course 1.0000000000000004 is not 1.0. So, even this doesn't comeback as a equal.

So you are really forced to check the confidence as an epsilon to say two numbers are equal.

Here's what I did to do that:

points_main(Xs):-
    findall(Xs1,equilateral_triangles(Xs1),Xs).

equilateral_triangles([P1,P2,P3]):-
    point(P1,X1,Y1),
    point(P2,X2,Y2),
    P1 \= P2,
    point(P3,X3,Y3),
    P1 \= P3,
    P2 \= P3,
    L12 is sqrt((X2 - X1)^2 + (Y2 - Y1)^2),
    L23 is sqrt((X3 - X2)^2 + (Y3 - Y2)^2),
    L31 is sqrt((X1 - X3)^2 + (Y1 - Y3)^2),
    D1223 is abs(L12 - L23),
    D1223<0.00000001,
    D2331 is abs(L23 - L31),
    D2331<0.00000001,
    D3112 is abs(L31 - L12),
    D3112<0.00000001.

Now, if I run that against my points above I get this:

?- points_main(Xs).
Xs = [[j1, j2, j3], [j1, j3, j2], [j2, j1, j3], [j2, j3, j1], [j3, j1, j2], [j3, j2, j1]].

All combinations of the above three points - so, yes, they are all equilateral triangles.

If I run against your original points, as expected, there are no equilateral triangles.

A few side notes.

(1) I removed all of the assert code. It wasn't helpful nor necessary to get your code working.

(2) you could define my j3 point as point(j3,0.5,X) :- X is sqrt(3)/2. and your original maths would work. However, this is just lucky. When dealing will floating-point numbers you can never be sure that two numbers are equal even though they should be.

(3) I introduced P1 \= P2, etc, to prevent the points unifying to themselves. That's why you got [a,a,a],... etc.

The reason you get odd results is that your

    ...
    point(P1,X1,Y1),
    point(P2,X2,Y2),
    point(P3,X3,Y3),
    ...

is producing the Cartesian Product of your set of points. If you had just 4 points defined, this would produced 4³ (64) possible results.

Since you are talking triangles here, where for a given 3 points A, B, and C, whether the triangle they describe is labelled ABC, ACB, BAC, BCA, CAB, CBA is irrelevant: they all describe exactly the same triangle, so what your are looking for are combinations: a selection of items from a set where order is not important, so {a,b,c] and {c,b,a] are the same combination.

In Prolog, getting combinations is easy:

combination( []    , []     ) .
combination( [H|T] , [H|T2] ) :- combination(T,T2).
combination( [_|T] ,    T2  ) :- combination(T,T2).

The only caveat is that the 2nd argument to combination/2 must be seeded as a list of the desired length. To draw combinations of 3 things from a set of 5 things is as easy as:

combination( [a,b,c,d,e] , [X,Y,Z] ).

Once you have that, getting the set of points from which to draw is easy, too, using findall/3 or setof/3:

setof( p(P,X:Y) , point(P,X,Y) , Ps ) ,

The one tricky thing is that floating point arithmetic leaves much to be desired, what with floating point jitter and all. This is how I'm computing distance:

%-----------------------------------------------------------------------
% computes the distance between two points on the Cartesian plane.
% the distance, D, is returned as an integer with an implied scale of 5,
% so the distance 1.23456789 is returned as 123457 (rounded up)
% ----------------------------------------------------------------------
distance( X1:Y1 , X2:Y2 , D ) :-
  V is sqrt( (X2-X1)^2 + (Y2-Y1)^2 ) ,
  D is round( 100000 * V )
  .

Putting it together, we get this (https://swish.swi-prolog.org/p/zHAfMTtA.pl):

equilateral_triangle( T ) :-
  setof( p(P,X:Y) , point(P,X,Y) , Ps ) ,
  combination(Ps,[A,B,C]),
  equilateral_triangle(A,B,C,T)
  .

%--------------------------------------------------------------------------------
% combination( +Set, +Combination )
%
% Takes a set (a list of distinct things) of M things and produces
% on backtracking all the combinations of M things taken N at a time.
%
% Combination, the 2nd argument MUST be initialize as a list of the desired
% length. For example, to get all the combination of 8 things taken 3 at a time,
% you'd say something like this:
%
% ?- length(C,3), combination([a,b,c,d,e],C).
%
% And get back
%
% C = [a, b, c]
% C = [a, b, d]
% C = [a, b, e]
% C = [a, c, d]
% C = [a, c, e]
% C = [a, d, e]
% C = [b, c, d]
% C = [b, c, e]
% C = [b, d, e]
% C = [c, d, e]
% 
%--------------------------------------------------------------------------------
combination( []    , []     ) .
combination( [H|T] , [H|T2] ) :- combination(T,T2).
combination( [_|T] ,    T2  ) :- combination(T,T2).

%--------------------------------------------------------------
% 3 points comprise an equilateral triangle if the triangle
% that they describe has legs of equal length
%--------------------------------------------------------------
equilateral_triangle( p(A,Pa), p(B,Pb), p(C,Pc) , t(A,B,C) ) :-
  distance(Pa,Pb,D),
  distance(Pb,Pc,D),
  distance(Pc,Pa,D).

%-----------------------------------------------------------------------
% computes the distance between two points on the Cartesian plane.
% the distance, D, is returned as an integer with an implied scale of 5,
% so the distance 1.23456789 is returned as 123457 (rounded up)
% ----------------------------------------------------------------------
distance( X1:Y1 , X2:Y2 , D ) :-
  V is sqrt( (X2-X1)^2 + (Y2-Y1)^2 ) ,
  D is round( 100000 * V )
  .

point( a , 0   , 0             ) .
point( b , 5   , 0             ) .
point( c , 2.5 , 4.33012701892 ) .
point( d , 1   , 1             ) .
point( e , 6   , 1             ) .
point( f , 3.5 , 5.33012701892 ) .