You try to write the command restart at the begin your document. 

> restart;

L:=[]:

for x1 from -5 to 5 do

for y1 from -5 to 5 do

for z1 from -5 to 5 do

for x2 from -5 to 5 do

for y2 from -5 to 5 do

for z2 from -5 to 5 do

M:=[1,-2,3]:

a:=(x1-M[1])^2+(y1-M[2])^2+(z1-M[3])^2:

if a>0 and a=(x2-M[1])^2+(y2-M[2])^2+(z2-M[3])^2 and a=(x2-x1)^2+(y2-y1)^2+(z2-z1)^2 then

L:=[op(L), [[M[1],M[2],M[3]], [x1,y1,z1], [x2,y2,z2]]] fi;

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

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

> with(geom3d):

L[10];

point(A,op(1,L[10])):

point(B,op(2,L[10])):

point(C,op(3,L[10])):

triangle(ABC,[A,B,C]):

IsEquilateral(ABC);

I don't understand, in the code

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];

i replaced [0,0,0] in the line L:=[op(L), [[0,0,0], [x1,y1,z1], [x2,y2,z2]]] fi; by [1,-2,-3], i received a triangle, for example A(1,-2,-3) , B(-10,-7,1), C(-5,-5,-10). But, this triangle is not a equilateral. Please help me. Thank you.

If i do as you said, then i don't get my required. This is my solution.

First step, 

> restart:

with(geom3d):

point(A,2,3,3):

point(B,1,1,1):

R:=distance(A,B):

Equation(sphere(S,[A,3*R],[x,y,z]));

f := (a,b,c) -> eval(lhs(Equation(S)), [x=a,y=b,z=c]);

Npoints:=proc(N, eq, m)

local L, i, a, b;

L:=[];

for i while nops(L)<N do

a:=rand(-m..m)(); b:=rand(-m..m)();

L:={op(L),[a, b, solve(eval(eq, [x=a, y=b]), z)]};

od;

L;

end proc:

Example:

Npoints(10, Equation(S),5);

next, i see in Npoints(10, Equation(S),5); and  choose solution (a, b, c, d) with a, b, c, d are integer numbers. 

Second step

> start:

with(LinearAlgebra):

A:=<2,3,3>:

B:=<1,1,1>:

C:=<6,1,7>:

CrossProduct(B-A,C-B);

M:=<x,y,z>:

AB:=Norm(B-A,2):

AC:=Norm(C-A,2):

k:=AB/AC:

eq1:=solve([seq((B-M)[i]= -k*(C-M)[i],i=1..3)],{x,y,z}): #bisector of the interior

T:=subs(eq1,M);

Norm(A-T,2);

eq2:=solve([seq((B-M)[i]= k*(C-M)[i],i=1..3)],{x,y,z}):

#bisector of the exterior

E:=subs(eq2,M);

Norm(A-E,2);

My computer runs too long. 

Thank you very much.

L:=[ ]:

for d from 1 to 100 do

for a from 1 to floor(d/sqrt(3)) do

for b from a to floor(sqrt(d^2-2*a^2)) do

for c from b to floor(sqrt(d^2-a^2-b^2)) do

if a^2+b^2+c^2=d^2 then L:=[op(L), [a, b, c, d]]; fi;

od: od: od: od:

If four numbers a, b, c, d are different from 0. How  i do? Please help me.

Thank you very much.

Please comment when i put point(M,a): and point(M,subs(sol,a)): 

I repaired my code

> restart:

with(geom3d):

point(A,1,-1,2):

a:=[2*t-1,t,t+2]:

point(M,a):

reflection(N,M,A):

T:=coordinates(N):

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

eq:=subs(x=T[1],y=T[2],z=T[3], Equation(P)):

sol:=solve(eq,{t}):

point(M,subs(sol,a)):

Equation(line(AM,[A,M],t));

Thank you very much. Because, i have a problem: Let A(1,-1,2) be a point, d,[2*t-1,t,t+2] be a line and (P): x+y-2*z+5=0 be a plane. Write the equation of the straigh line passing through the point A and cut the line d at M, the plane (P) at M so that A is midpoint of the segment MN.

This is my code, before i receive your answers

restart;with(geom3d):

point(A,1,-1,2):

line(d,[2*t-1,t,t+2],t):

point(M,op(Equation(d))):

reflection(N,M,A):

T:=coordinates(N):

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

f:=(a,b,c)->eval(lhs(Equation(P)), [x=a,y=b,z=c]): 

sig:=f(T[1],T[2],T[3]):

sol:=solve(sig=0,{t}):

point(M,subs(sol,coordinates(M))):

Equation(line(MA,[M,A],t));

Please comment to me. Thank you very much.

OK, i can answer my question. 

line(AM,3*x+5*y-4=0,[x,y]);

y:=unapply(solve(Equation(AM),y),x);

point(M,a,y(a));

coordinates(M);

This is my problem. Let ABC is a triangle and 

AM: 7*x-5*y=0, BD: x-2*y-30=0, CK: x-y+16=0 are the respective  equations of median, angle bisector, altitude drawn from A, B, C. Find the coordinates of vertices A, B, C of the triangle ABC.

My idea.

1. Let M be midpoint of the BC, M(a,7*a/5);
2. Let N be symmetric to the point M with respect to the point M;
3. Write the equation of the line AB passing through the point N and perpendicular to the line CK;
4. Find the point of intersection B of two line BD and AB;
5. Let C be symmetric to the point M with respect to the point M.

If the equation of the median AM is a a different equation, i must change the coordinates of the point M

For example, i have the line AM has equation 7*x-5*y=0. Now i want to choose the point on the line AM. 

line(AM,7*x-5*y=0,[x,y]):

point(M,a,7*a/5)

But, with this way, when i change equation the line AM, then i must change the coordinates the point M by hand. I want when

i change equation the line AM, i do not change coordinates the point M.

This is my code. 

> restart;with(geometry):

line(AM,7*x-5*y=0,[x,y]):

line(BD,x-2*y-30=0,[x,y]):

line(CK,x-y+16=0, [x,y]):

point(M,a,7*a/5): 

reflection(N,M,BD):

coordinates(N);

PerpendicularLine(AB,N,CK):

Equation(AB,[x,y]);

intersection(B,AB,BD):

coordinates(B);

reflection(C,B,M):

coordinates(C);

f := (a,b) -> eval(lhs(Equation(CK)), [x=a,y=b]):

sol:=solve(f(HorizontalCoord(C),VerticalCoord(C))=0,{a});

subs(op(1,sol),coordinates(B));

point('B',subs(op(1,sol),coordinates(B)));

subs(op(1,sol),coordinates(C));

PerpendicularLine(BN,B,CK);

Equation(BN,[x,y]);

intersection(A,AM,BN);

coordinates(A);