wu.nb

(* Sort list of polynomial functions (fs) based on degree of x in ascending order *)
SortByDegree[fs_, x_] := Module[{},
	Comp[f1_,f2_] := Exponent[f1,x] < Exponent[f2,x];
	Sort[fs,Comp]
];

(* Calculate pseud-remainder *)
Rem[f_, g_, x_] := Module[{q,r,u},
	{q,r} = Together[#]&/@PolynomialQuotientRemainder[f,g,x];
	u = PolynomialLCM[Denominator[q], Denominator[r]];
	Numerator[Together[r*u]]
];

(* Triangulation *)
setAttributes[ Triangulate, HoldAll];
Triangulate[fsImm_, xsImm_] := Module[{fs=fsImm, xs=xsImm, r, n, k, ts},
	r = Length[fs];
	For[n=r, n>0, nil,
		(* Step 3 *)
		ts = fs[[1;;n]];
		ts = SortByDegree[ts, xs[[n]]];
		k = 1;
		While[k<=n && Exponent[ts[[k]], xs[[n]] ]==0, k++];
		k--;
		If[ k==n-1,
			(* Step 4 *)
			fs[[1;;n]] = ts[[1;;n]];
			n--,
			(* Step 5 *)
			fs[[1;;k+1]] = ts[[1;;k+1]];
			fs[[k+2;;n]] = Rem[ts[[#]], ts[[k+1]], xs[[n]] ]& /@ Range[k+2, n];
		];
	];
	fs
];

(* Wu's Procedure *)
Wu[hyps_, conc_, xs_] := Module[{rems, subs, r, i},
	rems = {conc};
	subs = 1;
	r = Length[hyps];
	For[i=1, i<=r, i++,
		rems = Append[rems, Rem[rems[[i]], hyps[[r-i+1]], xs[[r-i+1]] ] ];
		subs *= PolynomialQuotient[hyps[[r-i+1]], xs[[r-i+1]]^(Exponent[hyps[r-i+1], xs[[r-i+1]] ]) , xs[[r-i+1]] ]
	];
	{rems, subs}
];

(* Solve and Print *)
SolveWithWu[hyps_, conc_, xs_] := Module[{tris,rems,subs,i},
	tris = Triangulate[hyps,xs];
	{rems, subs} = Wu[hyps, conc, xs];
	Do[Print[StringFom["Hyp``\t= ``", i, hyps[[i]] ]], {i,1,Length[hyps]}];
	Print[""];
	Do[Print[StringForm["Tri``\t= ``", i, tris[[i]] ]], {i,1,Length[tris]}];
	Print[StringForm["Conc\t= ``", conc]];
	Print[""];
	Do[Print[StringForm["Rem``\t= ``", i, rems[[i]] ]], {i,1,Length[rems]}];
	Print[StringForm["Subsidiaries = ``", subs]];
];

hyps = {ax^2+ay^2-bx^2-by^2, ax^2+ay^2-cx^2-cy^2, ax^2+ay^2-dx^2-dy^2, (fx-cx)(cy-by)-(fy-cy)(cx-bx), (fx-dx)(dy-ay)-(fy-dy)(dx-ax), (ex-bx)(by-ay)-(ey-by)(bx-ax), (ex-cx)(cy-dy)-(ey-cy)(cx-dx), gx-(ax+cx)/2, gy-(ay+cy)/2, hx-(bx+dx)/2, hy-(by+dy)/2};
conc = (fx-hx)(gx-ex)+(fy-hy)(gy-ey);
xs = {bx, cx, dx, ex, ey, fx, fy, gx, gy, hx, hy};
SolveWithWu[hyps, conc, xs];