Sun Aug 18 20:49:37 2002 run on Windows
Geometry 1.1 Last update Sept 6, 1998
% Author H.-G. Graebe | Univ. Leipzig | Version 6.9.1998
% graebe@informatik.uni-leipzig.de
comment
Test suite for the package GEOMETRY 1.1
end comment;
algebraic;
load cali,geometry;
off nat;
on echo;
showtime;
Time: 10 ms plus GC time: 200 ms
% #####################
% Some one line proofs
% #####################
% A generic triangle ABC
A:=Point(a1,a2);
a := {a1,a2}$
B:=Point(b1,b2);
b := {b1,b2}$
C:=Point(c1,c2);
c := {c1,c2}$
% Its midpoint perpendiculars have a point in common:
concurrent(mp(a,b),mp(b,c),mp(c,a));
0$
% This point
M:=intersection_point(mp(a,b),mp(b,c));
m := {(a1**2*b2 - a1**2*c2 + a2**2*b2 - a2**2*c2 - a2*b1**2 - a2*b2**2 + a2*c1**
2 + a2*c2**2 + b1**2*c2 + b2**2*c2 - b2*c1**2 - b2*c2**2)/(2*(a1*b2 - a1*c2 - a2
*b1 + a2*c1 + b1*c2 - b2*c1)),
( - a1**2*b1 + a1**2*c1 + a1*b1**2 + a1*b2**2 - a1*c1**2 - a1*c2**2 - a2**2*b1 +
a2**2*c1 - b1**2*c1 + b1*c1**2 + b1*c2**2 - b2**2*c1)/(2*(a1*b2 - a1*c2 - a2*b1
+ a2*c1 + b1*c2 - b2*c1))}$
% is the center of the circumscribed circle
sqrdist(M,A) - sqrdist(M,B);
0$
% The altitutes intersection theorem
concurrent(altitude(a,b,c),altitude(b,c,a),altitude(c,a,b));
0$
% The median intersection theorem
concurrent(median(a,b,c),median(b,c,a),median(c,a,b));
0$
% Euler's line
M:=intersection_point(mp(a,b),mp(b,c));
m := {(a1**2*b2 - a1**2*c2 + a2**2*b2 - a2**2*c2 - a2*b1**2 - a2*b2**2 + a2*c1**
2 + a2*c2**2 + b1**2*c2 + b2**2*c2 - b2*c1**2 - b2*c2**2)/(2*(a1*b2 - a1*c2 - a2
*b1 + a2*c1 + b1*c2 - b2*c1)),
( - a1**2*b1 + a1**2*c1 + a1*b1**2 + a1*b2**2 - a1*c1**2 - a1*c2**2 - a2**2*b1 +
a2**2*c1 - b1**2*c1 + b1*c1**2 + b1*c2**2 - b2**2*c1)/(2*(a1*b2 - a1*c2 - a2*b1
+ a2*c1 + b1*c2 - b2*c1))}$
H:=intersection_point(altitude(a,b,c),altitude(b,c,a));
h := {( - a1*a2*b1 + a1*a2*c1 + a1*b1*b2 - a1*c1*c2 - a2**2*b2 + a2**2*c2 + a2*
b2**2 - a2*c2**2 - b1*b2*c1 + b1*c1*c2 - b2**2*c2 + b2*c2**2)/(a1*b2 - a1*c2 -
a2*b1 + a2*c1 + b1*c2 - b2*c1),
(a1**2*b1 - a1**2*c1 + a1*a2*b2 - a1*a2*c2 - a1*b1**2 + a1*c1**2 - a2*b1*b2 + a2
*c1*c2 + b1**2*c1 + b1*b2*c2 - b1*c1**2 - b2*c1*c2)/(a1*b2 - a1*c2 - a2*b1 + a2*
c1 + b1*c2 - b2*c1)}$
S:=intersection_point(median(a,b,c),median(b,c,a));
s := {(a1 + b1 + c1)/3,(a2 + b2 + c2)/3}$
collinear(M,H,S);
0$
sqrdist(S,varpoint(M,H,2/3));
0$
% Feuerbach's circle
% Choose a special coordinate system
A:=Point(0,0);
a := {0,0}$
B:=Point(u1,0);
b := {u1,0}$
C:=Point(u2,u3);
c := {u2,u3}$
M:=intersection_point(mp(a,b),mp(b,c));
m := {u1/2,( - u1*u2 + u2**2 + u3**2)/(2*u3)}$
H:=intersection_point(altitude(a,b,c),altitude(b,c,a));
h := {u2,(u2*(u1 - u2))/u3}$
N:=midpoint(M,H);
n := {(u1 + 2*u2)/4,(u1*u2 - u2**2 + u3**2)/(4*u3)}$
sqrdist(N,midpoint(A,B))-sqrdist(N,midpoint(B,C));
0$
sqrdist(N,midpoint(A,B))-sqrdist(N,midpoint(H,C));
0$
D:=intersection_point(pp_line(A,B),pp_line(H,C));
d := {u2,0}$
sqrdist(N,midpoint(A,B))-sqrdist(N,D);
0$
clear_ndg();
{}$
clear(A,B,C,D,M,H,S,N);
% #############################
% Non-linear Geometric Objects
% #############################
% Bisector intersection theorem
A:=Point(0,0);
a := {0,0}$
B:=Point(1,0);
b := {1,0}$
C:=Point(u1,u2);
c := {u1,u2}$
P:=Point(x1,x2);
p := {x1,x2}$
polys:={
point_on_bisector(P,A,B,C),
point_on_bisector(P,B,C,A),
point_on_bisector(P,C,A,B)};
polys := { - 2*u1*x1*x2 + 2*u1*x2 + u2*x1**2 - 2*u2*x1 - u2*x2**2 + u2 + 2*x1*x2
- 2*x2,
- 2*u1**3*x2 + 2*u1**2*u2*x1 - u1**2*u2 + 2*u1**2*x1*x2 + 2*u1**2*x2 - 2*u1*u2
**2*x2 - 2*u1*u2*x1**2 + 2*u1*u2*x2**2 - 2*u1*x1*x2 + 2*u2**3*x1 - u2**3 - 2*u2
**2*x1*x2 + 2*u2**2*x2 + u2*x1**2 - u2*x2**2,
2*u1*x1*x2 - u2*x1**2 + u2*x2**2}$
con1:=num(sqrdist(P,pedalpoint(p,pp_line(A,C)))-x2^2);
con1 := u2*( - 2*u1**3*x1*x2 + u1**2*u2*x1**2 - u1**2*u2*x2**2 - 2*u1*u2**2*x1*
x2 + u2**3*x1**2 - u2**3*x2**2)$
con2:=num(sqrdist(p,pedalpoint(p,pp_line(B,C)))-x2^2);
con2 := u2*( - 2*u1**3*x1*x2 + 2*u1**3*x2 + u1**2*u2*x1**2 - 2*u1**2*u2*x1 - u1
**2*u2*x2**2 + u1**2*u2 + 6*u1**2*x1*x2 - 6*u1**2*x2 - 2*u1*u2**2*x1*x2 + 2*u1*
u2**2*x2 - 2*u1*u2*x1**2 + 4*u1*u2*x1 + 2*u1*u2*x2**2 - 2*u1*u2 - 6*u1*x1*x2 + 6
*u1*x2 + u2**3*x1**2 - 2*u2**3*x1 - u2**3*x2**2 + u2**3 + 2*u2**2*x1*x2 - 2*u2**
2*x2 + u2*x1**2 - 2*u2*x1 - u2*x2**2 + u2 + 2*x1*x2 - 2*x2)$
setring({x1,x2},{},lex);
{{x1,x2},{},lex,{1,1}}$
setideal(polys,polys);
{u2*x1**2 - (2*u1 - 2)*x1*x2 - (2*u2)*x1 - u2*x2**2 + (2*u1 - 2)*x2 + u2,
- (2*u1*u2 - u2)*x1**2 + (2*u1**2 - 2*u1 - 2*u2**2)*x1*x2 + (2*u1**2*u2 + 2*u2
**3)*x1 + (2*u1*u2 - u2)*x2**2 - (2*u1**3 - 2*u1**2 + 2*u1*u2**2 - 2*u2**2)*x2 -
(u1**2*u2 + u2**3),
- u2*x1**2 + (2*u1)*x1*x2 + u2*x2**2}$
gbasis polys;
{(4*u2)*x2**4 - (8*u1**2 - 8*u1 + 8*u2**2)*x2**3 + (4*u1**2*u2 - 4*u1*u2 + 4*u2
**3 - 4*u2)*x2**2 + (4*u2**2)*x2 - u2**3,
(2*u1*u2**2 - u2**2)*x1 + (2*u2)*x2**3 - (4*u1**2 - 4*u1 + 2*u2**2)*x2**2 - (2*
u1**2*u2 - 2*u1*u2 + 2*u2)*x2 - (u1*u2**2 - u2**2)}$
{con1,con2} mod gbasis polys;
{0,0}$
% Bisector intersection theorem. A constructive proof.
A:=Point(0,0);
a := {0,0}$
B:=Point(1,0);
b := {1,0}$
P:=Point(u1,u2);
p := {u1,u2}$
l1:=pp_line(A,B);
l1 := {0,-1,0}$
l2:=symline(l1,pp_line(A,P));
l2 := { - 2*u1*u2,u1**2 - u2**2,0}$
l3:=symline(l1,pp_line(B,P));
l3 := {2*u2*( - u1 + 1),
u1**2 - 2*u1 - u2**2 + 1,
2*u2*(u1 - 1)}$
point_on_bisector(P,A,B,intersection_point(l2,l3));
0$
clear_ndg();
{}$
clear(A,B,C,P,l1,l2,l3);
% Miquel's theorem
on gcd;
A:=Point(0,0);
a := {0,0}$
B:=Point(1,0);
b := {1,0}$
C:=Point(c1,c2);
c := {c1,c2}$
P:=choose_pl(pp_line(A,B),u1);
p := {u1,0}$
Q:=choose_pl(pp_line(B,C),u2);
q := {u2,(c2*(u2 - 1))/(c1 - 1)}$
R:=choose_pl(pp_line(A,C),u3);
r := {u3,(c2*u3)/c1}$
X:=other_cc_point(P,p3_circle(A,P,R),p3_circle(B,P,Q))$
point_on_circle(X,p3_circle(C,Q,R));
0$
off gcd;
clear_ndg();
{}$
clear(A,B,C,P,Q,R,X);
% ########################
% Theorems of linear type
% ########################
% Pappus' theorem
A:=Point(u1,u2);
a := {u1,u2}$
B:=Point(u3,u4);
b := {u3,u4}$
C:=Point(x1,u5);
c := {x1,u5}$
P:=Point(u6,u7);
p := {u6,u7}$
Q:=Point(u8,u9);
q := {u8,u9}$
R:=Point(u0,x2);
r := {u0,x2}$
polys:={collinear(A,B,C), collinear(P,Q,R)};
polys := {u1*u4 - u1*u5 - u2*u3 + u2*x1 + u3*u5 - u4*x1,
u0*u7 - u0*u9 + u6*u9 - u6*x2 - u7*u8 + u8*x2}$
con:=collinear(
intersection_point(pp_line(A,Q),pp_line(P,B)),
intersection_point(pp_line(A,R),pp_line(P,C)),
intersection_point(pp_line(B,R),pp_line(Q,C)))$
vars:={x1,x2};
vars := {x1,x2}$
sol:=solve(polys,vars);
sol := {{x1=( - u1*u4 + u1*u5 + u2*u3 - u3*u5)/(u2 - u4),
x2=(u0*u7 - u0*u9 + u6*u9 - u7*u8)/(u6 - u8)}}$
sub(sol,con);
0$
% Pappus' theorem. A constructive approach
A:=Point(u1,u2);
a := {u1,u2}$
B:=Point(u3,u4);
b := {u3,u4}$
P:=Point(u6,u7);
p := {u6,u7}$
Q:=Point(u8,u9);
q := {u8,u9}$
C:=choose_pl(pp_line(A,B),u5);
c := {u5,
(u1*u4 - u2*u3 + u2*u5 - u4*u5)/(u1 - u3)}$
R:=choose_pl(pp_line(P,Q),u0);
r := {u0,
(u0*u7 - u0*u9 + u6*u9 - u7*u8)/(u6 - u8)}$
con:=collinear(intersection_point(pp_line(A,Q),pp_line(P,B)),
intersection_point(pp_line(A,R),pp_line(P,C)),
intersection_point(pp_line(B,R),pp_line(Q,C)));
con := 0$
clear_ndg();
{}$
clear(A,B,C,P,Q,R);
% ###########################
% Theorems of non linear type
% ###########################
% Fermat Point
A:=Point(0,0);
a := {0,0}$
B:=Point(0,2);
b := {0,2}$
C:=Point(u1,u2);
c := {u1,u2}$
P:=Point(x1,x2);
p := {x1,x2}$
Q:=Point(x3,x4);
q := {x3,x4}$
R:=Point(x5,x6);
r := {x5,x6}$
polys1:={sqrdist(P,B)-sqrdist(B,C), sqrdist(P,C)-sqrdist(B,C),
sqrdist(Q,A)-sqrdist(A,C), sqrdist(Q,C)-sqrdist(A,C),
sqrdist(R,B)-sqrdist(A,B), sqrdist(R,A)-sqrdist(A,B)};
polys1 := { - u1**2 - u2**2 + 4*u2 + x1**2 + x2**2 - 4*x2,
- 2*u1*x1 - 2*u2*x2 + 4*u2 + x1**2 + x2**2 - 4,
- u1**2 - u2**2 + x3**2 + x4**2,
- 2*u1*x3 - 2*u2*x4 + x3**2 + x4**2,
x5**2 + x6**2 - 4*x6,
x5**2 + x6**2 - 4}$
con:=concurrent(pp_line(A,P), pp_line(B,Q), pp_line(C,R));
con := - u1*x1*x4*x6 + 2*u1*x1*x6 + u1*x2*x3*x6 - 2*u1*x2*x3 + 2*u2*x1*x3 + u2*
x1*x4*x5 - 2*u2*x1*x5 - u2*x2*x3*x5 - 2*x1*x3*x6 + 2*x2*x3*x5$
vars:={x1,x2,x3,x4,x5,x6};
vars := {x1,
x2,
x3,
x4,
x5,
x6}$
setring(vars,{},lex);
{{x1,x2,x3,x4,x5,x6},{},lex,{1,1,1,1,1,1}}$
iso:=isolatedprimes polys1;
iso := {{x5**2 - 3,
x6 - 1,
u1*x5 - u2 + 2*x4,
- u1*x5 - u2 + 2*x2 - 2,
- u1 - u2*x5 + 2*x3,
- u1 + u2*x5 + 2*x1 - 2*x5},
{x5**2 - 3,
x6 - 1,
- u1*x5 - u2 + 2*x4,
u1*x5 - u2 + 2*x2 - 2,
- u1 + u2*x5 + 2*x3,
- u1 - u2*x5 + 2*x1 + 2*x5},
{x5**2 - 3,
x6 - 1,
u1*x5 - u2 + 2*x4,
u1*x5 - u2 + 2*x2 - 2,
- u1 - u2*x5 + 2*x3,
- u1 - u2*x5 + 2*x1 + 2*x5},
{x5**2 - 3,
x6 - 1,
- u1*x5 - u2 + 2*x4,
- u1*x5 - u2 + 2*x2 - 2,
- u1 + u2*x5 + 2*x3,
- u1 + u2*x5 + 2*x1 - 2*x5}}$
for each u in iso collect con mod u;
{ - 3*u1**2*u2 + 3*u1**2 - 2*u1*u2*x5 + 2*u1*x5 - 3*u2**3 + 9*u2**2 - 6*u2,
0,
(u1**3*x5 + 3*u1**2*u2 - 6*u1**2 + u1*u2**2*x5 - 4*u1*u2*x5 + 3*u2**3 - 18*u2**2
+ 24*u2)/2,
( - u1**3*x5 + 3*u1**2*u2 - u1*u2**2*x5 + 4*u1*x5 + 3*u2**3 - 12*u2)/2}$
polys2:={sqrdist(P,B)-sqrdist(P,C),
sqrdist(Q,A)-sqrdist(Q,C),
sqrdist(R,A)-sqrdist(R,B),
num(p3_angle(R,A,B)-p3_angle(P,B,C)),
num(p3_angle(Q,C,A)-p3_angle(P,B,C))};
polys2 := { - u1**2 + 2*u1*x1 - u2**2 + 2*u2*x2 - 4*x2 + 4,
- u1**2 + 2*u1*x3 - u2**2 + 2*u2*x4,
4*(x6 - 1),
- u1*x1*x5 - u1*x2*x6 + 2*u1*x6 + u2*x1*x6 - u2*x2*x5 + 2*u2*x5 - 2*x1*x6 + 2*
x2*x5 - 4*x5,
u1**3*x2 - 2*u1**3 - u1**2*u2*x1 + u1**2*x1*x4 + 2*u1**2*x1 - u1**2*x2*x3 + 2*u1
**2*x3 + u1*u2**2*x2 - 2*u1*u2**2 - 2*u1*x1*x3 - 2*u1*x2*x4 + 4*u1*x4 - u2**3*x1
+ u2**2*x1*x4 + 2*u2**2*x1 - u2**2*x2*x3 + 2*u2**2*x3 - 2*u2*x1*x4 + 2*u2*x2*x3
- 4*u2*x3}$
sol:=solve(polys2,{x1,x2,x3,x4,x6});
sol := {{x1=(u1 + u2*x5 - 2*x5)/2,
x2=( - u1*x5 + u2 + 2)/2,
x3=(u1 - u2*x5)/2,
x4=(u1*x5 + u2)/2,
x6=1}}$
sub(sol,con);
0$
clear_ndg();
{}$
clear(A,B,C,P,Q,R);
% ####################
% Desargue's theorem
% ####################
% A constructive proof.
A:=Point(a1,a2);
a := {a1,a2}$
B:=Point(b1,b2);
b := {b1,b2}$
C:=Point(c1,c2);
c := {c1,c2}$
R:=Point(d1,d2);
r := {d1,d2}$
S:=choose_pl(par(R,pp_line(A,B)),u);
s := {u,
(a1*d2 - a2*d1 + a2*u - b1*d2 + b2*d1 - b2*u)/(a1 - b1)}$
T:=intersection_point(par(R,pp_line(A,C)),par(S,pp_line(B,C)));
t := {(a1*u - b1*d1 + c1*d1 - c1*u)/(a1 - b1),
(a1*d2 - a2*d1 + a2*u - b1*d2 + c2*d1 - c2*u)/(a1 - b1)}$
con:=concurrent(pp_line(A,R),pp_line(B,S),pp_line(C,T));
con := 0$
% Desargue's theorem as theorem of linear type.
A:=Point(u1,u2);
a := {u1,u2}$
B:=Point(u3,u4);
b := {u3,u4}$
C:=Point(u5,u6);
c := {u5,u6}$
R:=Point(u7,u8);
r := {u7,u8}$
S:=Point(u9,x1);
s := {u9,x1}$
T:=Point(x2,x3);
t := {x2,x3}$
polys:={parallel(pp_line(R,S),pp_line(A,B)),
parallel(pp_line(S,T),pp_line(B,C)),
parallel(pp_line(R,T),pp_line(A,C))};
polys := { - u1*u8 + u1*x1 + u2*u7 - u2*u9 + u3*u8 - u3*x1 - u4*u7 + u4*u9,
- u3*x1 + u3*x3 + u4*u9 - u4*x2 + u5*x1 - u5*x3 - u6*u9 + u6*x2,
- u1*u8 + u1*x3 + u2*u7 - u2*x2 + u5*u8 - u5*x3 - u6*u7 + u6*x2}$
con:=concurrent(pp_line(A,R),pp_line(B,S),pp_line(C,T));
con := - u1*u3*u6*u8 + u1*u3*u6*x1 + u1*u3*u8*x3 - u1*u3*x1*x3 + u1*u4*u5*u8 -
u1*u4*u5*x3 - u1*u4*u6*u9 + u1*u4*u6*x2 - u1*u4*u8*x2 + u1*u4*u9*x3 - u1*u5*u8*
x1 + u1*u5*x1*x3 + u1*u6*u8*u9 - u1*u6*x1*x2 - u1*u8*u9*x3 + u1*u8*x1*x2 - u2*u3
*u5*x1 + u2*u3*u5*x3 + u2*u3*u6*u7 - u2*u3*u6*x2 - u2*u3*u7*x3 + u2*u3*x1*x2 -
u2*u4*u5*u7 + u2*u4*u5*u9 + u2*u4*u7*x2 - u2*u4*u9*x2 + u2*u5*u7*x1 - u2*u5*u9*
x3 - u2*u6*u7*u9 + u2*u6*u9*x2 + u2*u7*u9*x3 - u2*u7*x1*x2 + u3*u5*u8*x1 - u3*u5
*u8*x3 - u3*u6*u7*x1 + u3*u6*u8*x2 + u3*u7*x1*x3 - u3*u8*x1*x2 + u4*u5*u7*x3 -
u4*u5*u8*u9 + u4*u6*u7*u9 - u4*u6*u7*x2 - u4*u7*u9*x3 + u4*u8*u9*x2 - u5*u7*x1*
x3 + u5*u8*u9*x3 + u6*u7*x1*x2 - u6*u8*u9*x2$
sol:=solve(polys,{x1,x2,x3});
sol := {{x1=(u1*u8 - u2*u7 + u2*u9 - u3*u8 + u4*u7 - u4*u9)/(u1 - u3),
x2=(u1*u9 - u3*u7 + u5*u7 - u5*u9)/(u1 - u3),
x3=(u1*u8 - u2*u7 + u2*u9 - u3*u8 + u6*u7 - u6*u9)/(u1 - u3)}}$
sub(sol,con);
0$
% The general theorem of Desargue.
A:=Point(0,0);
a := {0,0}$
B:=Point(0,1);
b := {0,1}$
C:=Point(u5,u6);
c := {u5,u6}$
R:=Point(u7,u8);
r := {u7,u8}$
S:=Point(u9,u1);
s := {u9,u1}$
T:=Point(u2,x1);
t := {u2,x1}$
con1:=collinear(intersection_point(pp_line(R,S),pp_line(A,B)),
intersection_point(pp_line(S,T),pp_line(B,C)),
intersection_point(pp_line(R,T),pp_line(A,C)));
con1 := (u5*( - u1**2*u2**2*u6*u7 + u1**2*u2*u5*u7*x1 + u1**2*u2*u6*u7**2 - u1**
2*u5*u7**2*x1 + u1*u2**2*u6*u7*u8 + u1*u2**2*u6*u7 + u1*u2**2*u6*u8*u9 - u1*u2**
2*u8*u9 - u1*u2*u5*u7*u8*x1 - u1*u2*u5*u7*x1 - u1*u2*u5*u8*u9*x1 + u1*u2*u5*u8*
u9 - u1*u2*u6*u7**2*x1 - u1*u2*u6*u7**2 - 2*u1*u2*u6*u7*u8*u9 + u1*u2*u6*u7*u9*
x1 - u1*u2*u6*u7*u9 + u1*u2*u7*u8*u9 + u1*u2*u7*u9*x1 + u1*u5*u7**2*x1**2 + u1*
u5*u7**2*x1 + 2*u1*u5*u7*u8*u9*x1 - u1*u5*u7*u8*u9 - u1*u5*u7*u9*x1**2 + u1*u6*
u7**2*u9 - u1*u7**2*u9*x1 - u2**2*u6*u7*u8 - u2**2*u6*u8**2*u9 + u2**2*u8**2*u9
+ u2*u5*u7*u8*x1 + u2*u5*u8**2*u9*x1 - u2*u5*u8**2*u9 + u2*u6*u7**2*x1 + u2*u6*
u7*u8*u9*x1 + 2*u2*u6*u7*u8*u9 - u2*u6*u7*u9*x1 + u2*u6*u8**2*u9**2 - u2*u6*u8*
u9**2*x1 - 2*u2*u7*u8*u9*x1 - u2*u8**2*u9**2 + u2*u8*u9**2*x1 - u5*u7**2*x1**2 -
u5*u7*u8*u9*x1**2 + u5*u7*u9*x1**2 - u5*u8**2*u9**2*x1 + u5*u8**2*u9**2 + u5*u8
*u9**2*x1**2 - u5*u8*u9**2*x1 - u6*u7**2*u9*x1 - u6*u7*u8*u9**2 + u6*u7*u9**2*x1
+ u7**2*u9*x1**2 + u7*u8*u9**2*x1 - u7*u9**2*x1**2))/(u1*u2*u5*u6*u7 - u1*u2*u5
*u6*u9 + u1*u5**2*u7*u8 - u1*u5**2*u7*x1 - u1*u5**2*u8*u9 + u1*u5**2*u9*x1 - u1*
u5*u6*u7**2 + u1*u5*u6*u7*u9 + u2**2*u6**2*u7 - u2**2*u6**2*u9 - u2**2*u6*u7 +
u2**2*u6*u9 + u2*u5*u6*u7*u8 - 2*u2*u5*u6*u7*x1 - u2*u5*u6*u8*u9 + 2*u2*u5*u6*u9
*x1 - u2*u5*u7*u8 + u2*u5*u7*x1 + u2*u5*u8*u9 - u2*u5*u9*x1 - u2*u6**2*u7**2 +
u2*u6**2*u9**2 + u2*u6*u7**2 - u2*u6*u9**2 - u5**2*u7*u8*x1 + u5**2*u7*x1**2 +
u5**2*u8*u9*x1 - u5**2*u9*x1**2 + u5*u6*u7**2*x1 - u5*u6*u7*u8*u9 + u5*u6*u8*u9
**2 - u5*u6*u9**2*x1 + u5*u7*u8*u9 - u5*u7*u9*x1 - u5*u8*u9**2 + u5*u9**2*x1 +
u6**2*u7**2*u9 - u6**2*u7*u9**2 - u6*u7**2*u9 + u6*u7*u9**2)$
con2:=concurrent(pp_line(A,R),pp_line(B,S),pp_line(C,T));
con2 := u1*u2*u6*u7 - u1*u5*u7*x1 - u2*u6*u7 - u2*u6*u8*u9 + u2*u8*u9 + u5*u7*x1
+ u5*u8*u9*x1 - u5*u8*u9 + u6*u7*u9 - u7*u9*x1$
sol:=solve(con2,x1);
sol := {x1=(u1*u2*u6*u7 - u2*u6*u7 - u2*u6*u8*u9 + u2*u8*u9 - u5*u8*u9 + u6*u7*
u9)/(u1*u5*u7 - u5*u7 - u5*u8*u9 + u7*u9)}$
sub(sol,con1);
0$
clear_ndg();
{}$
clear(A,B,C,R,S,T);
% #################
% Brocard points
% #################
A:=Point(0,0);
a := {0,0}$
B:=Point(1,0);
b := {1,0}$
C:=Point(u1,u2);
c := {u1,u2}$
c1:=Circle(1,x1,x2,x3);
c1 := {1,x1,x2,x3}$
c2:=Circle(1,x4,x5,x6);
c2 := {1,x4,x5,x6}$
c3:=Circle(1,x7,x8,x9);
c3 := {1,x7,x8,x9}$
polys:={
cl_tangent(c1,pp_line(A,C)),
point_on_circle(A,c1),
point_on_circle(B,c1),
cl_tangent(c2,pp_line(A,B)),
point_on_circle(B,c2),
point_on_circle(C,c2),
cl_tangent(c3,pp_line(B,C)),
point_on_circle(A,c3),
point_on_circle(C,c3)};
polys := {u1**2*x1**2 - 4*u1**2*x3 + 2*u1*u2*x1*x2 + u2**2*x2**2 - 4*u2**2*x3,
x3,
x1 + x3 + 1,
x4**2 - 4*x6,
x4 + x6 + 1,
u1**2 + u1*x4 + u2**2 + u2*x5 + x6,
u1**2*x7**2 - 4*u1**2*x9 + 2*u1*u2*x7*x8 + 4*u1*u2*x8 - 2*u1*x7**2 + 8*u1*x9 - 4
*u2**2*x7 + u2**2*x8**2 - 4*u2**2*x9 - 4*u2**2 - 2*u2*x7*x8 - 4*u2*x8 + x7**2 -
4*x9,
x9,
u1**2 + u1*x7 + u2**2 + u2*x8 + x9}$
vars:={x1,x2,x3,x4,x5,x6,x7,x8,x9};
vars := {x1,
x2,
x3,
x4,
x5,
x6,
x7,
x8,
x9}$
sol:=solve(polys,vars);
sol := {{x1=-1,
x2=u1/u2,
x3=0,
x4=-2,
x5=( - u1**2 + 2*u1 - u2**2 - 1)/u2,
x6=1,
x7=u1**2 - 2*u1 + u2**2,
x8=( - u1**3 + u1**2 - u1*u2**2 - u2**2)/u2,
x9=0}}$
P:=other_cc_point(B,sub(sol,c1),sub(sol,c2));
p := {(u1**3 - u1**2 + u1*u2**2 + u1 + u2**2)/(u1**4 - 2*u1**3 + 2*u1**2*u2**2 +
3*u1**2 - 2*u1*u2**2 - 2*u1 + u2**4 + 3*u2**2 + 1),
(u2*(u1**2 - 2*u1 + u2**2 + 1))/(u1**4 - 2*u1**3 + 2*u1**2*u2**2 + 3*u1**2 - 2*
u1*u2**2 - 2*u1 + u2**4 + 3*u2**2 + 1)}$
con:=point_on_circle(P,sub(sol,c3));
con := 0$
clear_ndg();
{}$
clear A,B,C,c1,c2,c3;
% ##################
% Simson's theorem
% ##################
% A constructive proof
M:=Point(0,0);
m := {0,0}$
A:=choose_pc(M,r,u1);
a := {(r*(u1**2 - 1))/(u1**2 + 1),(2*r*u1)/(u1**2 + 1)}$
B:=choose_pc(M,r,u2);
b := {(r*(u2**2 - 1))/(u2**2 + 1),(2*r*u2)/(u2**2 + 1)}$
C:=choose_pc(M,r,u3);
c := {(r*(u3**2 - 1))/(u3**2 + 1),(2*r*u3)/(u3**2 + 1)}$
P:=choose_pc(M,r,u4);
p := {(r*(u4**2 - 1))/(u4**2 + 1),(2*r*u4)/(u4**2 + 1)}$
X:=pedalpoint(P,pp_line(A,B))$
Y:=pedalpoint(P,pp_line(B,C))$
Z:=pedalpoint(P,pp_line(A,C))$
collinear(X,Y,Z);
0$
clear_ndg();
{}$
clear(M,A,B,C,P,X,Y,Z);
% Simson's theorem almost constructive
clear_ndg();
{}$
A:=Point(0,0);
a := {0,0}$
B:=Point(u1,u2);
b := {u1,u2}$
C:=Point(u3,u4);
c := {u3,u4}$
P:=Point(u5,x1);
p := {u5,x1}$
X:=pedalpoint(P,pp_line(A,B));
x := {(u1*(u1*u5 + u2*x1))/(u1**2 + u2**2),
(u2*(u1*u5 + u2*x1))/(u1**2 + u2**2)}$
Y:=pedalpoint(P,pp_line(B,C));
y := {(u1**2*u5 - u1*u2*u4 + u1*u2*x1 - 2*u1*u3*u5 + u1*u4**2 - u1*u4*x1 + u2**2
*u3 - u2*u3*u4 - u2*u3*x1 + u3**2*u5 + u3*u4*x1)/(u1**2 - 2*u1*u3 + u2**2 - 2*u2
*u4 + u3**2 + u4**2),
(u1**2*u4 - u1*u2*u3 + u1*u2*u5 - u1*u3*u4 - u1*u4*u5 + u2**2*x1 + u2*u3**2 - u2
*u3*u5 - 2*u2*u4*x1 + u3*u4*u5 + u4**2*x1)/(u1**2 - 2*u1*u3 + u2**2 - 2*u2*u4 +
u3**2 + u4**2)}$
Z:=pedalpoint(P,pp_line(A,C));
z := {(u3*(u3*u5 + u4*x1))/(u3**2 + u4**2),
(u4*(u3*u5 + u4*x1))/(u3**2 + u4**2)}$
poly:=p4_circle(A,B,C,P);
poly := u1**2*u3*x1 - u1**2*u4*u5 - u1*u3**2*x1 - u1*u4**2*x1 + u1*u4*u5**2 + u1
*u4*x1**2 + u2**2*u3*x1 - u2**2*u4*u5 + u2*u3**2*u5 - u2*u3*u5**2 - u2*u3*x1**2
+ u2*u4**2*u5$
con:=collinear(X,Y,Z);
con := ( - u1**4*u3*u4**2*x1 + u1**4*u4**3*u5 + 2*u1**3*u2*u3**2*u4*x1 - 2*u1**3
*u2*u3*u4**2*u5 + u1**3*u3**2*u4**2*x1 + u1**3*u4**4*x1 - u1**3*u4**3*u5**2 - u1
**3*u4**3*x1**2 - u1**2*u2**2*u3**3*x1 + u1**2*u2**2*u3**2*u4*u5 - u1**2*u2**2*
u3*u4**2*x1 + u1**2*u2**2*u4**3*u5 - 2*u1**2*u2*u3**3*u4*x1 - u1**2*u2*u3**2*u4
**2*u5 - 2*u1**2*u2*u3*u4**3*x1 + 3*u1**2*u2*u3*u4**2*u5**2 + 3*u1**2*u2*u3*u4**
2*x1**2 - u1**2*u2*u4**4*u5 + 2*u1*u2**3*u3**2*u4*x1 - 2*u1*u2**3*u3*u4**2*u5 +
u1*u2**2*u3**4*x1 + 2*u1*u2**2*u3**3*u4*u5 + u1*u2**2*u3**2*u4**2*x1 - 3*u1*u2**
2*u3**2*u4*u5**2 - 3*u1*u2**2*u3**2*u4*x1**2 + 2*u1*u2**2*u3*u4**3*u5 - u2**4*u3
**3*x1 + u2**4*u3**2*u4*u5 - u2**3*u3**4*u5 + u2**3*u3**3*u5**2 + u2**3*u3**3*x1
**2 - u2**3*u3**2*u4**2*u5)/(u1**4*u3**2 + u1**4*u4**2 - 2*u1**3*u3**3 - 2*u1**3
*u3*u4**2 + 2*u1**2*u2**2*u3**2 + 2*u1**2*u2**2*u4**2 - 2*u1**2*u2*u3**2*u4 - 2*
u1**2*u2*u4**3 + u1**2*u3**4 + 2*u1**2*u3**2*u4**2 + u1**2*u4**4 - 2*u1*u2**2*u3
**3 - 2*u1*u2**2*u3*u4**2 + u2**4*u3**2 + u2**4*u4**2 - 2*u2**3*u3**2*u4 - 2*u2
**3*u4**3 + u2**2*u3**4 + 2*u2**2*u3**2*u4**2 + u2**2*u4**4)$
remainder(num con,poly);
0$
print_ndg();
{u3**2 + u4**2,
u1**2 - 2*u1*u3 + u2**2 - 2*u2*u4 + u3**2 + u4**2,
u1**2 + u2**2}$
% Equational proof, first version:
M:=Point(0,0);
m := {0,0}$
A:=Point(0,1);
a := {0,1}$
B:=Point(u1,x1);
b := {u1,x1}$
C:=Point(u2,x2);
c := {u2,x2}$
P:=Point(u3,x3);
p := {u3,x3}$
X:=varpoint(A,B,x4);
x := {u1*( - x4 + 1), - x1*x4 + x1 + x4}$
Y:=varpoint(B,C,x5);
y := {u1*x5 - u2*x5 + u2,x1*x5 - x2*x5 + x2}$
Z:=varpoint(A,C,x6);
z := {u2*( - x6 + 1), - x2*x6 + x2 + x6}$
polys:={sqrdist(M,B)-1, sqrdist(M,C)-1, sqrdist(M,P)-1,
orthogonal(pp_line(A,B),pp_line(P,X)),
orthogonal(pp_line(A,C),pp_line(P,Z)),
orthogonal(pp_line(B,C),pp_line(P,Y))};
polys := {u1**2 + x1**2 - 1,
u2**2 + x2**2 - 1,
u3**2 + x3**2 - 1,
- u1**2*x4 + u1**2 - u1*u3 - x1**2*x4 + x1**2 - x1*x3 + 2*x1*x4 - x1 + x3 - x4,
- u2**2*x6 + u2**2 - u2*u3 - x2**2*x6 + x2**2 - x2*x3 + 2*x2*x6 - x2 + x3 - x6,
- u1**2*x5 + 2*u1*u2*x5 - u1*u2 + u1*u3 - u2**2*x5 + u2**2 - u2*u3 - x1**2*x5 +
2*x1*x2*x5 - x1*x2 + x1*x3 - x2**2*x5 + x2**2 - x2*x3}$
con:=collinear(X,Y,Z);
con := u1*x2*x4*x5 - u1*x2*x4*x6 - u1*x2*x5*x6 + u1*x2*x6 - u1*x4*x5 + u1*x4*x6
+ u1*x5*x6 - u1*x6 - u2*x1*x4*x5 + u2*x1*x4*x6 + u2*x1*x5*x6 - u2*x1*x6 + u2*x4*
x5 - u2*x4*x6 - u2*x5*x6 + u2*x6$
vars:={x4,x5,x6,x1,x2,x3};
vars := {x4,
x5,
x6,
x1,
x2,
x3}$
setring(vars,{},lex);
{{x4,x5,x6,x1,x2,x3},{},lex,{1,1,1,1,1,1}}$
setideal(polys,polys);
{x1**2 + (u1**2 - 1),
x2**2 + (u2**2 - 1),
x3**2 + (u3**2 - 1),
- x4*x1**2 + 2*x4*x1 - (u1**2 + 1)*x4 + x1**2 - x1*x3 - x1 + x3 + (u1**2 - u1*
u3),
- x6*x2**2 + 2*x6*x2 - (u2**2 + 1)*x6 + x2**2 - x2*x3 - x2 + x3 + (u2**2 - u2*
u3),
- x5*x1**2 + 2*x5*x1*x2 - x5*x2**2 - (u1**2 - 2*u1*u2 + u2**2)*x5 - x1*x2 + x1*
x3 + x2**2 - x2*x3 - (u1*u2 - u1*u3 - u2**2 + u2*u3)}$
con mod gbasis polys;
0$
% Second version:
A:=Point(0,0);
a := {0,0}$
B:=Point(1,0);
b := {1,0}$
C:=Point(u1,u2);
c := {u1,u2}$
P:=Point(u3,x1);
p := {u3,x1}$
X:=Point(x2,0);
x := {x2,0}$
% => on the line AB
Y:=varpoint(B,C,x3);
y := { - u1*x3 + u1 + x3,u2*( - x3 + 1)}$
Z:=varpoint(A,C,x4);
z := {u1*( - x4 + 1),u2*( - x4 + 1)}$
polys:={orthogonal(pp_line(A,C),pp_line(P,Z)),
orthogonal(pp_line(B,C),pp_line(P,Y)),
orthogonal(pp_line(A,B),pp_line(P,X)),
p4_circle(A,B,C,P)};
polys := { - u1**2*x4 + u1**2 - u1*u3 - u2**2*x4 + u2**2 - u2*x1,
- u1**2*x3 + u1**2 - u1*u3 + 2*u1*x3 - u1 - u2**2*x3 + u2**2 - u2*x1 + u3 - x3,
- u3 + x2,
- u1**2*x1 + u1*x1 - u2**2*x1 + u2*u3**2 - u2*u3 + u2*x1**2}$
con:=collinear(X,Y,Z);
con := u2*( - x2*x3 + x2*x4 - x3*x4 + x3)$
vars:={x2,x3,x4,x1};
vars := {x2,x3,x4,x1}$
setring(vars,{},lex);
{{x2,x3,x4,x1},{},lex,{1,1,1,1}}$
con mod interreduce polys;
0$
% The inverse theorem
polys:={orthogonal(pp_line(A,C),pp_line(P,Z)),
orthogonal(pp_line(B,C),pp_line(P,Y)),
orthogonal(pp_line(A,B),pp_line(P,X)),
collinear(X,Y,Z)};
polys := { - u1**2*x4 + u1**2 - u1*u3 - u2**2*x4 + u2**2 - u2*x1,
- u1**2*x3 + u1**2 - u1*u3 + 2*u1*x3 - u1 - u2**2*x3 + u2**2 - u2*x1 + u3 - x3,
- u3 + x2,
u2*( - x2*x3 + x2*x4 - x3*x4 + x3)}$
con:=p4_circle(A,B,C,P);
con := - u1**2*x1 + u1*x1 - u2**2*x1 + u2*u3**2 - u2*u3 + u2*x1**2$
con mod interreduce polys;
0$
clear_ndg();
{}$
clear(M,A,B,C,P,Y,Z);
% ########################
% The butterfly theorem
% ########################
% An equational proof with groebner factorizer and constraints.
P:=Point(0,0);
p := {0,0}$
O:=Point(u1,0);
o := {u1,0}$
A:=Point(u2,u3);
a := {u2,u3}$
B:=Point(u4,x1);
b := {u4,x1}$
C:=Point(x2,x3);
c := {x2,x3}$
D:=Point(x4,x5);
d := {x4,x5}$
F:=Point(0,x6);
f := {0,x6}$
G:=Point(0,x7);
g := {0,x7}$
polys:={sqrdist(O,B)-sqrdist(O,A),
sqrdist(O,C)-sqrdist(O,A),
sqrdist(O,D)-sqrdist(O,A),
point_on_line(P,pp_line(A,C)),
point_on_line(P,pp_line(B,D)),
point_on_line(F,pp_line(A,D)),
point_on_line(G,pp_line(B,C))
};
polys := {2*u1*u2 - 2*u1*u4 - u2**2 - u3**2 + u4**2 + x1**2,
2*u1*u2 - 2*u1*x2 - u2**2 - u3**2 + x2**2 + x3**2,
2*u1*u2 - 2*u1*x4 - u2**2 - u3**2 + x4**2 + x5**2,
- u2*x3 + u3*x2,
- u4*x5 + x1*x4,
- u2*x5 + u2*x6 + u3*x4 - x4*x6,
- u4*x3 + u4*x7 + x1*x2 - x2*x7}$
con:=num sqrdist(P,midpoint(F,G));
con := x6**2 + 2*x6*x7 + x7**2$
vars:={x6,x7,x3,x5,x1,x2,x4};
vars := {x6,
x7,
x3,
x5,
x1,
x2,
x4}$
setring(vars,{},lex);
{{x6,x7,x3,x5,x1,x2,x4},{},lex,{1,1,1,1,1,1,1}}$
sol:=groebfactor(polys,{sqrdist(A,C),sqrdist(B,D)});
sol := {{x1**2 + (2*u1*u2 - 2*u1*u4 - u2**2 - u3**2 + u4**2),
(u2**2 + u3**2)*x3 - (2*u1*u2*u3 - u2**2*u3 - u3**3),
(2*u1*u2 - 2*u1*u4 - u2**2 - u3**2)*x5 + (2*u1*u2 - u2**2 - u3**2)*x1,
(2*u1*u2 - 2*u1*u4 - u2**2 - u3**2)*x4 + (2*u1*u2*u4 - u2**2*u4 - u3**2*u4),
(u2**2 + u3**2)*x2 - (2*u1*u2**2 - u2**3 - u2*u3**2),
(2*u1*u2**2 - u2**3 - u2**2*u4 - u2*u3**2 - u3**2*u4)*x7 - (2*u1*u2**2 - u2**3 -
u2*u3**2)*x1 + (2*u1*u2*u3*u4 - u2**2*u3*u4 - u3**3*u4),
(2*u1*u2**2 - u2**3 - u2**2*u4 - u2*u3**2 - u3**2*u4)*x6 + (2*u1*u2**2 - u2**3 -
u2*u3**2)*x1 - (2*u1*u2*u3*u4 - u2**2*u3*u4 - u3**3*u4)}}$
for each u in sol collect con mod u;
{0}$
% A constructive proof
on gcd;
O:=Point(0,0);
o := {0,0}$
A:=Point(1,0);
a := {1,0}$
B:=choose_pc(O,1,u1);
b := {(u1**2 - 1)/(u1**2 + 1),(2*u1)/(u1**2 + 1)}$
C:=choose_pc(O,1,u2);
c := {(u2**2 - 1)/(u2**2 + 1),(2*u2)/(u2**2 + 1)}$
D:=choose_pc(O,1,u3);
d := {(u3**2 - 1)/(u3**2 + 1),(2*u3)/(u3**2 + 1)}$
P:=intersection_point(pp_line(A,C),pp_line(B,D));
p := {(u1*u2 - u1*u3 + u2*u3 - 1)/(u1*u2 - u1*u3 + u2*u3 + 1),
(2*u2)/(u1*u2 - u1*u3 + u2*u3 + 1)}$
h:=lot(P,pp_line(O,P));
h := {( - u1*u2 + u1*u3 - u2*u3 + 1)/(u1*u2 - u1*u3 + u2*u3 + 1),
( - 2*u2)/(u1*u2 - u1*u3 + u2*u3 + 1),
(u1**2*u2**2 - 2*u1**2*u2*u3 + u1**2*u3**2 + 2*u1*u2**2*u3 - 2*u1*u2*u3**2 - 2*
u1*u2 + 2*u1*u3 + u2**2*u3**2 + 4*u2**2 - 2*u2*u3 + 1)/(u1**2*u2**2 - 2*u1**2*u2
*u3 + u1**2*u3**2 + 2*u1*u2**2*u3 - 2*u1*u2*u3**2 + 2*u1*u2 - 2*u1*u3 + u2**2*u3
**2 + 2*u2*u3 + 1)}$
F:=intersection_point(h,pp_line(A,D));
f := {(u1**2*u2**2 - 2*u1**2*u2*u3 + u1**2*u3**2 - 2*u1*u2 + 2*u1*u3 - u2**2*u3
**2 + 4*u2**2 - 4*u2*u3 + 1)/(u1**2*u2**2 - 2*u1**2*u2*u3 + u1**2*u3**2 - u2**2*
u3**2 - 2*u2*u3 - 1),
(2*u3*(u1*u2 - u1*u3 - 2*u2**2 + u2*u3 - 1))/(u1**2*u2**2 - 2*u1**2*u2*u3 + u1**
2*u3**2 - u2**2*u3**2 - 2*u2*u3 - 1)}$
G:=intersection_point(h,pp_line(B,C));
g := {(u1**2*u2**2 - 2*u1**2*u2*u3 + u1**2*u3**2 - 2*u1*u2 + 2*u1*u3 - u2**2*u3
**2 - 4*u2**2 + 4*u2*u3 + 1)/(u1**2*u2**2 - 2*u1**2*u2*u3 + u1**2*u3**2 - u2**2*
u3**2 - 2*u2*u3 - 1),
(2*(2*u1*u2**2 - 3*u1*u2*u3 + u1*u3**2 - u2*u3**2 - 2*u2 + u3))/(u1**2*u2**2 - 2
*u1**2*u2*u3 + u1**2*u3**2 - u2**2*u3**2 - 2*u2*u3 - 1)}$
con:=sqrdist(P,midpoint(F,G));
con := 0$
off gcd;
clear_ndg();
{}$
clear(O,A,B,C,D,P,h,F,G);
% ################################
% Tangency of Feuerbach's circle
% ################################
A:=Point(0,0);
a := {0,0}$
B:=Point(2,0);
b := {2,0}$
C:=Point(u1,u2);
c := {u1,u2}$
M:=intersection_point(mp(A,B),mp(B,C));
m := {1,(u1**2 - 2*u1 + u2**2)/(2*u2)}$
H:=intersection_point(altitude(A,B,C),altitude(B,C,A));
h := {u1,(u1*( - u1 + 2))/u2}$
N:=midpoint(M,H);
n := {(u1 + 1)/2,( - u1**2 + 2*u1 + u2**2)/(4*u2)}$
c1:=c1_circle(N,sqrdist(N,midpoint(A,B)));
c1 := {1, - (u1 + 1),(u1**2 - 2*u1 - u2**2)/(2*u2),u1}$
% Feuerbach's circle
P:=Point(x1,x2);
p := {x1,x2}$
% => x2 is the radius of the inscribed circle.
polys:={point_on_bisector(P,A,B,C), point_on_bisector(P,B,C,A)};
polys := {2*( - 2*u1*x1*x2 + 4*u1*x2 + u2*x1**2 - 4*u2*x1 - u2*x2**2 + 4*u2 + 4*
x1*x2 - 8*x2),
2*( - u1**3*x2 + u1**2*u2*x1 - u1**2*u2 + u1**2*x1*x2 + 2*u1**2*x2 - u1*u2**2*x2
- u1*u2*x1**2 + u1*u2*x2**2 - 2*u1*x1*x2 + u2**3*x1 - u2**3 - u2**2*x1*x2 + 2*
u2**2*x2 + u2*x1**2 - u2*x2**2)}$
con:=cc_tangent(c1_circle(P,x2^2),c1);
con := (4*( - u1**3*x1*x2 + u1**3*x2 + u1**2*u2*x1**2 - 2*u1**2*u2*x1 - u1**2*u2
*x2**2 + u1**2*u2 + u1**2*x1**2*x2 + u1**2*x1*x2 - 2*u1**2*x2 + u1*u2**2*x1*x2 -
u1*u2**2*x2 - 2*u1*u2*x1**3 + 4*u1*u2*x1**2 - 2*u1*u2*x1 + 2*u1*u2*x2**2 - 2*u1
*x1**2*x2 + 2*u1*x1*x2 - u2**2*x1**2*x2 + u2**2*x1*x2 + u2*x1**4 - 2*u2*x1**3 +
u2*x1**2 - u2*x2**2))/u2$
vars:={x1,x2};
vars := {x1,x2}$
setring(vars,{},lex);
{{x1,x2},{},lex,{1,1}}$
setideal(polys,polys);
{(2*u2)*x1**2 - (4*u1 - 8)*x1*x2 - (8*u2)*x1 - (2*u2)*x2**2 + (8*u1 - 16)*x2 + 8
*u2,
- (2*u1*u2 - 2*u2)*x1**2 + (2*u1**2 - 4*u1 - 2*u2**2)*x1*x2 + (2*u1**2*u2 + 2*
u2**3)*x1 + (2*u1*u2 - 2*u2)*x2**2 - (2*u1**3 - 4*u1**2 + 2*u1*u2**2 - 4*u2**2)*
x2 - (2*u1**2*u2 + 2*u2**3)}$
num con mod gbasis polys;
0$
% Now let P be the incenter of the triangle ABH
polys1:={point_on_bisector(P,A,B,H), point_on_bisector(P,B,H,A)};
polys1 := {(2*( - u1**2*x1**2 + 4*u1**2*x1 + u1**2*x2**2 - 4*u1**2 - 2*u1*u2*x1*
x2 + 4*u1*u2*x2 + 2*u1*x1**2 - 8*u1*x1 - 2*u1*x2**2 + 8*u1 + 4*u2*x1*x2 - 8*u2*
x2))/u2,
(2*u1*( - u1**5*x1 + u1**5 - u1**4*u2*x2 + 6*u1**4*x1 - 6*u1**4 - u1**3*u2**2*x1
+ u1**3*u2**2 - u1**3*u2*x1*x2 + 6*u1**3*u2*x2 - 12*u1**3*x1 + 12*u1**3 - u1**2
*u2**3*x2 + u1**2*u2**2*x1**2 + 2*u1**2*u2**2*x1 - u1**2*u2**2*x2**2 - 2*u1**2*
u2**2 + 4*u1**2*u2*x1*x2 - 12*u1**2*u2*x2 + 8*u1**2*x1 - 8*u1**2 + u1*u2**3*x1*
x2 + 2*u1*u2**3*x2 - 3*u1*u2**2*x1**2 + 3*u1*u2**2*x2**2 - 4*u1*u2*x1*x2 + 8*u1*
u2*x2 - 2*u2**3*x1*x2 + 2*u2**2*x1**2 - 2*u2**2*x2**2))/u2**3}$
con1:=cc_tangent(c1_circle(P,x2^2),c1);
con1 := (4*( - u1**3*x1*x2 + u1**3*x2 + u1**2*u2*x1**2 - 2*u1**2*u2*x1 - u1**2*
u2*x2**2 + u1**2*u2 + u1**2*x1**2*x2 + u1**2*x1*x2 - 2*u1**2*x2 + u1*u2**2*x1*x2
- u1*u2**2*x2 - 2*u1*u2*x1**3 + 4*u1*u2*x1**2 - 2*u1*u2*x1 + 2*u1*u2*x2**2 - 2*
u1*x1**2*x2 + 2*u1*x1*x2 - u2**2*x1**2*x2 + u2**2*x1*x2 + u2*x1**4 - 2*u2*x1**3
+ u2*x1**2 - u2*x2**2))/u2$
setideal(polys1,polys1);
{ - (2*u1**2 - 4*u1)*x1**2 - (4*u1*u2 - 8*u2)*x1*x2 + (8*u1**2 - 16*u1)*x1 + (2*
u1**2 - 4*u1)*x2**2 + (8*u1*u2 - 16*u2)*x2 - (8*u1**2 - 16*u1),
(2*u1**3*u2**2 - 6*u1**2*u2**2 + 4*u1*u2**2)*x1**2 - (2*u1**4*u2 - 8*u1**3*u2 -
2*u1**2*u2**3 + 8*u1**2*u2 + 4*u1*u2**3)*x1*x2 - (2*u1**6 - 12*u1**5 + 2*u1**4*
u2**2 + 24*u1**4 - 4*u1**3*u2**2 - 16*u1**3)*x1 - (2*u1**3*u2**2 - 6*u1**2*u2**2
+ 4*u1*u2**2)*x2**2 - (2*u1**5*u2 - 12*u1**4*u2 + 2*u1**3*u2**3 + 24*u1**3*u2 -
4*u1**2*u2**3 - 16*u1**2*u2)*x2 + (2*u1**6 - 12*u1**5 + 2*u1**4*u2**2 + 24*u1**
4 - 4*u1**3*u2**2 - 16*u1**3)}$
num con1 mod gbasis polys1;
0$
clear_ndg();
{}$
clear A,B,C,P,M,N,H,c1;
% #############################
% Solutions to the exercises
% #############################
% 1)
A:=Point(0,0);
a := {0,0}$
B:=Point(1,0);
b := {1,0}$
C:=Point(1,1);
c := {1,1}$
D:=Point(0,1);
d := {0,1}$
P:=Point(x1,x2);
p := {x1,x2}$
Q:=Point(x3,1);
q := {x3,1}$
polys:={point_on_line(P,par(C,pp_line(B,D))),
sqrdist(B,D)-sqrdist(B,P),
point_on_line(Q,pp_line(B,P))};
polys := {x1 + x2 - 2,
- x1**2 + 2*x1 - x2**2 + 1,
- x1 + x2*x3 - x2 + 1}$
con:=sqrdist(D,P)-sqrdist(D,Q);
con := x1**2 + x2**2 - 2*x2 - x3**2 + 1$
setring({x1,x2,x3},{},lex);
{{x1,x2,x3},{},lex,{1,1,1}}$
setideal(polys,polys);
{x1 + x2 - 2,
- x1**2 + 2*x1 - x2**2 + 1,
- x1 + x2*x3 - x2 + 1}$
con mod gbasis polys;
0$
clear_ndg();
{}$
clear(A,B,C,D,P,Q);
% 2)
A:=Point(u1,0);
a := {u1,0}$
B:=Point(u2,0);
b := {u2,0}$
C:=Point(0,u3);
c := {0,u3}$
Q:=Point(0,0);
q := {0,0}$
% the pedal point on AB
R:=pedalpoint(B,pp_line(A,C));
r := {(u1*(u1*u2 + u3**2))/(u1**2 + u3**2),
(u1*u3*(u1 - u2))/(u1**2 + u3**2)}$
P:=pedalpoint(A,pp_line(B,C));
p := {(u2*(u1*u2 + u3**2))/(u2**2 + u3**2),
(u2*u3*( - u1 + u2))/(u2**2 + u3**2)}$
con1:=point_on_bisector(C,P,Q,R);
con1 := 0$
con2:=angle_sum(p3_angle(P,Q,C),p3_angle(R,Q,C));
con2 := 0$
clear_ndg();
{}$
clear(A,B,C,P,Q,R);
% 3)
A:=Point(u1,0);
a := {u1,0}$
B:=Point(u2,0);
b := {u2,0}$
C:=Point(0,u3);
c := {0,u3}$
P:=pedalpoint(A,pp_line(B,C));
p := {(u2*(u1*u2 + u3**2))/(u2**2 + u3**2),
(u2*u3*( - u1 + u2))/(u2**2 + u3**2)}$
Q:=pedalpoint(B,pp_line(A,C));
q := {(u1*(u1*u2 + u3**2))/(u1**2 + u3**2),
(u1*u3*(u1 - u2))/(u1**2 + u3**2)}$
R:=pedalpoint(C,pp_line(A,B));
r := {0,0}$
P1:=pedalpoint(P,pp_line(A,B));
p1 := {(u2*(u1*u2 + u3**2))/(u2**2 + u3**2),0}$
P2:=pedalpoint(P,pp_line(A,C));
p2 := {(u1*(u1**2*u2**2 + 2*u1*u2*u3**2 + u3**4))/(u1**2*u2**2 + u1**2*u3**2 +
u2**2*u3**2 + u3**4),
(u3**3*(u1**2 - 2*u1*u2 + u2**2))/(u1**2*u2**2 + u1**2*u3**2 + u2**2*u3**2 + u3
**4)}$
Q1:=pedalpoint(Q,pp_line(A,B));
q1 := {(u1*(u1*u2 + u3**2))/(u1**2 + u3**2),0}$
Q2:=pedalpoint(Q,pp_line(B,C));
q2 := {(u2*(u1**2*u2**2 + 2*u1*u2*u3**2 + u3**4))/(u1**2*u2**2 + u1**2*u3**2 +
u2**2*u3**2 + u3**4),
(u3**3*(u1**2 - 2*u1*u2 + u2**2))/(u1**2*u2**2 + u1**2*u3**2 + u2**2*u3**2 + u3
**4)}$
R1:=pedalpoint(R,pp_line(A,C));
r1 := {(u1*u3**2)/(u1**2 + u3**2),(u1**2*u3)/(u1**2 + u3**2)}$
R2:=pedalpoint(R,pp_line(B,C));
r2 := {(u2*u3**2)/(u2**2 + u3**2),(u2**2*u3)/(u2**2 + u3**2)}$
con:=for each X in {Q2,R1,R2} collect p4_circle(P1,P2,Q1,X);
con := {0,0,0}$
clear_ndg();
{}$
clear(O,A,B,C,P,Q,R,P1,P2,Q1,Q2,R1,R2);
% 4)
A:=Point(u1,0);
a := {u1,0}$
B:=Point(u2,0);
b := {u2,0}$
C:=Point(0,u3);
c := {0,u3}$
% => Pedalpoint from C is (0,0)
M:=intersection_point(mp(A,B),mp(B,C));
m := {(u1 + u2)/2,(u1*u2 + u3**2)/(2*u3)}$
% Prove (2*h_c*R = a*b)^2
con:=4*u3^2*sqrdist(M,A)-sqrdist(C,B)*sqrdist(A,C);
con := 0$
clear_ndg();
{}$
clear(A,B,C,M);
% 5. A solution of constructive type.
on gcd;
O:=Point(0,u1);
o := {0,u1}$
A:=Point(0,0);
a := {0,0}$
% hence k has radius u1.
B:=Point(u2,0);
b := {u2,0}$
M:=midpoint(A,B);
m := {u2/2,0}$
D:=choose_pc(O,u1,u3);
d := {(u1*(u3**2 - 1))/(u3**2 + 1),(u1*(u3**2 + 2*u3 + 1))/(u3**2 + 1)}$
k:=c1_circle(O,u1^2);
k := {1,0, - 2*u1,0}$
C:=other_cl_point(D,k,pp_line(M,D));
c := {(u1*u2*(4*u1*u3**2 + 8*u1*u3 + 4*u1 - u2*u3**2 + u2))/(8*u1**2*u3**2 + 16*
u1**2*u3 + 8*u1**2 - 4*u1*u2*u3**2 + 4*u1*u2 + u2**2*u3**2 + u2**2),
(u1*u2**2*(u3**2 + 2*u3 + 1))/(8*u1**2*u3**2 + 16*u1**2*u3 + 8*u1**2 - 4*u1*u2*
u3**2 + 4*u1*u2 + u2**2*u3**2 + u2**2)}$
Eh:=other_cl_point(D,k,pp_line(B,D));
eh := {(u1*u2*(2*u1*u3**2 + 4*u1*u3 + 2*u1 - u2*u3**2 + u2))/(2*u1**2*u3**2 + 4*
u1**2*u3 + 2*u1**2 - 2*u1*u2*u3**2 + 2*u1*u2 + u2**2*u3**2 + u2**2),
(u1*u2**2*(u3**2 + 2*u3 + 1))/(2*u1**2*u3**2 + 4*u1**2*u3 + 2*u1**2 - 2*u1*u2*u3
**2 + 2*u1*u2 + u2**2*u3**2 + u2**2)}$
F:=other_cl_point(C,k,pp_line(B,C));
f := {(u1*u2*( - 2*u1*u3**2 - 4*u1*u3 - 2*u1 + u2*u3**2 - u2))/(2*u1**2*u3**2 +
4*u1**2*u3 + 2*u1**2 - 2*u1*u2*u3**2 + 2*u1*u2 + u2**2*u3**2 + u2**2),
(u1*u2**2*(u3**2 + 2*u3 + 1))/(2*u1**2*u3**2 + 4*u1**2*u3 + 2*u1**2 - 2*u1*u2*u3
**2 + 2*u1*u2 + u2**2*u3**2 + u2**2)}$
con:=parallel(pp_line(A,B),pp_line(Eh,F));
con := 0$
off gcd;
clear_ndg();
{}$
clear(O,A,B,C,D,Eh,F,M,k);
% 6)
Z:=Point(0,0);
z := {0,0}$
X:=Point(0,1);
x := {0,1}$
Y:=Point(0,-1);
y := {0,-1}$
B:=Point(u1,0);
b := {u1,0}$
C:=Point(u2,0);
c := {u2,0}$
P:=Point(0,u3);
p := {0,u3}$
M:=Point(x1,x2);
m := {x1,x2}$
N:=Point(x3,x4);
n := {x3,x4}$
A:=Point(x5,0);
a := {x5,0}$
D:=Point(x6,0);
d := {x6,0}$
polys:={p4_circle(X,Y,B,N), p4_circle(X,Y,C,M),
p4_circle(X,Y,B,D), p4_circle(X,Y,C,A),
collinear(B,P,N), collinear(C,P,M)};
polys := {2*( - u1**2*x3 + u1*x3**2 + u1*x4**2 - u1 + x3),
2*( - u2**2*x1 + u2*x1**2 + u2*x2**2 - u2 + x1),
2*( - u1**2*x6 + u1*x6**2 - u1 + x6),
2*( - u2**2*x5 + u2*x5**2 - u2 + x5),
u1*u3 - u1*x4 - u3*x3,
u2*u3 - u2*x2 - u3*x1}$
con:=concurrent(pp_line(A,M),pp_line(D,N),pp_line(X,Y));
con := 2*( - x1*x4*x6 + x2*x3*x5 - x2*x5*x6 + x4*x5*x6)$
vars:={x1,x2,x3,x4,x5,x6};
vars := {x1,
x2,
x3,
x4,
x5,
x6}$
setring(vars,{},lex);
{{x1,x2,x3,x4,x5,x6},{},lex,{1,1,1,1,1,1}}$
res:=groebfactor(polys,{x5-u2,x1-u2,x6-u1,x3-u1});
res := {{u1*x6 + 1,
(u2**2 + u3**2)*x2 - (u2**2*u3 + u3),
(u2**2 + u3**2)*x1 - (u2*u3**2 - u2),
(u1**2 + u3**2)*x4 - (u1**2*u3 + u3),
(u1**2 + u3**2)*x3 - (u1*u3**2 - u1),
u2*x5 + 1}}$
% constraints A\neq C, M\neq C, D\neq B, N\neq B
for each u in res collect con mod u;
{0}$
clear_ndg();
{}$
clear(Z,X,Y,B,C,P,M,N,A,D);
% 7)
M:=Point(0,0);
m := {0,0}$
A:=Point(0,u1);
a := {0,u1}$
B:=Point(-1,0);
b := {-1,0}$
C:=Point(1,0);
c := {1,0}$
Eh:=varpoint(A,B,x1);
eh := {x1 - 1,u1*x1}$
F:=varpoint(A,C,x2);
f := { - x2 + 1,u1*x2}$
O:=intersection_point(pp_line(A,M),lot(B,pp_line(A,B)));
o := {0,( - 1)/u1}$
Q:=intersection_point(pp_line(Eh,F),pp_line(B,C));
q := {( - 2*x1*x2 + x1 + x2)/(x1 - x2),0}$
con1:=num orthogonal(pp_line(O,Q),pp_line(Eh,Q));
con1 := 2*x1*(x1**2*x2 - x1**2 + x1*x2**2 - 2*x1*x2 + x1 - x2**2 + x2)$
con2:=num sqrdist(Q,midpoint(Eh,F));
con2 := u1**2*x1**4 - 2*u1**2*x1**2*x2**2 + u1**2*x2**4 + x1**4 + 4*x1**3*x2 - 4
*x1**3 + 6*x1**2*x2**2 - 12*x1**2*x2 + 4*x1**2 + 4*x1*x2**3 - 12*x1*x2**2 + 8*x1
*x2 + x2**4 - 4*x2**3 + 4*x2**2$
vars:={x1,x2};
vars := {x1,x2}$
setring(vars,{},lex);
{{x1,x2},{},lex,{1,1}}$
p1:=groebfactor({con1},{x1-1,x2-1,x1,x2});
p1 := {{x1 + x2}}$
p2:=groebfactor({con2},{x1-1,x2-1,x1,x2});
p2 := {{x1 + x2},
{(u1**2 + 1)*x1**2 - (2*u1**2 - 2)*x1*x2 - 4*x1 + (u1**2 + 1)*x2**2 - 4*x2 + 4}}
$
% constraint A,C\neq Eh, B,C\neq F
for each u in p1 collect con2 mod u;
{0}$
for each u in p2 collect con1 mod u;
{0,
(2*(5*u1**4*x1*x2**3 - 8*u1**4*x1*x2**2 + 3*u1**4*x1*x2 - 3*u1**4*x2**4 + 4*u1**
4*x2**3 - u1**4*x2**2 - 10*u1**2*x1*x2**3 + 32*u1**2*x1*x2**2 - 30*u1**2*x1*x2 +
8*u1**2*x1 - 2*u1**2*x2**4 + 12*u1**2*x2**3 - 26*u1**2*x2**2 + 20*u1**2*x2 - 4*
u1**2 + x1*x2**3 - 8*x1*x2**2 + 15*x1*x2 - 8*x1 + x2**4 - 8*x2**3 + 23*x2**2 -
28*x2 + 12))/(u1**4 + 2*u1**2 + 1)}$
% Note that the second component of p2 has no relevant *real* roots,
% since it factors as u1^2 * (x1 - x2)^2 + (x1 + x2 -2)^2 :
u1^2 * (x1 - x2)^2 + (x1 + x2 -2)^2 mod second p2;
0$
clear_ndg();
{}$
clear(M,A,B,C,O,Eh,F,Q);
% 8)
on gcd;
A:=Point(u1,0);
a := {u1,0}$
B:=Point(u2,0);
b := {u2,0}$
l1:=pp_line(A,B);
l1 := {0,u1 - u2,0}$
M:=Point(0,u3);
m := {0,u3}$
% the incenter, hence u3 = incircle radius
C:=intersection_point(symline(l1,pp_line(A,M)),
symline(l1,pp_line(B,M)));
c := {(u3**2*(u1 + u2))/(u1*u2 + u3**2),
(2*u1*u2*u3)/(u1*u2 + u3**2)}$
N:=intersection_point(mp(A,B),mp(B,C));
n := {(u1 + u2)/2,
(u1**2*u2**2 - u1**2*u3**2 + 4*u1*u2*u3**2 - u2**2*u3**2 + u3**4)/(4*u3*(u1*u2 +
u3**2))}$
% the outcenter
sqr_rad:=sqrdist(A,N);
sqr_rad := (u1**4*u2**4 + 2*u1**4*u2**2*u3**2 + u1**4*u3**4 + 2*u1**2*u2**4*u3**
2 + 4*u1**2*u2**2*u3**4 + 2*u1**2*u3**6 + u2**4*u3**4 + 2*u2**2*u3**6 + u3**8)/(
16*u3**2*(u1**2*u2**2 + 2*u1*u2*u3**2 + u3**4))$
% the outcircle sqradius.
(sqr_rad-sqrdist(M,N))^2-4*u3^2*sqr_rad;
0$
off gcd;
clear_ndg();
{}$
clear A,B,C,M,N,l1,sqr_rad;
% 9)
on gcd;
A:=Point(0,0);
a := {0,0}$
B:=Point(1,0);
b := {1,0}$
M:=Point(u1,0);
m := {u1,0}$
C:=Point(u1,u1);
c := {u1,u1}$
F:=Point(u1,1-u1);
f := {u1, - u1 + 1}$
c1:=red_hom_coords p3_circle(A,M,C);
c1 := {1, - u1, - u1,0}$
c2:=red_hom_coords p3_circle(B,M,F);
c2 := {-1,u1 + 1, - u1 + 1, - u1}$
N:=other_cc_point(M,c1,c2);
n := {u1**2/(2*u1**2 - 2*u1 + 1),(u1*( - u1 + 1))/(2*u1**2 - 2*u1 + 1)}$
point_on_line(N,pp_line(A,F));
0$
point_on_line(N,pp_line(B,C));
0$
l1:=red_hom_coords pp_line(M,N);
l1 := {-1,2*u1 - 1,u1}$
l2:=sub(u1=u2,l1);
l2 := {-1,2*u2 - 1,u2}$
intersection_point(l1,l2);
{1/2,( - 1)/2}$
% = (1/2,-1/2)
off gcd;
clear_ndg();
{}$
clear A,B,C,F,M,N,c1,c2,l1,l2;
% ####################
% Some more examples
% ####################
% Origin: D. Wang at
% http://cosmos.imag.fr/ATINF/Dongming.Wang/geother.html
% --------------------------
% Given triangle ABC, H orthocenter, O circumcenter, A1 circumcenter
% of BHC, B1 circumcenter of AHC.
%
% Claim: OH, AA1, BB1 are concurrent.
% --------------------------
A:=Point(u1,0);
a := {u1,0}$
B:=Point(u2,0);
b := {u2,0}$
C:=Point(0,u3);
c := {0,u3}$
H:=intersection_point(altitude(C,A,B),altitude(A,B,C));
h := {0,( - u1*u2)/u3}$
O:=circle_center(p3_circle(A,B,C));
o := {(u1 + u2)/2,(u1*u2 + u3**2)/(2*u3)}$
A1:=circle_center(p3_circle(H,B,C));
a1 := {( - u1 + u2)/2,( - u1*u2 + u3**2)/(2*u3)}$
B1:=circle_center(p3_circle(H,A,C));
b1 := {(u1 - u2)/2,( - u1*u2 + u3**2)/(2*u3)}$
con:=concurrent(pp_line(O,H),pp_line(A,A1),pp_line(B,B1));
con := 0$
end;
Time for test: 114295 ms, plus GC time: 5438 ms