k: (x-m)^2+(y-1)^2=16;

p: x-3*y+6=0;

solve([k,p],[x,y]);

solve(p,y);

rhs(%[1]);

subst(%,y,k);

expand(%);

%-16;

lhs(%);

ratcoef(%,x,2);

let(a,%o10);

ratcoef(%o9,x,1);

let(b,%o12);

ratcoef(%o9,x,0);

let(c,%o14);

letsimp(b^2-4*a*c);

%=0;

solve(%);

rhs(%[1]);

subst(%,m,k);

solve([%,p],[x,y]);

rhs(%o18[2]);

subst(%,m,k);

solve([%,p],[x,y]);