ClawPositionMatricesO3[angles_,lengths_]:=(
Module[{P01,P02,P03,P04,P12,P23,P24,P34},
P01=RotationMatrix[angles[[1]],{0,0,1}];
P12=RotationMatrix[angles[[2]],{1,0,0}];
P23=RotationMatrix[angles[[3]],{1,0,0}];

P02=P01.P12;
P03=P02.P23;
EEPos=P01.{0,0,lengths[[1]]}+P02.{0,0,lengths[[2]]}+P03.{0,0,lengths[[3]]}
]
)

Manipulate[
ClawPositionMatricesO3[{θ1,θ2,θ3,β},{1,1,1}],
{{θ1,0},0,2Pi},
{{θ2,3Pi/2},0,2Pi},
{{θ3,Pi/2},0,2Pi},
{{β,0},0,2Pi}
]

Clear[θ1,θ2,θ3]

<<Quaternions`;

ClawPositionQuaternions[angles_,lengths_]:=(
q01=Quaternion[Cos[angles[[1]]/2],0,0,Sin[angles[[1]]/2]];
q01C=Conjugate[q01];

q12=Quaternion[Cos[angles[[2]]/2],Sin[angles[[2]]/2],0,0];
q12C=Conjugate[q12];

q23=Quaternion[Cos[angles[[3]]/2],Sin[angles[[3]]/2],0,0];
q23C=Conjugate[q23];

q34=Quaternion[Cos[angles[[4]]/2],0,0,Sin[angles[[4]]/2]];
q34C=Conjugate[q34];

l1q=Quaternion[0,0,0,lengths[[1]]];
l2q=Quaternion[0,0,0,lengths[[2]]];
l3q=Quaternion[0,0,0,lengths[[3]]];
q02=q01**q12;
q02C=Conjugate[q02];
q03=q02**q23;
q03C=Conjugate[q03];
q04=q03**q34;
q04=Conjugate[q04];
N4Pos=q01**l1q**q01C+q02**l2q**q02C+q03**l3q**q03C
)

Manipulate[
ClawPositionQuaternions[{θ1,θ2,θ3,β},{1,1,1}],
{{θ1,0},0,2Pi},
{{θ2,3Pi/2},0,2Pi},
{{θ3,Pi/2},0,2Pi},
{{β,0},0,2Pi}
]

Clear[θ1,θ2,θ3]

ClawPositionMatricesSE3[angles_,lengths_]:=(
Module[{R01,R12,R23},
{θ1,θ2,θ3}={angles[[1]],angles[[2]],angles[[3]]};
R01={{Cos[θ1],-Sin[θ1],0,0},{Sin[θ1],Cos[θ1],0,0},{0,0,1,0},{0,0,0,1}};
R12={{1,0,0,0},{0,Cos[θ2],-Sin[θ2],0},{0,Sin[θ2],Cos[θ2],lengths[[1]]},{0,0,0,1}};
R23={{1,0,0,0},{0,Cos[θ3],-Sin[θ3],0},{0,Sin[θ3],Cos[θ3],lengths[[2]]},{0,0,0,1}};

EEPos=R01.R12.R23.{0,0,lengths[[3]],1}
]
)

Manipulate[
ClawPositionMatricesSE3[{θ1,θ2,θ3},{1,1,1}],
{{θ1,0},0,2Pi},
{{θ2,3Pi/2},0,2Pi},
{{θ3,Pi/2},0,2Pi}
]

Clear[θ1,θ2,θ3]

DQP[Q1_,Q2_]:=(
Return[{
{Q1[[1,1]]**Q2[[1,1]],Q1[[1,1]]**Q2[[1,2]]+Q1[[1,2]]**Q2[[1,1]]},
{Quaternion[0,0,0,0],Q1[[1,1]]**Q2[[1,1]]}
}]
)

ClawPositionDualQuaternions[angles_,lengths_]:=(
Module[{q01M,q01CM,q12M,q12CM,q23M,q23CM,l1,l2,l3,ql3,θ1,θ2,θ3},
{θ1,θ2,θ3}={angles[[1]],angles[[2]],angles[[3]]};
{l1,l2,l3}={lengths[[1]]/2,lengths[[2]]/2,lengths[[3]]};
q01M={{Quaternion[Cos[((θ1)/2)],0,0,Sin[((θ1)/2)]],Quaternion[0,0,0,0]},{Quaternion[0,0,0,0],Quaternion[Cos[((θ1)/2)],0,0,Sin[((θ1)/2)]]}};
q01CM=Conjugate[q01M];
q12M={{Quaternion[Cos[((θ2)/2)],Sin[((θ2)/2)],0,0],Quaternion[0,0,l1 Sin[((θ2)/2)],l1 Cos[((θ2)/2)]]},{Quaternion[0,0,0,0],Quaternion[Cos[((θ2)/2)],Sin[((θ2)/2)],0,0]}};
q12CM={{Quaternion[Cos[((θ2)/2)],-Sin[((θ2)/2)],0,0],-Quaternion[0,0,-l1 Sin[((θ2)/2)],-l1 Cos[((θ2)/2)]]},{Quaternion[0,0,0,0],Quaternion[Cos[((θ2)/2)],-Sin[((θ2)/2)],0,0]}};
q23M={{Quaternion[Cos[((θ3)/2)],Sin[((θ3)/2)],0,0],Quaternion[0,0,l2 Sin[((θ3)/2)],l2 Cos[((θ3)/2)]]},{Quaternion[0,0,0,0],Quaternion[Cos[((θ3)/2)],Sin[((θ3)/2)],0,0]}};
q23CM={{Quaternion[Cos[((θ3)/2)],-Sin[((θ3)/2)],0,0],-Quaternion[0,0,-l2 Sin[((θ3)/2)],-l2 Cos[((θ3)/2)]]},{Quaternion[0,0,0,0],Quaternion[Cos[((θ3)/2)],-Sin[((θ3)/2)],0,0]}};
ql3={{Quaternion[1,0,0,0],Quaternion[0,0,0,l3]},{Quaternion[0,0,0,0],Quaternion[1,0,0,0]}};
Return[Simplify[DQP[DQP[DQP[DQP[DQP[DQP[q01M,q12M],q23M],ql3],q23CM],q12CM],q01CM]][[1,2]]];
]
)

Manipulate[
ClawPositionDualQuaternions[{θ1,θ2,θ3},{1,1,1}],
{{θ1,0},0,2Pi},
{{θ2,3Pi/2},0,2Pi},
{{θ3,Pi/2},0,2Pi}
]

Manipulate[
claw1=Cuboid[{0.2,-1,0},{-0.2,1,0.4}];
claw2=Cuboid[{0.2,-1,0},{-0.2,-0.8,0.8}];
claw3=Cuboid[{0.2,1,0},{-0.2,0.8,0.8}];
claw={claw1,claw2,claw3};
clawR=Rotate[claw,0,{0,0,1}];
clawRT=Translate[clawR,{0,0,3}];
unit4={clawRT};
unit4T=Translate[unit4,{0,0,l3*10-3}];
segment4=Cylinder[{{0,0,0},{0,0,l3*10}},1/4];
unit4TR=Rotate[{unit4T,segment4},θ3,{1,0,0}];
pivot3=Cylinder[{{-1/2,0,0},{1/2,0,0}},1/2];
unit3={unit4TR,pivot3};
Segment2=Cylinder[{{0,0,0},{0,0,l2*10}},1/4];
unit3T=Translate[unit3,{0,0,l2*10}];
unit2={unit3T,Segment2};
unit2R=Rotate[unit2,θ2,{1,0,0}];
pivot2=Cylinder[{{-1/2,0,0},{1/2,0,0}},1/2];
unit2RT=Translate[{unit2R,pivot2},{0,0,l1*10}];
Segment1=Cylinder[{{0,0,0},{0,0,l1*10}},1/4];
unit1={unit2RT,Segment1};
pivot1=Cylinder[{{0,0,0},{0,0,3}},1/2];
unit1R=Rotate[{unit1,pivot1},θ1,{0,0,1}];
base=Cylinder[{{0,0,0},{0,0,-1/2}},2];
Graphics3D[{unit1R,base,Text[Style[ClawPositionMatricesO3[{θ1,θ2,θ3},{l1,l2,l3}],White],{0,0,-3},FormatType->StandardForm]},Background->RGBColor[0,0.2,0.3]
,Boxed->False,ImageSize->{350,350}],{{θ1,0},0,2Pi},{{θ2,3Pi/2},0,2Pi},{{θ3,Pi/2},0,2Pi},{{l1,1},0.4,2},{{l2,1},0.4,2},{{l3,1},0.4,2},TrackedSymbols->{θ3,θ2,θ1,l1,l2,l3},
ControlPlacement->Left]

InverseKinematicsExponentialMap[FinalPosition_,lengths_]:=(
Module[{a,b,x,y,z,s1,s2,s3,s4,c1,c2,c3,l1,l2,l3,v1,v2,v3,r1,r2,r3,p0,p,angles,f},
{l1,l2,l3}={lengths[[1]],lengths[[2]],lengths[[3]]};
{x,y,z}={FinalPosition[[1]],FinalPosition[[2]],FinalPosition[[3]]};
If[(x^2+y^2+(z-l1)^2>(l1+l2)^2||x^2+y^2+(z-l1)^2<(l2-l1)^2)||z<0,Print["Desired position is out of configuration space."];Return[{0,0,0}],
If[{x,y,z}=={0,0,l1},Return[{{0,Pi/2,Pi},{0,Pi/2,3Pi},{0,3Pi/2,Pi},{0,3Pi/2,3Pi}}];,
If[{x,y,z}=={0,0,z},
b=l2;c=l3;
If[z>l1,
a=z-l1;
Return[{
{0,2Pi-ArcCos[(-c^2+a^2+b^2)/(2a b)],Pi-ArcCos[(-a^2+b^2+c^2)/(2b c)]},{0,ArcCos[(-c^2+a^2+b^2)/(2a b)],Pi+ArcCos[(-a^2+b^2+c^2)/(2b c)]}}]];
If[z<l1,
a=l1-z;
Return[{
{0,Pi+ArcCos[(-c^2+a^2+b^2)/(2a b)],Pi+ArcCos[(-a^2+b^2+c^2)/(2b c)]},{0,Pi-ArcCos[(-c^2+a^2+b^2)/(2a b)],Pi-ArcCos[(-a^2+b^2+c^2)/(2b c)]}}]];,
v1={0,0,1};v2={1,0,0};v3={1,0,0};r1={0,0,0};r2={0,0,l1};r3={0,0,l1+l2};p0={0,0,l1+l2+l3};p={x,y,z};
a=p0+Cross[Sin[θ3]v3,p0-r3]+Cross[(1-Cos[θ3])v3,Cross[v3,p0-r3]];
b=p-Cross[Sin[θ1]v1,p-r1]+Cross[(1-Cos[θ1])v1,Cross[v1,p-r1]];
s1=Simplify[a.v2-b.v2];
s3=Simplify[(a-r2).(a-r2)-(b-r2).(b-r2)];
s2=Simplify[1-((b-a).Cross[v2,Cross[v2,(a-r2)]])/(Cross[v2,Cross[v2,(a-r2)]].Cross[v2,Cross[v2,(a-r2)]])-Cos[θ2]];
s4=Simplify[((b-a).(Cross[v2,a-r2]))/(Cross[v2,a-r2].Cross[v2,a-r2])-Sin[θ2]];angles=List[ToRules[Reduce[s1==0&&s3==0&&s2==0&&s4==0,{θ1,θ2,θ3}]/.{C[1]->1,C[2]->1,C[3]->1}]];
f=Dimensions[angles][[1]];
Return[Table[{angles[[i,1,2]],angles[[i,2,2]],angles[[i,3,2]]},{i,f}]];
];];];]
)

Panel[
DynamicModule[
{x=0,y=1,z=2},
Column[{
InputField[Dynamic[{x,y,z}]],
Dynamic[InverseKinematicsExponentialMap[{x,y,z},{1,1,1}]]
}]]]

LM[polys_,variablesOrder_]:=(
dim=Dimensions[polys][[1]];
LMSet={};
For[i=1,i≤dim,i++,
LMCoef=CoefficientRules[polys[[i]],variablesOrder][[1,1]];
LMSet=Append[LMSet,FromCoefficientRules[{LMCoef->1},variablesOrder]];
];
Return[LMSet];
)

LC[polys_,variablesOrder_]:=(
dim=Dimensions[polys][[1]];
LCSet={};
For[i=1,i≤dim,i++,
LCCoef=CoefficientRules[polys[[i]],variablesOrder][[1,2]];
LCSet=Append[LCSet,LCCoef];
];
Return[LCSet];
)

ReducedGroebnerBasis[GB_,variablesOrder_]:=(
Module[{F,f0},
F=GB;
H={};
n=Dimensions[F][[1]];
Table[
f0=F[[1]];
F=Delete[F,1];
If[
PolynomialReduce[LM[{f0},variablesOrder][[1]],LM[F,variablesOrder],variablesOrder][[2]]=!=0&&
PolynomialReduce[LM[{f0},variablesOrder][[1]],LM[H,variablesOrder],variablesOrder][[2]]=!=0,
f0=PolynomialReduce[f0,H,variablesOrder][[2]];
f0=Expand[f0/LC[{f0},variablesOrder]][[1]];
f0=Simplify[f0];
H=Append[H,f0]
],{i,1,n}];
Return[H]
]
)

polys={s1(l2 s2+l3(s2 c3+c2 s3))-x,-c1(l2 s2+l3(s2 c3+c2 s3))-y,l1+l2 c2+l3(c2 c3-s2 s3)-z,c1^2+s1^2-1,c2^2+s2^2-1,c3^2+s3^2-1}

GB=GroebnerBasis[polys,{c3,s3,c2,s2,c1,s1}];

ReducedGroebnerBasis[GB,{c3,s3,c2,s2,c1,s1}]

InverseKinematicsGroebnerBasis[FinalPosition_,lengths_]:=(
Module[{a,b,c,x,y,z,s1,s2,s3,c1,c2,c3,l1,l2,l3,angles,f},
{l1,l2,l3}={lengths[[1]],lengths[[2]],lengths[[3]]};
{x,y,z}={FinalPosition[[1]],FinalPosition[[2]],FinalPosition[[3]]};
If[(x^2+y^2+(z-l1)^2>(l1+l2)^2||x^2+y^2+(z-l1)^2<(l2-l1)^2)||z<0,Print["Desired position is out of configuration space."];Return[{0,0,0}],
If[{x,y,z}=={0,0,l1},Return[{{0,Pi/2,Pi},{0,Pi/2,3Pi},{0,3Pi/2,Pi},{0,3Pi/2,3Pi}}];,
If[{x,y,z}=={0,0,z},
b=l2;c=l3;
If[z>l1,
a=z-l1;
Return[{
{0,2Pi-ArcCos[(-c^2+a^2+b^2)/(2a b)],Pi-ArcCos[(-a^2+b^2+c^2)/(2b c)]},{0,ArcCos[(-c^2+a^2+b^2)/(2a b)],Pi+ArcCos[(-a^2+b^2+c^2)/(2b c)]}}]];
If[z<l1,
a=l1-z;
Return[{
{0,Pi+ArcCos[(-c^2+a^2+b^2)/(2a b)],Pi+ArcCos[(-a^2+b^2+c^2)/(2b c)]},{0,Pi-ArcCos[(-c^2+a^2+b^2)/(2a b)],Pi-ArcCos[(-a^2+b^2+c^2)/(2b c)]}}]];,
If[{x,y,z}=={0,y,l1},
a=Abs[y];b=l2;c=l3;
If[y>0,Return[{
{0,3Pi/2+ArcCos[(-c^2+a^2+b^2)/(2a b)],Pi+ArcCos[(-a^2+b^2+c^2)/(2c b)]},{0,3Pi/2-ArcCos[(-c^2+a^2+b^2)/(2a b)],Pi-ArcCos[(-a^2+b^2+c^2)/(2c b)]},{Pi,Pi/2-ArcCos[(-c^2+a^2+b^2)/(2a b)],Pi-ArcCos[(-a^2+b^2+c^2)/(2c b)]},{Pi,Pi/2+ArcCos[(-c^2+a^2+b^2)/(2a b)],Pi+ArcCos[(-a^2+b^2+c^2)/(2c b)]}}]];
If[y<0,Return[{
{Pi,3Pi/2+ArcCos[(-c^2+a^2+b^2)/(2a b)],Pi+ArcCos[(-a^2+b^2+c^2)/(2c b)]},{Pi,3Pi/2-ArcCos[(-c^2+a^2+b^2)/(2a b)],Pi-ArcCos[(-a^2+b^2+c^2)/(2c b)]},{0,Pi/2-ArcCos[(-c^2+a^2+b^2)/(2a b)],Pi-ArcCos[(-a^2+b^2+c^2)/(2c b)]},{0,Pi/2+ArcCos[(-c^2+a^2+b^2)/(2a b)],Pi+ArcCos[(-a^2+b^2+c^2)/(2c b)]}}]];,
If[{x,y,z}=={0,y,z},
sol=Solve[s1==0&&-1+(c1^2)==0&&((l1^4)+(l2^4)+4 c1 (l2^3) s2 y-4 (l1^3) z+4 c1 l2 s2 y (-(l3^2)+(y^2)+(z^2))+(((-(l3^2)+(y^2)+(z^2)))^2)-4 l1 z (-(l3^2)+(l2^2) (-1+2 (s2^2))+2 c1 l2 s2 y+(y^2)+(z^2))+2 (l1^2) (-(l3^2)+(l2^2) (-1+2 (s2^2))+2 c1 l2 s2 y+(y^2)+3 (z^2))-2 (l2^2) ((l3^2)-(1+2 (s2^2)) (y^2)+(1-2 (s2^2)) (z^2)))/(4 (l2^2) ((l1^2)+(y^2)-2 l1 z+(z^2)))==0&&(1/(2 l1 l2-2 l2 z))((l1^2)+2 c2 l1 l2+(l2^2)-(l3^2)+2 c1 l2 s2 y+(y^2)-2 l1 z-2 c2 l2 z+(z^2))==0&&-(1/(2 l2 l3 (l1-z)))((l1^2) (2 l2 s2+c1 y)-2 l1 (l2 l3 s3+2 l2 s2 z+c1 y z)+c1 y ((l2^2)-(l3^2)+(y^2)+(z^2))+2 l2 (l3 s3 z+s2 ((y^2)+(z^2))))==0&&-(((l1^2)-(l2^2)-2 c3 l2 l3-(l3^2)+(y^2)-2 l1 z+(z^2))/(2 l2 l3))==0,{s1,s2,s3,c1,c2,c3}];
solution={s1,s2,s3,c1,c2,c3}/.sol;
f=Dimensions[solution][[1]];
Return[Table[angles=Reduce[Sin[θ1]==solution[[i,1]]&&Sin[θ2]==solution[[i,2]]&&Sin[θ3]==solution[[i,3]]&&Cos[θ1]==solution[[i,4]]&&Cos[θ2]==solution[[i,5]]&&Cos[θ3]==solution[[i,6]],{θ1,θ2,θ3},Reals]/.C[1]->1;
{angles[[1,2]],angles[[2,2]],angles[[3,2]]},{i,f}]],
If[{x,y,z}=={x,y,l1},
sol=Solve[(((-1+(s1^2)) (x^2)+(s1^2) (y^2))/((x^2)+(y^2)))==0&&c1+((s1 y)/x)==0&&-(((l2^2) s1-2 l2 s2 x+s1 (-(l3^2)+(x^2)+(y^2)))/(2 l2 x))==0&&((l2^4)+(((-(l3^2)+(x^2)+(y^2)))^2)-2 (l2^2) ((l3^2)-(-1+2 (c2^2)) ((x^2)+(y^2))))/(4 (l2^2) ((x^2)+(y^2)))==0&&((l3 s3 x-c2 s1 ((x^2)+(y^2)))/(l3 x))==0&&(((l2^2)+2 c3 l2 l3+(l3^2)-(x^2)-(y^2))/(2 l2 l3))==0,{s1,s2,s3,c1,c2,c3}];
solution={s1,s2,s3,c1,c2,c3}/.sol;
f=Dimensions[solution][[1]];
Return[Table[angles=Reduce[Sin[θ1]==solution[[i,1]]&&Sin[θ2]==solution[[i,2]]&&Sin[θ3]==solution[[i,3]]&&Cos[θ1]==solution[[i,4]]&&Cos[θ2]==solution[[i,5]]&&Cos[θ3]==solution[[i,6]],{θ1,θ2,θ3},Reals]/.C[1]->1;
{angles[[1,2]],angles[[2,2]],angles[[3,2]]},{i,f}]],
sol=Solve[(((-1+(s1^2)) (x^2)+(s1^2) (y^2))/((x^2)+(y^2)))==0&&c1+((s1 y)/x)==0&&((l1^4) x+(l2^4) x-4 (l2^3) s1 s2 ((x^2)+(y^2))-4 (l1^3) x z-4 l2 s1 s2 ((x^2)+(y^2)) (-(l3^2)+(x^2)+(y^2)+(z^2))+x (((-(l3^2)+(x^2)+(y^2)+(z^2)))^2)+2 (l2^2) x (-(l3^2)+(x^2)+2 (s2^2) (x^2)+(y^2)+2 (s2^2) (y^2)-(z^2)+2 (s2^2) (z^2))-4 l1 z ((l2^2) (-1+2 (s2^2)) x-2 l2 s1 s2 ((x^2)+(y^2))+x (-(l3^2)+(x^2)+(y^2)+(z^2)))+2 (l1^2) ((l2^2) (-1+2 (s2^2)) x-2 l2 s1 s2 ((x^2)+(y^2))+x (-(l3^2)+(x^2)+(y^2)+3 (z^2))))/(4 (l2^2) x ((l1^2)+(x^2)+(y^2)-2 l1 z+(z^2)))==0&&(1/(2 l2 x (l1-z)))((l1^2) x+(l2^2) x+2 l1 x (c2 l2-z)-2 l2 (s1 s2 ((x^2)+(y^2))+c2 x z)+x (-(l3^2)+(x^2)+(y^2)+(z^2)))==0&&(1/(2 l2 l3 x (l1-z)))((l2^2) s1 ((x^2)+(y^2))+(l1^2) (-2 l2 s2 x+s1 ((x^2)+(y^2)))+s1 ((x^2)+(y^2)) (-(l3^2)+(x^2)+(y^2)+(z^2))+2 l1 (-s1 ((x^2)+(y^2)) z+l2 x (l3 s3+2 s2 z))-2 l2 x (l3 s3 z+s2 ((x^2)+(y^2)+(z^2))))==0&&-(((l1^2)-(l2^2)-2 c3 l2 l3-(l3^2)+(x^2)+(y^2)-2 l1 z+(z^2))/(2 l2 l3))==0,{s1,s2,s3,c1,c2,c3}];
solution={s1,s2,s3,c1,c2,c3}/.sol;
f=Dimensions[solution][[1]];
Return[Table[angles=Reduce[Sin[θ1]==solution[[i,1]]&&Sin[θ2]==solution[[i,2]]&&Sin[θ3]==solution[[i,3]]&&Cos[θ1]==solution[[i,4]]&&Cos[θ2]==solution[[i,5]]&&Cos[θ3]==solution[[i,6]],{θ1,θ2,θ3},Reals]/.C[1]->1;
{angles[[1,2]],angles[[2,2]],angles[[3,2]]},{i,f}]]];
];];];];];
]
)

Panel[
DynamicModule[
{x=0,y=1,z=2},
Column[{
InputField[Dynamic[{x,y,z}]],
Dynamic[InverseKinematicsGroebnerBasis[{x,y,z},{1,1,1}]]
}]]]

InverseKinematicsCircles[FinalPosition_,lengths_]:=(
Module[{a,b,c,x,y,z,u,v,d,tmp,l1,l2,l3,α,γ},
{l1,l2,l3}={lengths[[1]],lengths[[2]],lengths[[3]]};
{x,y,z}={FinalPosition[[1]],FinalPosition[[2]],FinalPosition[[3]]};
If[(x^2+y^2+(z-l1)^2>(l1+l2)^2||x^2+y^2+(z-l1)^2<(l2-l1)^2)(*||z<0*),Print["Desired position is out of configuration space."];Return[{0,0,0}],
If[{x,y,z}=={0,0,l1},Return[{{0,Pi/2,Pi},{0,Pi/2,3Pi},{0,3Pi/2,Pi},{0,3Pi/2,3Pi}}];,
If[{x,y,z}=={0,0,z},
b=l2;c=l3;
If[z>l1,
a=z-l1;
Return[{
{0,2Pi-ArcCos[(-c^2+a^2+b^2)/(2a b)],Pi-ArcCos[(-a^2+b^2+c^2)/(2b c)]},{0,ArcCos[(-c^2+a^2+b^2)/(2a b)],Pi+ArcCos[(-a^2+b^2+c^2)/(2b c)]}}]];
If[z<l1,
a=l1-z;
Return[{
{0,Pi+ArcCos[(-c^2+a^2+b^2)/(2a b)],Pi+ArcCos[(-a^2+b^2+c^2)/(2b c)]},{0,Pi-ArcCos[(-c^2+a^2+b^2)/(2a b)],Pi-ArcCos[(-a^2+b^2+c^2)/(2b c)]}}]];,
tmp={x,y,z}-{0,0,l1};
d=Norm[tmp];
v=(l2^2-l3^2+d^2)/(2d);
u=d-v;
α=ArcCos[tmp[[3]]/d];
γ=ArcCos[-y/Sqrt[x^2+y^2]];
If[x>0,
Return[{
{γ,α-ArcCos[v/l2],Pi-ArcSin[u/l3]-ArcSin[v/l2]},
{γ,α+ArcCos[v/l2],Pi+ArcSin[u/l3]+ArcSin[v/l3]},
{γ+Pi,2Pi-α+ArcCos[v/l2],Pi+ArcSin[u/l3]+ArcSin[v/l3]},
{γ+Pi,2Pi-α-ArcCos[v/l2],Pi-ArcSin[u/l3]-ArcSin[v/l2]}
}],
Return[{
{2Pi-γ,α-ArcCos[v/l2],Pi-ArcSin[u/l3]-ArcSin[v/l2]},
{2Pi-γ,α+ArcCos[v/l2],Pi+ArcSin[u/l3]+ArcSin[v/l3]},
{Pi-γ,2Pi-α+ArcCos[v/l2],Pi+ArcSin[u/l3]+ArcSin[v/l3]},
{Pi-γ,2Pi-α-ArcCos[v/l2],Pi-ArcSin[u/l3]-ArcSin[v/l2]}
}]
];];];];]
)

Panel[
DynamicModule[
{x=0,y=1,z=2},
Column[{
InputField[Dynamic[{x,y,z}]],
Dynamic[InverseKinematicsCircles[{x,y,z},{1,1,1}]]
}]]]

RobotModel[angles_,lengths_]:=(
Module[{θ1,θ2,θ3,β,l1,l2,l3},
{l1,l2,l3}={lengths[[1]],lengths[[2]],lengths[[3]]};
{θ1,θ2,θ3}={angles[[1]],angles[[2]],angles[[3]]};
claw1=Cuboid[{0.2,-1,0},{-0.2,1,0.4}];
claw2=Cuboid[{0.2,-1,0},{-0.2,-0.8,0.8}];
claw3=Cuboid[{0.2,1,0},{-0.2,0.8,0.8}];
claw={claw1,claw2,claw3};
clawR=Rotate[claw,0,{0,0,1}];
clawRT=Translate[clawR,{0,0,3}];
unit4={clawRT};
unit4T=Translate[unit4,{0,0,l3*10-3}];
segment4=Cylinder[{{0,0,0},{0,0,l3*10}},1/4];
unit4TR=Rotate[{unit4T,segment4},θ3,{1,0,0}];
pivot3=Cylinder[{{-1/2,0,0},{1/2,0,0}},1/2];
unit3={unit4TR,pivot3};
Segment2=Cylinder[{{0,0,0},{0,0,l2*10}},1/4];
unit3T=Translate[unit3,{0,0,l2*10}];
unit2={unit3T,Segment2};
unit2R=Rotate[unit2,θ2,{1,0,0}];
pivot2=Cylinder[{{-1/2,0,0},{1/2,0,0}},1/2];
unit2RT=Translate[{unit2R,pivot2},{0,0,l1*10}];
Segment1=Cylinder[{{0,0,0},{0,0,l1*10}},1/4];
unit1={unit2RT,Segment1};
pivot1=Cylinder[{{0,0,0},{0,0,3}},1/2];
unit1R=Rotate[{unit1,pivot1},θ1,{0,0,1}];
base=Cylinder[{{0,0,0},{0,0,-1/2}},2];
Graphics3D[{unit1R,base,Text[Style[ClawPositionMatricesO3[{θ1,θ2,θ3,0},{l1,l2,l3}],12,White],{0,0,-3},FormatType->StandardForm],Text[Style[{θ1 180/Pi,θ2 180/Pi,θ3 180/Pi},12,White]
,{0,0,-5},FormatType->StandardForm]},Background->RGBColor[0,0.2,0.3]
,Boxed->False,ImageSize->{350,350}]])

finalConf={0,0,0};

Manipulate[
{t1,t2,t3}=finalConf;
sol=InverseKinematicsCircles[Position,{1,1,1}];
n=Dimensions[sol][[1]];
For[i=1,i<n+1,i++,
If[sol[[i,1]]≥2Pi,sol[[i,1]]=sol[[i,1]]-Floor[sol[[i,1]]/(2Pi)]2Pi];
If[sol[[i,2]]≥2Pi,sol[[i,2]]=sol[[i,2]]-Floor[sol[[i,2]]/(2Pi)]2Pi];
If[sol[[i,3]]≥2Pi,sol[[i,3]]=sol[[i,3]]-Floor[sol[[i,3]]/(2Pi)]2Pi];
];
finalConf=sol[[1]];
For[i=2,i<n+1,i++,
If[Norm[Subtract[{t1,t2,t3},finalConf]]>Norm[Subtract[{t1,t2,t3},sol[[i]]]],finalConf=sol[[i]]];
];
Animate[
RobotModel[{θ1,θ2,θ3},{1,1,1}],
{θ1,t1,finalConf[[1]]},
{θ2,t2,finalConf[[2]]},
{θ3,t3,finalConf[[3]]},
AnimationRepetitions->1
],
{Position,{List[1,1,1]}},ControlType->InputField,TrackedSymbols:>{Position}
]

tC=Timing[InverseKinematicsCircles[{1,1,1},{1,1,1}]];

tGB=Timing[InverseKinematicsGroebnerBasis[{1,1,1},{1,1,1}]];

tEM=Timing[InverseKinematicsExponentialMap[{1,1,1},{1,1,1}]];

BarChart[{tC[[1]],tGB[[1]],tEM[[1]]}]