source = DirectoryName[
    ToFileName["FileName"/.NotebookInformation[EvaluationNotebook[]]]];
SetDirectory[source];
Needs["GeoBeweiser`"]
clearndg[];
Clear[u1,u2,u3,x1,x2,x3,x4];

vars = {x1, x2, x3, x4};
parameters = {u1, u2, u3};

Unprotect[C, D, N];
A = point[0, 0];
B = point[u1, 0];
C = point[u2, u3];
D = point[x1, x2];
N = point[x3, x4];

d1 = isparallel[ppline[A, B], ppline[C, D]]

d2 = isparallel[ppline[A, C], ppline[B, D]]

d3 = iscollinear[N, A, D]

d4 = iscollinear[N, B, C]

polys = {d1, d2, d3, d4}

con1 = sqrdist[A, N] - sqrdist[D, N]

con2 = sqrdist[B, N] - sqrdist[C, N]

zw1 = Union[polys, {1 - y con1}]

zw2 = Union[polys, {1 - y con2}]

GroebnerBasis[zw1, vars, CoefficientDomain -> RationalFunctions] // 
      Timing

GroebnerBasis[zw2, vars, CoefficientDomain -> RationalFunctions] // 
      Timing

dis[x_] := {x} // Transpose // Factor

h1 = x2 - u3;
h2 = (x1 - u1) u3 - u2 x2;
h3 = x4 x1 - x3 u3;
h4 = 
        x4 (u2 - u1) - (x3 - u1) u3;
g = x1^2 - 2  x1 x3 - 2  x4 x2 + x2^2;
GroebnerBasis[{h1, h2, h3, h4, 1 - y g}, {u1, u2, u3, x1, x2, x3, x4, 
        y}, MonomialOrder -> Lexicographic]

dis[
      f = GroebnerBasis[{h1, h2, h3, h4}, {x1, x2, x3, x4, u1, u2, u3}, 
          MonomialOrder -> Lexicographic]]

f[[5]] = u3

{v2 = GroebnerBasis[f, {x1, x2, x3, x4, u1, u2, u3}]} // Transpose // 
      Factor

f[[5]] = x1 - u1 - u2

{v1 = GroebnerBasis[f, {x1, x2, x3, x4, u1, u2, u3}]} // Transpose // 
      Factor

v2[[2]] = -x4;
dis[v21 = GroebnerBasis[v2, {x1, x2, x3, x4, u1, u2, u3}]]

v2[[2]] = u1 - u2;
dis[v22 = GroebnerBasis[v2, {x1, x2, x3, x4, u1, u2, u3}]]

v22[[4]] = x1;
dis[
      v221 = GroebnerBasis[v22, {x1, x2, x3, x4, u1, u2, u3}]]

v22[[4]] = x4;
dis[
      v222 = GroebnerBasis[v22, {x1, x2, x3, x4, u1, u2, u3}]]

nun V1

v1[[1]] = u1;
dis[v11 = GroebnerBasis[v1, {x1, x2, x3, x4, u1, u2, u3}]]

v1[[1]] = u3 - 2  x4;
dis[v12 = GroebnerBasis[v1, {x1, x2, x3, x4, u1, u2, u3}]]

v12[[2]] = u3;
dis[
      v121 = GroebnerBasis[v12, {x1, x2, x3, x4, u1, u2, u3}]]

v12[[2]] = u1 + u2 - 2  x3;
dis[
      v122 = GroebnerBasis[v12, {x1, x2, x3, x4, u1, u2, u3}]]

{v11, v121, v122, v21, v221, v222}

v122

g1 = -con

GroebnerBasis[{v122, 1 - y g1}, {u1, u2, u3, x1, x2, x3, 
        x4}]