1) If you want to [0,0,0] always been in the first position, then replace the braces on the brackets.

2) It's an obvious Maple bug.

Your main mistake: you look for only in the first octant in which there are no such triangles, as if two vectors out of the origin and lie entirely in the first octant, the angle between them does not exceed 90 degrees.

 Write code as follows:

N:=5:

L:=[]:

for x1 from -N to N do

for y1 from x1 to N do

for z1 from y1 to N do

for x2 from -N to N do

for y2 from -N to N do

for z2 from -N to N do

a:=x1^2+y1^2+z1^2:

b:=x2^2+y2^2+z2^2:

c:=a+b-2*sqrt(a*b)*cos(2*Pi/3):

if type(a, positive) and type(b, positive) and c=(x2-x1)^2 + (y2-y1)^2 + (z2-z1)^2 then L:=[op(L), {[0, 0, 0], [x1, y1, z1], [x2, y2, z2]}]: fi:

od: od: od: od: od: od:

nops(L);

L;

Code in Maple 13:

line1:=20-0.05*x:

line2:=2+0.0002*x^2:

A:=plot([line1,line2], x=0..450, color=[red, blue], thickness=3):

B:=plots[implicitplot]((y-line1)*(y-line2),x=0..450, y=-5..50, coloring=[yellow, white], filledregions=true, numpoints=50000):

plots[display](A, B);

v:=Vector([5, -1, 6, 3]):
op(1, v);

4

Trivial examples:  Pythagorean triangles [[0, 0, 0], [2ab, 0, 0], [0, 0, a^2-b^2]] , where a and b - any unequal positive integers.

Nontrivial examples can be found using a simple code:

N:=10:

L:=[]:

for x1 to N do

for y1 from x1 to N do

for z1 from y1 to N do

for x2 to N do

for y2 to N do

for z2 to N do

a:=sqrt((x2-x1)^2+(y2-y1)^2+(z2-z1)^2):

if type(a, positive) and type(sqrt(x1^2+y1^2+z1^2)+ sqrt(x2^2+y2^2+z2^2)+a, integer) then L:=[op(L), {[0, 0, 0], [x1, y1, z1], [x2, y2, z2]}]: fi:

od: od: od: od: od: od:

nops(L);

L;

Yes, of course! See help for the  addcoords  command.

Your code can be significantly improved:

1) Since always  x2=x1+x3 and y2=y1+y3 and z2=z1+z3 , then variables x2, y2, z2  should not be written as a variables in a loop. This will speed up program execution.

2) You can exclude the  different trivial variants. If the points (x1, y1, z1) and (x3, y3, z3) lie in a coordinate plane, eg, xOy, then it will always be

(x1, y1, z1) = (x, y, 0) , (x3, y3, z3) = (-y, x, 0)  or  (x1, y1, z1) = (x, y, 0) , (x3, y3, z3) = (y, -x, 0)

(x and y - any numbers)

3) You can also foreclose variants, which are obtained from one another by simple transformations, such as reflections in the coordinate planes and axes.

4) Lists  [[0,0,0], [x1, y1, z1], [x2, y2, z2], [x3, y3, z3]]  and [[0,0,0], [x3, y3, z3], [x2, y2, z2], [x1, y1, z1]] give the same square, so from two lists, you can leave one list.

To implement  3) and 4) can be used  ListTools[Categorize]  command.

Write your 5 variables   (a,b,c,d,e)  in the list.

Example:

Case:= [ [1, 2, 3, 1.5, 6], [3, 3, 1, 3, 8], [5, 3, 3, 5, 3] ]:
Case[1], Case[2], Case[3];

       [1, 2, 3, 1.5, 6], [3, 3, 1, 3, 8], [5, 3, 3, 5, 3]

This problem can be solved similar to your problem of an equilateral triangle. Take advantage of the fact that if [[0,0,0], [x1, y1, z1], [x2, y2, z2], [x3, y3, z3]] are the vertices of a square, then x2=x1+x3, y2=y1+y3, z2=z1+z3. It is therefore sufficient to seek only two vertices [x1, y1, z1], [x3, y3, z3] .

Let us the solution to be non trivial if none of the vertices of other than [0, 0, 0] , does not lie in the coordinate plane, and the absolute values ​​of the coordinates of each vertex are unequal and integer numbers. The simplest nontrivial solution is  [[0, 0, 0], [6, 3, -2], [9, 1, 4], [3, -2, 6]] . Interestingly, the length of the side of the square is also an integer equal to 7.

If I understand your problem, it is simple. Should be in the loop to multiply all the elements of each row to the least common multiple of their denominators.

Code of the procedure:

DoIntegers:=proc(A)

local m, n, i, j, d, B, f;

uses LinearAlgebra;

m:=RowDimension(A);

n:=ColumnDimension(A);

for i to m do

d:=ilcm(seq(denom(A[i,j]), j=1..n));

for j to n do

B[i,j]:=A[i,j]*d;

od; od;

f:=(i,j)->B[i,j];

B:=Matrix(m, n, f);

end proc:

Example:

A:=Matrix([[1, 1/2, 1/3], [1/4, 1/5, 1/6], [7, 8, 9]]);

B=DoIntegers(A);

1) Your code is correct.

2) Since the triangle with the vertex at an arbitrary point can be obtained from the triangle with the vertex at the origin of the usual translation, there is no need to write new code.

3) Rewrote my previous code. Now the triangles, differing only in the numbering of the vertices are not distinct. In this code, the user himself determines the required number of triangles (parameter N). The search range is not fixed from the beginning, but gradually expands as new triangles. Therefore, if the number N is not too large, then the procedure is fast.

N_equil_triangles:=proc(N)

local x1, y1, z1, x2, y2, z2, k, n, a;

global L;

k:=0:

L:=[]:

for n while k<N do

for x1 from -n to n while k<N do

for y1 from -n to n while k<N do

for z1 from -n to n while k<N do

for x2 from -n to n while k<N do

for y2 from -n to n while k<N do

for z2 from -n to n while k<N do

a:=x1^2+y1^2+z1^2:

if a>0 and a=x2^2+y2^2+z2^2 and a=(x2-x1)^2+(y2-y1)^2+(z2-z1)^2 then

L:=[op(L),{[0,0,0], [x1,y1,z1], [x2,y2,z2]}]: k:=nops(L): fi:

od: od: od: od: od: od: od:

end proc:

Example:

t:=time():

N_equil_triangles(1000):

L[500..510];

Time=cat(convert(time()-t, symbol),`  sec`);

I'm not sure that all colors are identified correctly. Details can be improved:

R:=3: r:=1: h:=0.25:  # h specifies the width of the band

Colors:=["Red", "HotPink",  "Blue", "BlueViolet", "DeepSkyBlue", "LightBlue", "LightSteelBlue", "LimeGreen", "Lime", "Yellow", "Orange", "Gold"]:

E:=seq(seq(plottools[line]([(R+r*cos(11/3*t))*cos(t), (R+r*cos(11/3*t))*sin(t), r*sin(11/3*t)],[(R+r*cos(11/3*t))*cos(t+h), (R+r*cos(11/3*t))*sin(t+h), r*sin(11/3*t)], color=Colors[k+1], thickness=5), t=6*k/12*3.14..6*(k+1)/12*3.14, 0.01), k=0..11):

plots[display]( E, scaling=constrained);

Points A and B are taken arbitrarily, and point C is taken so that  vector C-A is 3 times longer than  B-A.

The set of points in space equidistant from two given planes is the union of two perpendicular planes, each of which bisects the corresponding dihedral angle between the given planes. In fact you have found this set. Its equation is E . To select from it each of the planes you can as follows:

with(geom3d):

plane(P,2*x+3*y-6*z+1=0, [x, y, z]):

plane(Q,x -4*y-8*z+1=0, [x,y,z]):

point(M, x, y, z):

E:=distance(M,P)= distance(M,Q);

simplify(E) assuming op([2,1], lhs(E))>0 and op([2,1], rhs(E))>0 ;  # equation of the first plane

simplify(E) assuming op([2,1], lhs(E))>0 and op([2,1], rhs(E))<0 ;  # equation of the second plane

If there is no  DirectSearch package, then the problem can be easily solved by exhaustive search . To reduce the amount of enumeration, we assume that one vertex locates at the origin:

restart;

L:=[]:

for x1 from -10 to 10 do

for y1 from -10 to 10 do

for z1 from -10 to 10 do

for x2 from -10 to 10 do

for y2 from -10 to 10 do

for z2 from -10 to 10 do

a:=x1^2+y1^2+z1^2:

if a>0 and a=x2^2+y2^2+z2^2 and a=(x2-x1)^2+(y2-y1)^2+(z2-z1)^2 then

L:=[op(L), [[0,0,0], [x1,y1,z1], [x2,y2,z2]]] fi;

od: od: od: od: od: od:

nops(L);  L[1..20];

We found 3936 solutions and displayed on the screen the first 20 decisions.

Let compare the solution obtained with the previous one by DirectSearch package. For correctness  we put one vertex at the origin (assume [х3, y3, z3]=[0, 0, 0]). Also we find the total number of solutions and display the first 20 ones:

restart;

interface(rtablesize=infinity):

L:=DirectSearch:-SolveEquations([x1^2+y1^2+z1^2 = (x2-x1)^2+(y2-y1)^2+(z2-z1)^2,

x1^2+y1^2+z1^2 = x2^2+y2^2+z2^2],  {x1^2+y1^2+z1^2 > 0, abs(x1) <= 10,

abs(y1) <= 10, abs(z1)<=10, abs(x2) <= 10, abs(y2) <= 10, abs(z2) <= 10}, assume = integer, AllSolutions):

LinearAlgebra[RowDimension](L);

seq(L[k,3], k=1..20);

We see that the program found only 98 solutions from 3936 solutions .