reduce-historical/r37/log/geometry.rlg

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