Skip to content

Instantly share code, notes, and snippets.

@takuma7
Created November 20, 2014 06:08
Show Gist options
  • Save takuma7/3fbb199a02d6421d7773 to your computer and use it in GitHub Desktop.
Save takuma7/3fbb199a02d6421d7773 to your computer and use it in GitHub Desktop.
(* 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];
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment