% GeoProver | Version 1.3a | Jan 20 2003
% Author: H.-G. Graebe, Univ. Leipzig, Germany
% http://www.informatik.uni-leipzig.de/~graebe
module geoprover;
comment
GeoProver inline part. Version 1.3
Data structures:
Point A :== {a1,a2} <=> A=(a1,a2)
Line a :== {a1,a2,a3} <=> a1*x+a2*y+a3 = 0
Circle c :== {c0,c1,c2,c3} <=> c0*(x^2+y^2)+c1*x+c2*y+c3 = 0
end comment;
put ('geoprover,'name," GeoProver ")$
put ('geoprover,'version," 1.3a ")$
put ('geoprover,'date," December 30, 2002")$
algebraic(write(" Geoprover ", get('geoprover,'version),
" Last update ",get('geoprover,'date)));
% ============= Handling non degeneracy conditions ===============
algebraic procedure clear_ndg; !*ndg!*:={};
algebraic procedure print_ndg; !*ndg!*;
algebraic procedure add_ndg(d);
if not member(d,!*ndg!*) then !*ndg!*:=d . !*ndg!*;
clear_ndg(); % Initialization
% ================= elementary geometric constructions ===============
% Generators:
algebraic procedure is_equal(a,b); a-b;
%algebraic procedure Normal(a); a;
algebraic procedure Point(a,b); {a,b};
algebraic procedure Line(a,b,c); reduce_coords({a,b,c});
algebraic procedure par_point(a,b,c);
Point(part(a,1)-part(b,1)+part(c,1),
part(a,2)-part(b,2)+part(c,2));
algebraic procedure pp_line(a,b);
% The line through A and B.
Line(part(b,2)-part(a,2),part(a,1)-part(b,1),
part(a,2)*part(b,1)-part(a,1)*part(b,2));
algebraic procedure intersection_point(a,b);
% The intersection point of the lines a,b.
begin scalar d,d1,d2;
d:=part(a,1)*part(b,2)-part(b,1)*part(a,2);
d1:=part(a,3)*part(b,2)-part(b,3)*part(a,2);
d2:=part(a,1)*part(b,3)-part(b,1)*part(a,3);
if d=0 then rederr"Lines are parallel";
add_ndg(num d);
return Point(-d1/d,-d2/d);
end;
algebraic procedure ortho_line(p,a);
% The line through P orthogonal to the line a.
begin scalar u,v; u:=first a; v:=second a;
return Line(v,-u,u*second p-v*first p);
end;
algebraic procedure par_line(p,a);
% The parallel to line a through P.
Line(part(a,1),part(a,2),
-(part(a,1)*part(p,1)+part(a,2) *part(p,2)));
algebraic procedure varpoint(b,a,l);
% The point D=l*A+(1-l)*B.
Point(l*part(a,1)+(1-l)*part(b,1),l*part(a,2)+(1-l)*part(b,2));
algebraic procedure line_slider(a,u);
% Slider on the line a using parameter u.
begin scalar p,d;
if part(a,2)=0 then
<< p:=Point(-part(a,3)/part(a,1),u); d:=part(a,1); >>
else
<< p:=Point(u,-(part(a,3)+part(a,1)*u)/part(a,2));
d:=part(a,2);
>>;
add_ndg(num d);
return p;
end;
algebraic procedure circle_slider(M,A,u);
% Slider on the circle with center M and circumfere point A using
% parameter u.
begin scalar a1,a2,m1,m2,d;
a1:=part(A,1); a2:=part(A,2); d:= u^2 + 1;
m1:=part(M,1); m2:=part(M,2);
add_ndg(num d);
return Point((a1*(u^2-1) + 2*m1 + 2*(m2-a2)*u)/d,
(a2 + 2*(m1-a1)*u + (2*m2-a2)*u^2)/d);
end;
algebraic procedure sqrdist(a,b);
% The square of the distance between the points A and B.
(part(b,1)-part(a,1))^2+(part(b,2)-part(a,2))^2;
% ================= elementary geometric properties ===============
algebraic procedure is_collinear(a,b,c);
% A,B,C are on a common line.
det mat((part(a,1),part(a,2),1),
(part(b,1),part(b,2),1),
(part(c,1),part(c,2),1));
algebraic procedure is_concurrent(a,b,c);
% Lines a,b,c have a common point.
det mat((part(a,1),part(a,2),part(a,3)),
(part(b,1),part(b,2),part(b,3)),
(part(c,1),part(c,2),part(c,3)));
algebraic procedure is_parallel(a,b);
% 0 <=> the lines a,b are parallel.
part(a,1)*part(b,2)-part(b,1)*part(a,2);
algebraic procedure is_orthogonal(a,b);
% 0 <=> the lines a,b are orthogonal.
part(a,1)*part(b,1)+part(a,2)*part(b,2);
algebraic procedure on_line(p,a);
% Substitute point P into the line a.
part(p,1)*part(a,1)+part(p,2)*part(a,2)+part(a,3);
algebraic procedure eq_dist(a,b,c,d); sqrdist(a,b)-sqrdist(c,d);
% #########################################
% # #
% # Non linear geometric objects #
% # #
% #########################################
% ===================== angles
algebraic procedure l2_angle(a,b);
% tan of the angle between the lines a and b.
begin scalar d; d:=(part(a,1)*part(b,1)+part(a,2)*part(b,2));
add_ndg(num(d));
return (part(a,2)*part(b,1)-part(b,2)*part(a,1))/d;
end;
algebraic procedure angle_sum(a,b);
% a=tan(\alpha), b=tan(\beta). Returns tan(\alpha+\beta)
begin scalar d; d:=(1-a*b); add_ndg(num d);
return (a+b)/d;
end;
algebraic procedure eq_angle(a,b,c,d,e,f);
p3_angle(a,b,c)-p3_angle(d,e,f);
algebraic procedure on_bisector(P,A,B,C);
% P is a point on the bisector of the angle ABC.
% Returns num(u)*den(v)-num(v)*den(u) with
% u:=angle(pp_line(A,B),pp_line(P,B))
% v:=angle(pp_line(P,B),pp_line(C,B))
begin scalar a1,a2,b1,b2,c1,c2,p1,p2;
a1:=part(A,1); a2:=part(A,2);
b1:=part(b,1); b2:=part(b,2);
c1:=part(c,1); c2:=part(c,2);
p1:=part(p,1); p2:=part(p,2);
return ( - a1*b2 + a1*p2 + a2*b1 - a2*p1 - b1*p2 + b2*p1)*(b1^2 -
b1*c1 - b1*p1 + b2^2 - b2*c2 - b2*p2 + c1*p1 + c2*p2) - (a1*b1 -
a1*p1 + a2*b2 - a2*p2 - b1^2 + b1*p1 - b2^2 + b2*p2)*(b1*c2 -
b1*p2 - b2*c1 + b2*p1 + c1*p2 - c2*p1)
end;
algebraic procedure rotate(C, A, angle);
begin scalar ac1,ac2;
ac1:=part(A,1)-part(C,1); ac2:=part(A,2)-part(C,2);
return Point(part(C,1)+ac1*cos(angle*pi)-ac2*sin(angle*pi),
part(C,2)+ac1*sin(angle*pi)+ac2*cos(angle*pi));
end;
% ========== symmetric lines and points
algebraic procedure sym_line(a,l);
% The line symmetric to a wrt. the line l.
begin scalar a1,a2,a3,l1,l2,l3,u;
a1:=part(a,1); a2:=part(a,2); a3:=part(a,3);
l1:=part(l,1); l2:=part(l,2); l3:=part(l,3);
u:=l1^2 - l2^2;
return Line(- a1*u - 2*a2*l1*l2, - 2*a1*l1*l2 + a2*u,
- 2*(a1*l1 + a2*l2)*l3 + a3*(l1^2 + l2^2));
end;
% ===================== circles
algebraic procedure Circle(c1,c2,c3,c4); reduce_coords({c1,c2,c3,c4});
algebraic procedure pc_circle(M,A);
% Circle with center M and Point A on circumference.
Circle(1, -2*part(M,1), -2*part(M,2),
part(A,1)*(2*part(M,1)-part(A,1)) +
part(A,2)*(2*part(M,2)-part(A,2)));
algebraic procedure circle_center c;
% The center of the circle c.
begin add_ndg(num part(c,1));
return Point(-part(c,2)/2/part(c,1) ,-part(c,3)/(2*part(c,1)));
end;
algebraic procedure circle_sqradius c;
% The squared radius of the circle c.
begin add_ndg(num part(c,1));
return
((part(c,2)^2+part(c,3)^2) - 4*part(c,4)*part(c,1)) /
(2*part(c,1))^2;
end;
algebraic procedure p3_circle(A,B,C);
% The circle through three given points
begin scalar a1,a2,a3,b1,b2,b3,c1,c2,c3;
a1:=part(A,1); a2:=part(A,2); a3:=a1^2+a2^2;
b1:=part(b,1); b2:=part(b,2); b3:=b1^2+b2^2;
c1:=part(c,1); c2:=part(c,2); c3:=c1^2+c2^2;
return Circle(a1*(b2-c2) + (a2-b2)*c1 + b1*(c2-a2),
a3*(c2-b2) + (a2-c2)*b3 + (b2-a2)*c3,
a3*(b1-c1) + (c1-a1)*b3 + (a1-b1)*c3,
a3*(b2*c1-b1*c2) + (a1*c2-a2*c1)*b3 + (a2*b1-a1*b2)*c3)
end;
algebraic procedure on_circle(P,c);
begin scalar p1,p2; p1:=part(P,1); p2:=part(P,2);
return part(c,1)*(p1^2+p2^2)+part(c,2)*p1+part(c,3)*p2+part(c,4);
end;
% Intersecting with circles
algebraic procedure other_cl_point(P,c,l);
% circle c and line l intersect at P. The procedure returns their
% second intersection point.
if on_line(P,l) neq 0 then rederr "Point not on the line"
else if on_circle(P,c) neq 0 then
rederr "Point not on the circle"
else begin scalar c1,c2,c3,l1,l2,d,d1,p1,p2;
c1:=part(c,1); c2:=part(c,2); c3:=part(c,3);
l1:=part(l,1); l2:=part(l,2); p1:=part(P,1); p2:=part(P,2);
d:=c1*(l1^2 + l2^2); add_ndg(num d); d1:=c1*(l1^2-l2^2);
return {(d1*p1+((2*c1*p2 + c3)*l1-c2*l2)*l2)/d,
(- d1*p2+((2*c1*p1 + c2)*l2-c3*l1)*l1)/d};
end;
algebraic procedure radical_axis(c1,c2);
% Radical axis of the circles c1 and c2, i.e. the line through the
% intersection points of the two circles if they intersect.
for i:=2:4 collect
(part(c1,1)*part(c2,i)-part(c1,i)*part(c2,1));
algebraic procedure other_cc_point(P,c1,c2);
% Circles c1 and c2 intersect at P. The procedure returns their
% second intersection point.
other_cl_point(P,c1,radical_axis(c1,c2));
algebraic procedure is_cl_tangent(c,l);
% Line l is tangent to the circle c.
begin scalar c1,c2,c3,c4,l1,l2,l3;
c1:=part(c,1); c2:=part(c,2); c3:=part(c,3); c4:=part(c,4);
l1:=part(l,1); l2:=part(l,2); l3:=part(l,3);
return - 4*c1^2*l3^2 + 4*c1*c2*l1*l3 + 4*c1*c3*l2*l3 -
4*c1*c4*l1^2 - 4*c1*c4*l2^2 + c2^2*l2^2 - 2*c2*c3*l1*l2 +
c3^2*l1^2
end;
algebraic procedure is_cc_tangent(c,d);
% Two circles c,d are tangent.
begin scalar c1,c2,c3,c4,d1,d2,d3,d4;
c1:=part(c,1); c2:=part(c,2); c3:=part(c,3); c4:=part(c,4);
d1:=part(d,1); d2:=part(d,2); d3:=part(d,3); d4:=part(d,4);
return
4*c1^2*d4^2 - 4*c1*c2*d2*d4 - 4*c1*c3*d3*d4 - 8*c1*c4*d1*d4 +
4*c1*c4*d2^2 + 4*c1*c4*d3^2 + 4*c2^2*d1*d4 - c2^2*d3^2 + 2*c2*c3*d2*d3
- 4*c2*c4*d1*d2 + 4*c3^2*d1*d4 - c3^2*d2^2 - 4*c3*c4*d1*d3 +
4*c4^2*d1^2
end;
% ============= some additional tools ===============
symbolic operator list2mat;
symbolic procedure list2mat u;
'mat. for each x in cdr reval u collect cdr x;
algebraic procedure extractmat(polys,vars);
% extract the coefficient matrix from the linear system polys.
begin
if length polys neq length vars then
rederr"Number of variables doesn't match";
for each p in polys do for each x in vars do
if deg(p,x)>1 then rederr"Equations not of linear type";
return list2mat
for each x in vars collect
for each p in polys collect coeffn(p,x,1);
end;
algebraic procedure reduce_coords u;
% Divide out the content of homogeneous coordinates.
begin scalar l,g;
l:=den first u; g:=num first u;
for each x in rest u do <<l:=lcm(l,den x); g:=gcd(g,num x) >>;
add_ndg(g);
return for each x in u collect (x*l/g);
end;
% ================ new
algebraic procedure circle_inverse(M,R,P);
% compute the inverse of P wrt. the circle pc_circle(M,R)
begin scalar m1,m2,r1,r2,p1,p2,d;
m1:=part(M,1); m2:=part(M,2);
r1:=part(R,1); r2:=part(R,2);
p1:=part(P,1); p2:=part(P,2);
d:=(m1-p1)^2+(m2-p2)^2;
add_ndg(d);
return ((m1-p1)^2+(m2-p2)^2+(m1-r1)^2+(m2-r2)^2)/d;
end;
% GeoProver code generated from database
algebraic procedure altitude(A__,B__,C__);
ortho_line(A__,pp_line(B__,C__));
algebraic procedure centroid(A__,B__,C__);
intersection_point(median(A__,B__,C__),median(B__,C__,A__));
algebraic procedure circumcenter(A__,B__,C__);
intersection_point(p_bisector(A__,B__), p_bisector(B__,C__));
algebraic procedure csym_point(P__,Q__);
varpoint(Q__,P__,-1);
algebraic procedure fixedpoint(A__,B__,u_);
varpoint(A__,B__,u_);
algebraic procedure is_concyclic(A__,B__,C__,D__);
on_circle(D__,p3_circle(A__,B__,C__));
algebraic procedure median(A__,B__,C__);
pp_line(A__,midpoint(B__,C__));
algebraic procedure midpoint(A__,B__);
fixedpoint(A__,B__,1/2);
algebraic procedure orthocenter(A__,B__,C__);
intersection_point(altitude(A__,B__,C__),altitude(B__,C__,A__));
algebraic procedure other_incenter(M__,A__,B__);
intersection_point(ortho_line(A__,pp_line(M__,A__)),
ortho_line(B__,pp_line(M__,B__)));
algebraic procedure p3_angle(A__,B__,C__);
l2_angle(pp_line(B__,A__),pp_line(B__,C__));
algebraic procedure p9_center(A__,B__,C__);
circle_center(p9_circle(A__,B__,C__));
algebraic procedure p9_circle(A__,B__,C__);
p3_circle(midpoint(A__,B__),midpoint(A__,C__),midpoint(B__,C__));
algebraic procedure p_bisector(B__,C__);
ortho_line(midpoint(B__,C__),pp_line(B__,C__));
algebraic procedure pappus_line(A__,B__,C__,D__,E__,F__);
pp_line(intersection_point(pp_line(A__,E__),pp_line(B__,D__)),
intersection_point(pp_line(A__,F__),pp_line(C__,D__)));
algebraic procedure pedalpoint(P__,a_);
intersection_point(ortho_line(P__,a_),a_);
algebraic procedure sqrdist_pl(A__,l_);
sqrdist(A__,pedalpoint(A__,l_));
algebraic procedure sym_point(P__,l_);
fixedpoint(P__,pedalpoint(P__,l_),2);
algebraic procedure triangle_area(A__,B__,C__);
1/2*is_collinear(A__,B__,C__);
endmodule; % GeoProver
end;