This Notebook contains supplementary material related to the paper
Stanley Rabinowitz, “Relationships Between Six Incircles”, Sangaku Journal of Mathematics, 3(2019)51–66
which can be found at http://www.sangaku-journal.eu/2019/SJM_2019_51-66_Rabinowitz.pdf.

This Notebook is Copyright © August 2019 by Stanley Rabinowitz (stan.rabinowitz@comcast.net) and is published with open access. This Notebook is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author and the source are credited.
------------------------------------------------------------------------------

Proof of Theorem 7.3

(* ================ Set up equations ================= *)

AB=Sin[c+d];
AC=Sin[a+b];
BC=Sin[a+b+c+d];
AF=Sin[a+b]Sin[d]/Sin[a+b+c];
BF=BC Sin[c]/Sin[a+b+c];
AE=Sin[c+d]Sin[a]/Sin[b+c+d];
CE=BC Sin[b]/Sin[b+c+d];
PB=BC Sin[c]/Sin[b+c];
PC=BC Sin[b]/Sin[b+c];
PE=CE Sin[d]/Sin[b+c];
PF=BF Sin[a]/Sin[b+c];
BD=BF AE BC/(BF AE+AF CE);
CD=AF CE BC/(BF AE+AF CE);
PA=Sqrt[AB^2+PB^2-2AB PB Cos[a]];
PD=PA BD CE / (BC AE);
K1=(1/2)PB*BD*Sin[b];
K2=(1/2)PC*CD*Sin[c];
K3=(1/2)PC*CE*Sin[d];
K4=(1/2)AE*PE*Sin[b+c+d];
K5=(1/2)AF*PF*Sin[a+b+c];
K6=(1/2)PB*BF*Sin[a];
s1=(PB+BD+PD)/2;
s2=(PD+CD+PC)/2;
s3=(PC+CE+PE)/2;
s4=(PE+AE+PA)/2;
s5=(PA+AF+PF)/2;
s6=(PF+BF+PB)/2;
r1=K1/s1;
r2=K2/s2;
r3=K3/s3;
r4=K4/s4;
r5=K5/s5;
r6=K6/s6;
rad=Pi/180;

(* Proof of the relationship: *)

{nr1,nr2,nr3,nr4,nr5,nr6}={r1,r2,r3,r4,r5,r6}/.{a->t,b->2t,c->60rad-2t,d->30rad-t};
expr=1/nr1-1/nr2-1/nr3+1/nr4-1/nr5+1/nr6;
Simplify[expr]

(* Proof of the relationship: *)

{nr1,nr2,nr3,nr4,nr5,nr6}={r1,r2,r3,r4,r5,r6}/.{a->90rad-t,b->60rad-t,c->t,d->2t};
expr=1/nr1-1/nr2-1/nr3+1/nr4-1/nr5+1/nr6;
Simplify[expr]

(* Proof of the relationship: *)

{nr1,nr2,nr3,nr4,nr5,nr6}={r1,r2,r3,r4,r5,r6}/.{a->2t,b->t,c->3t,d->90rad-t};
expr=1/nr1-1/nr2-1/nr3+1/nr4-1/nr5+1/nr6;
Simplify[expr]

.

(* Proof of the relationship: *)

{nr1,nr2,nr3,nr4,nr5,nr6}={r1,r2,r3,r4,r5,r6}/.{a->t,b->t,c->90rad-t,d->45rad};
expr=1/nr1-1/nr2-1/nr3+1/nr4-1/nr5+1/nr6;
Simplify[expr]

(* Proof of the relationship: *)

{nr1,nr2,nr3,nr4,nr5,nr6}={r1,r2,r3,r4,r5,r6}/.{a->3t,b->30rad-t,c->90rad-3t,d->6t};
expr=1/nr1-1/nr2-1/nr3+1/nr4-1/nr5+1/nr6;
Simplify[expr]

(* ====================================================================== *)

(* How these relationships were discovered. *)

(* Look for inverse relationship *)
FindInverseRelationships[step_]:=
(
count=0;
pcount=0;
type1=0; (* 60 degree angle and incenter *)
solutions={};
paramlist={};
max=155;
For[ia=step,ia<max,ia+=step,
For[ib=step,ib<max,ib+=step,
(* want c≥b *)
For[ic=ib,ic<max,ic+=step,
If[ia+ib+ic≥180,Break[]];
For[id=step,id<max,id+=step,
If[ia+ib+ic+id≥180,Break[]];
(* omit isosceles triangles *)
If[ia==id&&ib==ic,Continue[]];
(* omit where P is the incenter and one angle is 60º *)
If[ia==ib&&ic==30&&id==30,type1++;Continue[]];
If[ic==id&&ia==30&&ib==30,type1++;Continue[]];
If[ia==ib&&ic==id&&ib+ic==60,type1++;Continue[]];
{nr1,nr2,nr3,nr4,nr5,nr6}={r1,r2,r3,r4,r5,r6}/.{a->ia rad,b->ib rad,c->ic rad,d->id rad};
ans=Quiet[FindIntegerNullVector[N[{1/nr1,1/nr2,1/nr3,1/nr4,1/nr5,1/nr6},100],50,WorkingPrecision->100]];
If[Head[ans]=!=FindIntegerNullVector,
(* omit relationships that do not involve all 6 ri *)
If[Count[ans,1|-1]≠6,Continue[]];
theta=N[ArcSin[PF Sin[ia rad+ib rad+ic rad]/PA/.{a->ia rad,b->ib rad,c->ic rad,d->id rad}]];
itheta=Chop[N[theta/rad]];
itheta2=180-(ia+ib+ic+id+itheta);
(* omit cases where P is the circumcenter *)
If[ia==itheta&&ib==ic&&id==itheta2,Continue[]];
(* Normalize so first coefficient is positive *)
If[ans[[1]]<0,ans=-ans];
Print["a=",ia," b=",ib," c=",ic," d=",id," theta=",itheta," theta2=",itheta2," A=",itheta+itheta2,"   relationship=",ans];
count++;
AppendTo[solutions,{ia,ib,ic,id,itheta,itheta2,ans}];
];
]
]
]
];
Print["count=",count,"  type1=",type1];
(* Now go through pairs of solutions and see if any produce a third soltuion *)
For[i=1, i<count, i++,
For[j=i+1,j≤count,j++,
sol1full=solutions[[i]];
sol2full=solutions[[j]];
vector=sol1full[[7]];
If[sol2full[[7]]≠vector,Continue[]];
sol1=Drop[sol1full,-3];
sol2=Drop[sol2full,-3];
(* Assume the 2 vectors start an arithmetic progression *)
(* 2sol2-sol1 is the 3rd vector in the arithmetic progression *)
dif=sol1-sol2;
(* Find the arithmetic progression vector *)
sol3=sol1-tt*dif;
{ax,bx,cx,dx}=sol3;
{nr1,nr2,nr3,nr4,nr5,nr6}={r1,r2,r3,r4,r5,r6}/.{a->ax rad,b->bx rad,c->cx rad,d->dx rad};
expr=(1/{nr1,nr2,nr3,nr4,nr5,nr6}).vector;
nexpr=Chop[N[expr/.tt->1.34]];
If[nexpr==0,
paramsol=sol1+tt*dif;
(* normalize by making the first entry a monomial *)
If[dif[[1]]≠0,
paramsol=Expand[paramsol/.tt->(t-sol1[[1]])/dif[[1]]],
(* else *)
paramsol=Expand[paramsol/.tt->(t-sol1[[2]])/dif[[2]]],
];
paramsol=Chop[paramsol];
(* Normalize by removing fractions *)
lcm=LCM@@Denominator[Together[paramsol]];
If[lcm>1,paramsol=paramsol=paramsol/.t->lcm t];
If[MemberQ[paramlist,paramsol],dup=True,
dup=False;
AppendTo[paramlist,paramsol];
pcount++;
];
hres=0;
If[!dup,
Print["Parametric solution: ",paramsol ,"  relationship=",vector];
]
]
]
];
Print[pcount, " parametric solutions found."];
)

FindInverseRelationships[5]