Your solution is right.

I suggest another variant to solve your problem without solve command by direct computation (firstly find the points of tangancy):

restart:  with(geom3d):  with(linalg):  

line(d1,[-5+2*t,-3*t+1,-13+2*t],t):  line(d2,[-7+3*t,-1-2*t,8],t):  

a:=ParallelVector(d1):  b:=ParallelVector(d2):  

c:=crossprod(a,b):

sphere(S, x^2+  y^2 + z^2 -10*x +2*y +26*z -113=0, [x,y,z]):

T:=coordinates(center(S)):  R:=radius(S): 

H:=[T + [seq(c[i]/norm(c,2)*R,i=1..3)], T - [seq(c[i]/norm(c,2)*R,i=1..3)]]:

for i to 2 do

sort(Equation(plane(P,c[1]*(x-H[i,1]) +c[2]*(y-H[i,2]) +c[3]*(z-H[i,3]) =0,[x,y,z])), {x,y,z}):

od;

I have not found a suitable way to sort, so I suggest a way to manually isolate complete squares.

The code:

restart:  with(geom3d):

a:=[-t, -t, t]: point(T, a): point(A, 1, -1, -1): point(B, 2, 1, 2): point(C, 1, 3, 1): plane(ABC, [A,B,C], [x, y, z]):

d:= distance(T, ABC):  R:= distance(T, A):

r:= R^2 - d^2;  simplify(r) assuming real;

sol:= minimize(r, location);  Sol:= op(sol[2])[1];

point(T, subs(Sol, a)):  R:= subs(Sol, R):

Eq:=lhs(Equation(sphere(S, [T,R], [x, y, z])));

(x+coeff(Eq,x)/2)^2+(y+coeff(Eq,y)/2)^2+(z+coeff(Eq,z)/2)^2=sqrt((coeff(Eq,x)^2+coeff(Eq,y)^2+coeff(Eq,z)^2)/4-tcoeff(Eq))^`2`;  # Canonical equation of S

You can use the parametric equations of the surface.

The example:

theta:=Pi/3:

plot3d( [u*cos(t)*cos(theta), u*cos(t)*sin(theta), u*sin(t)], u=0..1, t=0..2*Pi, color=red, axes=normal, scaling=constrained, style=surface, view=[-1..1, -1..1, -1.2..1.2] );

For example:

Sol:= 9.356887917*10^5 , { [{P = 330.1448618, W = 3.339008977}, 9.356887917*10^5 ] }:

P:=rhs(Sol[2,1,1,1]);  W:=rhs(Sol[2,1,1,2]);

P:=330.1448618

W:= 3.339008977

I would write your code so:

restart:  with(geom3d):

a:=[-t, -t, t]: point(T, a): point(A, 1, -1, -1): point(B, 2, 1, 2): point(C, 1, 3, 1):

plane(ABC, [A,B,C], [x, y, z]):

d:= distance(T, ABC):  R:= distance(T, A):

r:= R^2 - d^2;  simplify(r) assuming real;

sol:= minimize(r, location);  Sol:= op(sol[2])[1];

point(T, subs(Sol, a)):  R:= subs(Sol, R):

Equation(sphere(S, [T,R], [x, y, z]));

Eq:=Student[Precalculus][CompleteSquare](Equation(S));

op(1, Eq) - op(4, op(1, Eq)) = sqrt(-op(4, op(1, Eq)))^`2`;  # Canonical equation of S

Here you will find answers to your questions!

Sorry! I made a slip of the pen. Instead  2  must be  d .

After  v:=crossprod(w,a)  insert two lines into your code:

d:=igcd(v[1],v[2],v[3]):

v:=map(x->x/2,v):

The lower boundary your region is a common segment  [-1, 1] . So look for  maximum and minimum of the function of one variable  f(x,0)=x^2+x  on this segment by the commands  maximize  and  minimize . In this decision, the derivatives are not needed.

Computed plane passes through the line AB so that AB and its projection onto given plane  are perpendicular to the line of intersection of computed plane with given plane. Draw a picture and everything will be clear!

restart:

with(LinearAlgebra):

A:=<1, -2, 4>:  B:=<3, 5, -1>:  n1:=<1, 1, 1>:  M:=<x, y, z>:

n2:=CrossProduct(n1, B - A):  n:=CrossProduct(n2, B - A):  d:=igcd(seq(n[i],i=1..3)):

sort((n/d).(M - A)) = 0;   # Equation of the required plane

nops(x^3+x+1);

3

nops(expand((x+1)^100));

101

If this statement, i.e.  (abcabcabc...abc‾)  Ξ (abc‾)(mod91) , is true for any number of repeats of abc, then the number abc itself  must be divisible by 91. This follows from the fact that  abcabc - abc=abc*1000, but 1000 and 91 are coprimes.

The example:

15 mod(4);

3

None sequence is  uniquely determined by several members without additional conditions! If we look for a formula in a class of polynomials of degree <=5, then the decision will be

a[n]=CurveFitting[PolynomialInterpolation]([1, 2, 3, 4, 5, 6],[1, 3, 6, 12, 33, 51],n);

I got the same answer. In my decision  s  and  t  are the parameter values ​​for those points on the lines, for which is realized a distance equal to 1.

The code of my decision:

restart:

A:=<1, -1, 3>:   u:=<a, b, c>:   v:=<2, -5, 6>:   M:=A+s*u:   N:=t*v:

solve([a^2+b^2+c^2=1, DotProduct(<1, 1, 5>,u,conjugate = false)=0, DotProduct(M-N,u,conjugate = false)=0, DotProduct(M-N,v,conjugate = false)=0, Norm(M-N,2)=1],[a,b,c,s,t]);