Language: C++

#include <cstdio>
#include <cmath>
#include <vector>
#include <algorithm>
using namespace std;

#define INF 1e10
#define EPS 1e-8

class Point3
{
public:
	double x,y,z;
	Point3()
	{}
	
	Point3(double x, double y, double z)
	:x(x),y(y),z(z)
	{}
	
	Point3 operator-(const Point3& p) const
	{
		return Point3(x-p.x, y-p.y, z-p.z);
	}
	
	Point3 operator+(const Point3& p) const
	{
		return Point3(x+p.x, y+p.y, z+p.z);
	}
	
	Point3 operator*(double s) const
	{
		return Point3(s*x, s*y, s*z);
	}
};

class Sphere
{
public:
	Point3 p;
	double r;
	
	Sphere(const Point3& p, double r)
	:p(p),r(r)
	{}
};

double dist3(Point3& a, Point3& b)
{
	return pow(a.x-b.x, 2)+pow(a.y-b.y, 2)+pow(a.z-b.z, 2);
}

double dot(Point3& a, Point3& b)
{
	return (a.x*b.x + a.y*b.y + a.z*b.z);
}

Point3 cross(Point3& a, Point3& b)
{
	Point3 r;
	r.x = (a.y*b.z - a.z*b.y);
	r.y = (a.z*b.x - a.x*b.z);
	r.z = (a.x*b.y - a.y*b.x);
	return r;
}

Point3 reflection(Point3& r, Point3& n)
{
	return (n*dot(n,r) - (r-n*dot(n,r)))*-1;
}

Point3 unit(Point3& p)
{
	double r=sqrt(p.x*p.x+p.y*p.y+p.z*p.z);
	return Point3(p.x/r, p.y/r, p.z/r);
}

vector<Point3> crosspointSpL(const Sphere& s, const Point3& l, const Point3& LV)
{
	vector<Point3> ret;
	
	double A=(LV.x*LV.x)+(LV.y*LV.y)+(LV.z*LV.z);
	double B=(LV.x*l.x)-(LV.x*s.p.x)+(LV.y*l.y)-(LV.y*s.p.y)+(LV.z*l.z)-(LV.z*s.p.z);
	double C=(l.x-s.p.x)*(l.x-s.p.x)+(l.y-s.p.y)*(l.y-s.p.y)+(l.z-s.p.z)*(l.z-s.p.z)-(s.r*s.r);
	
	double t;
	t=(-B+sqrt(B*B-A*C))/A;
	ret.push_back(Point3(l.x+t*LV.x, l.y+t*LV.y, l.z+t*LV.z));
	
	t=(-B-sqrt(B*B-A*C))/A;
	ret.push_back(Point3(l.x+t*LV.x, l.y+t*LV.y, l.z+t*LV.z));
	
	return ret;
}

bool intersectSpL(const Sphere& s, const Point3& l, const Point3& LV)
{
	double A=(LV.x*LV.x)+(LV.y*LV.y)+(LV.z*LV.z);
	double B=(LV.x*l.x)-(LV.x*s.p.x)+(LV.y*l.y)-(LV.y*s.p.y)+(LV.z*l.z)-(LV.z*s.p.z);
	double C=(l.x-s.p.x)*(l.x-s.p.x)+(l.y-s.p.y)*(l.y-s.p.y)+(l.z-s.p.z)*(l.z-s.p.z)-(s.r*s.r);

	return (B*B-A*C) > EPS;
}

bool intersectSpS(const Point3& r, const Point3& p, const Sphere& s)
{
	vector<Point3> l;
	l.push_back(r);
	l.push_back(p);
	
	double vc=dot(Point3(l[1]-l[0]), Point3(s.p-l[0]));
	if(vc+EPS<0)
	{
		return dist3(Point3(l[0]), Point3(s.p))<s.r+EPS;
	}
	
	double vv=dot(Point3(l[1]-l[0]), Point3(l[1]-l[0]));
	
	if(vc > vv)
	{
		double d=dist3(l[1], Point3(s.p));
		
		return d*d<s.r*s.r+EPS;
	}
	
	double cc=dot(Point3(s.p-l[0]), Point3(s.p-l[0]));
	
	return cc-(vc/vv)*vc < s.r*s.r+EPS;
}


int main()
{
	int N;
	while(scanf("%d", &N), N)
	{
		double u,v,w;
		scanf("%lf%lf%lf", &u,&v,&w);
		
		vector<Sphere> Sp;
		while(N--)
		{
			double x,y,z,r;
			scanf("%lf%lf%lf%lf", &x,&y,&z,&r);
			Sp.push_back(Sphere(Point3(x,y,z),r));
		}
		
		int prev=-1;
		Point3 vt(u,v,w), nt(0,0,0), at(u,v,w), ft(0,0,0);
		while(1)
		{
			int next=-1;
			double d=INF;
			for(int i=0; i<Sp.size(); i++)
			{
				if(i==prev) continue;
				if(intersectSpL(Sp[i], nt, vt))
				{
					vector<Point3> xp(crosspointSpL(Sp[i], nt, vt));
					if(dist3(xp[0], at) > dist3(xp[0], ft) && !intersectSpS(at, ft, Sp[i])) continue;
					
					double td=min(dist3(nt, xp[0]), dist3(nt,xp[1]));
					if(d>td)
					{
						next=i;
						d=td;
					}
				}
			}
			
			if(next==-1)
			{
				printf("%.3lf %.3lf %.3lf\n", nt.x, nt.y, nt.z);
				break;
			}
			
			vector<Point3> xp(crosspointSpL(Sp[next], nt, vt));
			nt=(dist3(nt, xp[0])<dist3(nt, xp[1])?xp[0]:xp[1]);
			Point3 vs=nt-Sp[next].p;
			vs=unit(vs);
			vt=reflection(vt, vs);
			prev=next;
			
			xp=crosspointSpL(Sp[next], nt, vt);
			if(dist3(nt, xp[0]) < dist3(nt, xp[1]))
			{
				at=xp[0];
				ft=xp[1];
			}
			else
			{
				at=xp[1];
				ft=xp[0];
			}

		}
	}
}

Compiler output:
a.cc: In function 'bool intersectSpS(const Point3&, const Point3&, const Sphere&)':
a.cc:111:23: error: cannot bind non-const lvalue reference of type 'Point3&' to an rvalue of type 'Point3'
  111 |         double vc=dot(Point3(l[1]-l[0]), Point3(s.p-l[0]));
      |                       ^~~~~~~~~~~~~~~~~
a.cc:53:20: note:   initializing argument 1 of 'double dot(Point3&, Point3&)'
   53 | double dot(Point3& a, Point3& b)
      |            ~~~~~~~~^
a.cc:114:30: error: cannot bind non-const lvalue reference of type 'Point3&' to an rvalue of type 'Point3'
  114 |                 return dist3(Point3(l[0]), Point3(s.p))<s.r+EPS;
      |                              ^~~~~~~~~~~~
a.cc:48:22: note:   initializing argument 1 of 'double dist3(Point3&, Point3&)'
   48 | double dist3(Point3& a, Point3& b)
      |              ~~~~~~~~^
a.cc:117:23: error: cannot bind non-const lvalue reference of type 'Point3&' to an rvalue of type 'Point3'
  117 |         double vv=dot(Point3(l[1]-l[0]), Point3(l[1]-l[0]));
      |                       ^~~~~~~~~~~~~~~~~
a.cc:53:20: note:   initializing argument 1 of 'double dot(Point3&, Point3&)'
   53 | double dot(Point3& a, Point3& b)
      |            ~~~~~~~~^
a.cc:121:38: error: cannot bind non-const lvalue reference of type 'Point3&' to an rvalue of type 'Point3'
  121 |                 double d=dist3(l[1], Point3(s.p));
      |                                      ^~~~~~~~~~~
a.cc:48:33: note:   initializing argument 2 of 'double dist3(Point3&, Point3&)'
   48 | double dist3(Point3& a, Point3& b)
      |                         ~~~~~~~~^
a.cc:126:23: error: cannot bind non-const lvalue reference of type 'Point3&' to an rvalue of type 'Point3'
  126 |         double cc=dot(Point3(s.p-l[0]), Point3(s.p-l[0]));
      |                       ^~~~~~~~~~~~~~~~
a.cc:53:20: note:   initializing argument 1 of 'double dot(Point3&, Point3&)'
   53 | double dot(Point3& a, Point3& b)
      |            ~~~~~~~~^