(*  GeoBeweiser von Volker Messerschmidt (Juli 2004)*)
(*  http://www.vm4.de/diplom/ *)
(*  abgeleitet von *)
(*  GeoProver | Version 1.3a | Jan 20 2003 *)
(*  Author: H.-G. Graebe, Univ. Leipzig, Germany *)
(*  http://www.informatik.uni-leipzig.de/~graebe *)

BeginPackage["GeoBeweiser`"]

(* GeoBeweiser *)

point::usage:="Punktkonstruktor";
line::usage:="Linienkonstruktor";
circle::usage:="Kreiskonstruktor";

altitude::usage:="Entfernung von A zur Geraden BC";
anglesum::usage:="Addition der Winkel tan(alpha) und tan(beta) zu tan(alpha+beta)";
centroid::usage:="Schwerpunkt des Dreiecks ABC, Schnittpunkt der Mittellinien";
circlecenter::usage:="Mittelpunkt des Kreises";
circleslider::usage:="Wähle einen Punkt auf dem Kreis mit Mittelpunkt M und Randpunkt A";
circlesqradius::usage:="Quadrat des Radius des Kreises";
circumcenter::usage:="Mittelpunkt des Umkreis des Dreiecks ABC";
csympoint::usage:="Punkt symmetrisch zu P mit Q als Symmetriepunkt";
eqangle::usage:="Test auf gleichgroße Winkel w(ABC) = w(DEF)";
eqdist::usage:="Test auf gleiche Distanz d(AB) = d(CD)";
fixedpoint::usage:="Fester Punkt auf der Linie AB. u ist ein fester Wert";
intersectionpoint::usage:="Schnittpunkt zweier Linien";
iscctangent::usage:="Test, ob zwei Kreise tangential liegen";
iscltangent::usage:="Test, ob eine Linie tangential zum Kreis liegt";
iscollinear::usage:="Test, ob drei Punkte auf einer Geraden liegen";
isconcurrent::usage:="Test, ob drei Geraden durch einen Punkt gehen";
isconcyclic::usage:="Test, ob vier Punkte auf einem Kreis liegen";
isequal::usage:="Test auf Gleichheit von zwei Objekten";
isorthogonal::usage:="Test, ob zwei Geraden orthogonal sind";
isparallel::usage:="Test, ob zwei Geraden parallel sind";
l2angle::usage:="Tangens des Winkels zwischen den Geraden a und b";
lineslider::usage:="Wählt einen Punkt auf einer Geraden";
median::usage:="Gerade von A zum Mittelpunkt von BC";
midpoint::usage:="Mittelpunkt zwischen A und B";
onbisector::usage:="Testet, ob P ein Punkt auf der (inneren oder äußeren) Winkelhalbierenden von ABC ist";
oncircle::usage:="Test, ob ein Punkt auf einem Kreis liegt";
online::usage:="Test, ob ein Punkt auf der Geraden liegt";
ortholine::usage:="Gerade durch P orthogonal zu a";
orthocenter::usage:="Höhenschnittpunkt des Dreiecks ABC";
otherccpoint::usage:="Zwei Kreise schneiden sich im Punkt P. Ergibt den zweiten Schnittpunkt";
otherclpoint::usage:="Gerade und Kreis schneiden sich im Punkt P. Ergibt den zweiten Schnittpunkt";
otherincenter::usage:="M ist der Inkreismittelpunkt von ABC. Gibt den Ankreismittelpunkt von ABC auf der Winkelhalbierenden CM aus.";
p3angle::usage:="Tangens des Winkels zwischen BA und BC";
p3circle::usage:="Kreis durch drei gegebene Punkte";
p9center::usage:="Mittelpunkt des Neunpunktekreises des Dreiecks ABC";
p9circle::usage:="Neunpunktekreis des Dreiecks ABC";
pappusline::usage:="Pappuslinie von einem kegelförmigen 6-Tupel von Punkten";
parline::usage:="Gerade durch den Punkt P parallel zur Geraden a";
parpoint::usage:="Punkt D, der ABCD zu einem einem Parallelogramm macht";
pbisector::usage:="Mittelsenkrechte von BC";
pccircle::usage:="Kreis mit Mittelpunkt M und Randpunkt A";
pedalpoint::usage:="Fusspunkt der Senkrechten von Punkt P auf die Gerade a";
ppline::usage:="Gerade durch A und B";
radicalaxis::usage:="Die Radikalachse von zwei Kreisen";
rotate::usage:="Rotiere Punkt A um den Punkt C, mit angle*Pi (für angle nur Konstanten erlaubt)";
sqrdist::usage:="Quadrat der Entfernung von Punkt A zu Punkt B";
sqrdistpl::usage:="Quadrat der Entfernung von Punkt A zu der Gerade l";
symline::usage:="Spiegelung der Geraden a an l";
sympoint::usage:="Punkt P gespiegelt an der Geraden l";
trianglearea::usage:="Flächeninhalt des Dreiecks ABC";
varpoint::usage:="Variabler Punkt auf der Geraden AB mit u als Parameter";

clearndg::usage:="degeneracy tracing";
printndg::usage:="degeneracy tracing";
addndg::usage:="degeneracy tracing";

(* Neue Funktionen (nicht in GeoProver enthalten)*)
dreiteilung::usage="Winkeldreiteilung (Vollerwinkel, Drittelwinkel)";
IdealMember::usage="Test, ob ein Polynom im Radikal eines Ideals enthalten ist";
AddGleichung::usage="Sammelt die Gleichung ein für die Funktion Gleichungen";
Gleichungen::usage="Gibt die gesammelten Gleichungen zurück";
LoeschGleichungen::usage="Löscht die Gleichungen"

(* GeoProver Inline part Version 1.3a *)

intersectionpoint::prll:="Die Geraden sind parallel."
otherclpoint::pnol:="Punkt ist nicht auf der Geraden."
otherclpoint::pnoc:="Punkt ist nicht auf dem Kreis."

Begin["Private`"] 

AddGleichungErg={};
AddGleichung[gl_]:=Module[{},
  AddGleichungErg = Append[AddGleichungErg, gl];
  gl
  ]
AddGleichung[gl_GPPoint]:=Module[{},
  AddGleichungErg = Append[AddGleichungErg, gl[[1]]];
  AddGleichungErg = Append[AddGleichungErg, gl[[2]]];
  gl
  ]
Gleichungen:=AddGleichungErg
LöschGleichungen:=Module[{},AddGleichungErg={};]

IdealMember[con_, polys_,vars_] := 
  GroebnerBasis[Flatten[{polys, 1 - y con}], vars, 
    CoefficientDomain -> RationalFunctions, MonomialOrder -> DegreeReverseLexicographic]=={1}

Unprotect[Plus]
GPPoint /: Plus[x_GPPoint, y_GPPoint] := 
  point[x[[1]] + y[[1]], x[[2]] + y[[2]]]
Protect[Plus]

Unprotect[Times]
GPPoint /: Times[x_GPPoint,-1] := 
  point[-x[[1]], -x[[2]]]
Protect[Times]

dreiteilung[ganzerw_, drittelw_]:=Together[ganzerw (3 * drittelw^2-1) - drittelw^3 + 3 * drittelw]
det3=Det[{{a11,a12,a13},{a21,a22,a23},{a31,a32,a33}}]
(*Accelarates computation of iscollinear and isconcurrent*) 
 
(*Handling non degeneracy conditions*)

clearndg[]:=ndglist={};


printndg[]:=Print[ndglist];

addndg[d_]:=
  Module[{tmplist},tmplist=ndglist;
    If[!MemberQ[tmplist,d],ndglist=Append[tmplist,d]]]

reduceCoords[u_]:=Module[{p,l,g,i,x},
  p=Part[u,1]; l=Denominator[p]; g=Numerator[p]; 
    For[i=2,i<=Length[u],i++,
      p=Part[u,i]; 
      l=PolynomialLCM[l,Denominator[p]]; 
      g=PolynomialGCD[g,Numerator[p]]];
  addndg[g];
  u/.x_->Together[x*l/g]]

clearndg[]  (* Initialization *)

(* How geometric objects should be formatted on output *)

Format[GPPoint[A___]]:={A}
Format[GPLine[A___]]:={A}
Format[GPCircle[A___]]:={A}

(* elementary geometric constructions *)

(* Generators *)

point[a_,b_]:=GPPoint@@Together/@{a,b};

line[a_,b_,c_]:=GPLine@@(Together/@{a,b,c}//reduceCoords); 

isequal[a_,b_]:=(a-b)//Together;

parpoint[a_GPPoint,b_GPPoint,c_GPPoint]:= 
  point[Part[a,1]-Part[b,1]+Part[c,1], Part[a,2]-Part[b,2]+Part[c,2]]

ppline[a_GPPoint,b_GPPoint]:=
  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]]

intersectionpoint[a_GPLine,b_GPLine]:= Module[{d,d1,d2,p},
  d=Together[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, Message[intersectionpoint::prll];Return[]];
  addndg[Numerator[d]];
  point[-d1/d,-d2/d]]                   

ortholine[p_GPPoint,a_GPLine]:= Module[{u,v},
  u=Part[a,1];v=Part[a,2];
  line[v,-u,u*Part[p,2]-v*Part[p,1]]]

parline[p_GPPoint,a_GPLine]:=
  line[Part[a,1],Part[a,2],-(Part[a,1]*Part[p,1]+Part[a,2]*Part[p,2])]

varpoint[b_GPPoint,a_GPPoint,l_]:=
  point[l*Part[a,1]+(1-l)*Part[b,1],l*Part[a,2]+(1-l)*Part[b,2]]

lineslider[a_GPLine,u_]:= Module[{p,d},
  If[Part[a,2]===0,p=point[-Part[a,3]/Part[a,1],u]; 
      d=Part[a,1], p=point[u,-(Part[a,3]+Part[a,1]*u)/Part[a,2]]; 
      d=Part[a,2]];
  addndg[Numerator[d]];
  p]

sqrdist[a_GPPoint,b_GPPoint]:=
  Together[(Part[b,1]-Part[a,1])^2+(Part[b,2]-Part[a,2])^2]

(*elementary geometric properties*)

iscollinear[a_GPPoint,b_GPPoint,c_GPPoint]:=
    Together[det3/.{a11->Part[a,1],a12->Part[a,2],a13->1,
    a21->Part[b,1],a22->Part[b,2],a23->1,
    a31->Part[c,1],a32->Part[c,2],a33->1}]

isconcurrent[a_GPLine,b_GPLine,c_GPLine]:=
    Together[det3/.{a11->Part[a,1],a12->Part[a,2],a13->Part[a,3],
    a21->Part[b,1],a22->Part[b,2],a23->Part[b,3],
    a31->Part[c,1],a32->Part[c,2],a33->Part[c,3]}]

isparallel[a_GPLine,b_GPLine]:=
  Together[Part[a,1]*Part[b,2]-Part[b,1]*Part[a,2]]

isorthogonal[a_GPLine,b_GPLine]:=
  Together[Part[a,1]*Part[b,1]+Part[a,2]*Part[b,2]]

online[p_GPPoint,a_GPLine]:=
  Together[Part[p,1]*Part[a,1]+Part[p,2]*Part[a,2]+Part[a,3]]

eqdist[a_GPPoint,b_GPPoint,c_GPPoint,d_GPPoint]:=
  Together[sqrdist[a,b]-sqrdist[c,d]]

(*      Non linear geometric objects                 *)

(*angles*)

l2angle[a_GPLine,b_GPLine]:= Module[{d}, 
  d:=Together[Part[a,1]*Part[b,1]+Part[a,2]*Part[b,2]]; 
  addndg[Numerator[d]];
  Together[(Part[a,2]*Part[b,1]-Part[b,2]*Part[a,1])/d]] 

p3angle[a_GPPoint,b_GPPoint,c_GPPoint]:=l2angle[ppline[b,a],ppline[b,c]] 

anglesum[a_,b_]:= Module[{d},
  d:=Together[1-a*b]; addndg[Numerator[d]]; 
  Together[(a+b)/d]]

eqangle[a_GPPoint,b_GPPoint,c_GPPoint,d_GPPoint,e_GPPoint,f_GPPoint]:=
  Together[p3angle[a,b,c]-p3angle[d,e,f]]

onbisector[p_GPPoint,a_GPPoint,b_GPPoint,c_GPPoint]:=
Module[{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];
  Together[
    Together[ - a1*b2 + a1*p2 + a2*b1 - a2*p1 - b1*p2 + b2*p1]*
    Together[ b1^2 - b1*c1 - b1*p1 + b2^2 - b2*c2 - b2*p2 + c1*p1
        + c2*p2]  
    - Together[a1*b1 - a1*p1 + a2*b2 - a2*p2 - b1^2 + b1*p1 - b2^2
        + b2*p2] *
    Together[b1*c2 -  b1*p2 - b2*c1 + b2*p1 + c1*p2 - c2*p1]]]  

rotate[C_GPPoint, A_GPPoint, angle_]:= Module[{ac1,ac2},
  ac1:=Part[A,1]-Part[C,1]; ac2:=Part[A,2]-Part[C,2];
  point[Part[C,1]+ac1*Cos[angle*\[Pi]]-ac2*Sin[angle*\[Pi]],
    Part[C,2]+ac1*Sin[angle*\[Pi]]+ac2*Cos[angle*\[Pi]]]]

(*symmetric lines and points*)

symline[a_GPLine,l_GPLine]:= Module[{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:=Together[l1^2 - l2^2];
  line[- a1*u - 2*a2*l1*l2, - 2*a1*l1*l2 + a2*u, 
        - 2*(a1*l1 + a2*l2)*l3 + a3*(l1^2 + l2^2)]]

(* circles *)

circleslider[m_GPPoint,a_GPPoint,u_]:= Module[{d,a1,a2,m1,m2},
  a1:=Part[a,1]; a2:=Part[a,2];
  m1:=Part[m,1]; m2:=Part[m,2];
  d:=Together[u^2+1]; addndg[Numerator[d]]; 
  point[(a1*(u^2-1) + 2*m1 + 2*(m2-a2)*u)/d, 
     (a2 + 2*(m1-a1)*u + (2*m2-a2)*u^2)/d]]  

circle[c1_,c2_,c3_,c4_]:=
  GPCircle@@(Together/@{c1,c2,c3,c4}//reduceCoords); 

pccircle[m_GPPoint,a_GPPoint]:=Module[{d,a1,a2,m1,m2},
  a1:=Part[a,1]; a2:=Part[a,2];
  m1:=Part[m,1]; m2:=Part[m,2];
  circle[1, -2*m1, -2*m2, a1*(2*m1-a1) + a2*(2*m2-a2)]]

circlecenter[c_GPCircle]:=(addndg[Numerator[Part[c,1]]];
  point[-Part[c,2]/2/Part[c,1] ,-Part[c,3]/(2*Part[c,1])])

circlesqradius[c_GPCircle]:= (addndg[Numerator[Part[c,1]]];
  Together[((Part[c,2]^2+Part[c,3]^2) - 4*Part[c,4]*Part[c,1]) / 
    (2*Part[c,1])^2])

p3circle[a_GPPoint,b_GPPoint,c_GPPoint]:= 
Module[{a1,a2,a3,b1,b2,b3,c1,c2,c3},
  a1:=Part[a,1]; a2:=Part[a,2]; a3:=Together[a1^2+a2^2];
  b1:=Part[b,1]; b2:=Part[b,2]; b3:=Together[b1^2+b2^2];
  c1:=Part[c,1]; c2:=Part[c,2]; c3:=Together[c1^2+c2^2];
  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]]
 
oncircle[p_GPPoint,c_GPCircle]:=Module[{p1,p2}, 
  p1:=Part[p,1]; p2:=Part[p,2];
  Together[Part[c,1]*(p1^2+p2^2)+Part[c,2]*p1+Part[c,3]*p2+Part[c,4]]]  

(*Intersections with circles*)

otherclpoint[p_GPPoint,c_GPCircle,l_GPLine]:= 
Module[{c1,c2,c3,l1,l2,d,d1,p1,p2},
  If[ !(online[p,l] === 0), 
      Message[otherclpoint::pnol]; Return[],
      If[!(oncircle[p,c] === 0),
          Message[otherclpoint::pnoc]; Return[], 
          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:=Together[c1*(l1^2 + l2^2)]; addndg[Numerator[d]]; 
          d1:=Together[c1*(l1^2-l2^2)];
          point[(d1*p1+((2*c1*p2 + c3)*l1-c2*l2)*l2)/d, 
            (- d1*p2+((2*c1*p1 + c2)*l2-c3*l1)*l1)/d]]]]  

radicalaxis[c1_GPCircle,c2_GPCircle]:=
  Apply[line,Table[Together[
    Part[c1,1]*Part[c2,i]-Part[c1,i]*Part[c2,1]], {i,2,4}]]; 

otherccpoint[p_GPPoint,c1_GPCircle,c2_GPCircle]:= 
  otherclpoint[p,c1,radicalaxis[c1,c2]]

iscltangent[c_GPCircle,l_GPLine]:=Module[{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]; 
  Together[4*c1*( (- c1*l3 + c2*l1 + c3*l2)*l3 - c4*(l1^2 + l2^2))  
    + (c2*l2 - c3*l1)^2]]

iscctangent[c_GPCircle,d_GPCircle]:=Module[{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]; 
  Together[4*(c1*d4-c4*d1)^2 - (c2*d3-c3*d2)^2 
    - 4*(c1*d4+c4*d1)*(c2*d2+c3*d3) + 4*c1*c4*(d2^2 + d3^2) 
    + 4*(c2^2 + c3^2)*d1*d4]] 

(* During 'Build' some more generated code will be added  *)

(* GeoProver code generated from database *)

altitude[$A_GPPoint,$B_GPPoint,$C_GPPoint]:=
  ortholine[$A,ppline[$B,$C]];

centroid[$A_GPPoint,$B_GPPoint,$C_GPPoint]:=
  intersectionpoint[median[$A,$B,$C],median[$B,$C,$A]];

circumcenter[$A_GPPoint,$B_GPPoint,$C_GPPoint]:=
  intersectionpoint[pbisector[$A,$B], pbisector[$B,$C]];

csympoint[$P_GPPoint,$Q_GPPoint]:=
  varpoint[$Q,$P,-1];

fixedpoint[$A_GPPoint,$B_GPPoint,$u_]:=
  varpoint[$A,$B,$u];

isconcyclic[$A_GPPoint,$B_GPPoint,$C_GPPoint,$D_GPPoint]:=
  oncircle[$D,p3circle[$A,$B,$C]];

median[$A_GPPoint,$B_GPPoint,$C_GPPoint]:=
  ppline[$A,midpoint[$B,$C]];

midpoint[$A_GPPoint,$B_GPPoint]:=
  fixedpoint[$A,$B,1/2];

orthocenter[$A_GPPoint,$B_GPPoint,$C_GPPoint]:=
  intersectionpoint[altitude[$A,$B,$C],altitude[$B,$C,$A]];

otherincenter[$M_GPPoint,$A_GPPoint,$B_GPPoint]:=
  intersectionpoint[ortholine[$A,ppline[$M,$A]],
ortholine[$B,ppline[$M,$B]]];

p3angle[$A_GPPoint,$B_GPPoint,$C_GPPoint]:=
  l2angle[ppline[$B,$A],ppline[$B,$C]];

p9center[$A_GPPoint,$B_GPPoint,$C_GPPoint]:=
  circlecenter[p9circle[$A,$B,$C]];

p9circle[$A_GPPoint,$B_GPPoint,$C_GPPoint]:=
  p3circle[midpoint[$A,$B],midpoint[$A,$C],midpoint[$B,$C]];

pbisector[$B_GPPoint,$C_GPPoint]:=
  ortholine[midpoint[$B,$C],ppline[$B,$C]];

pappusline[$A_GPPoint,$B_GPPoint,$C_GPPoint,$D_GPPoint,$E_GPPoint,$F_GPPoint]:=
  ppline[intersectionpoint[ppline[$A,$E],ppline[$B,$D]],
    intersectionpoint[ppline[$A,$F],ppline[$C,$D]]];

pedalpoint[$P_GPPoint,$a_GPLine]:=
  intersectionpoint[ortholine[$P,$a],$a];

sqrdistpl[$A_GPPoint,$l_GPLine]:=
  sqrdist[$A,pedalpoint[$A,$l]];

sympoint[$P_GPPoint,$l_GPLine]:=
  fixedpoint[$P,pedalpoint[$P,$l],2];

trianglearea[$A_GPPoint,$B_GPPoint,$C_GPPoint]:=
  1/2*iscollinear[$A,$B,$C];

End[] (* Private *)
EndPackage[]
