Computations

> restart;

> A := [1,0];

> P:=[cos(theta),sin(theta)];

> with(Groebner):

> Id:=[(x-1)^2+y^2-6,(x-ct)^2+(y-st)^2-4,ct^2+st^2-1];

> G:=gbasis(Id,plex(x,y,ct,st));

> ylist:=[solve(G[3],y)];

> for i to 2 do xlist[i]:=solve(subs(y=ylist[i],G[4]),x): od:

> R:=[subs({ct = cos(theta),st = sin(theta)},xlist[1]),subs({ct = cos(theta),st = sin(theta)},ylist[1])]:

> S:=[subs({ct = cos(theta),st = sin(theta)},xlist[2]),subs({ct = cos(theta),st = sin(theta)},ylist[2])]:

> with(plots): with(linalg):

> OA:=scalarmul(A,t):

> LOA:=[OA[1],OA[2],t=0..1]:

> OP:=scalarmul(P,t):

> LOP:=[OP[1],OP[2],t=0..1]:

> AR:=matadd(scalarmul(A,t),scalarmul(R,1-t)):

> LAR:=[AR[1],AR[2],t=0..1]:

> AS:=matadd(scalarmul(A,t),scalarmul(S,1-t)):

> LAS:=[AS[1],AS[2],t=0..1]:

> PR:=matadd(scalarmul(P,t),scalarmul(R,1-t)):

> LPR:=[PR[1],PR[2],t=0..1]:

> PS:=matadd(scalarmul(P,t),scalarmul(S,1-t)):

> LPS:=[PS[1],PS[2],t=0..1]:

> Q:=matadd(S,R-P):

> QR:=matadd(scalarmul(Q,t),scalarmul(R,1-t)):

> LQR:=[QR[1],QR[2],t=0..1]:

> QS:=matadd(scalarmul(Q,t),scalarmul(S,1-t)):

> LQS:=[QS[1],QS[2],t=0..1]: