With an analytic surface like a sphere, it is possible to calculate the tangents of a point on the sphere with differential calculus. You have to choose first a parametric representation of the sphere S(u,v) = [ f(u,v) i,g(u,v) j,h(u,v) k ] (not unique). Then calculate the partial derivative of the parametric representation for u and v and you have your 2 tangent vectors.
But be careful, some parametric representation of a sphere (like the one given in the above link) can give tangent vector of length 0 at different place on the sphere, so choose another one.
Assuming this sphere is centered around (0,0,0), you could use also spherical coordinates (r, theta, phi) and just add a quarter turn to theta or phi to get a unit-length vector tangent to the original point.
r = sqrt(x^2 + y^2 + z^2)
theta = acos(z/r)
phi = atan2(y,x)
then add pi/2 to theta or phi
tangent.x = sin(theta)cos(phi)
tangent.y = sin(theta)sin(phi)
tangent.z = cos(theta)
The vector [0,0,1] is not a valid tangent for the point (0,0,1)
If I use theta = pi/2 and phi = pi/2 this give a second valid tangent vector [0,1,0] but now this doesn’t work for retrieving the tangents for the point (0,1,0) on the sphere.
So you need at least to check for all possibilities [(theta,phi+pi/2);(theta+pi/2,phi);(theta+pi/2),phi+pi/2)] and verify the dot product of the resulting vector with the normal of the sphere. If it is 0 then the resulting vector is a valid tangent. But maybe it exist a point on the sphere where this algorithm can give false tangents for all the combinaisons or two identics tangents.