计算几何-球和直线交点

const double eps=1e-8;
int sgn(double x) { return fabs(x) < eps ? 0 : (x > 0 ? 1 : -1); }
struct P{
    double x,y,z;
    P(double x=0,double y=0,double z=0):x(x),y(y),z(z){
    };
};
P operator-(const P &a,const P&b){
    return P(a.x-b.x,a.y-b.y,a.z-b.z);
}
P operator/(const P&a,double k){
    return P(a.x/k,a.y/k,a.z/k);
}
bool operator==(const P&a,const P&b){
    return a.x==b.x&&a.y==b.y&&a.z==b.z;
}
double dot(const P&a,const P&b){
    return a.x*b.x+a.y*b.y+a.z*b.z;
}
double dist(const P&p){
    return sqrt(p.x*p.x+p.y*p.y+p.z*p.z);
}
vector<P> c_l_intersection(const P&o,const P&s,const P&t,double r){
    vector<P>ret;
    if(s==t) return ret;
    P vec=(t-s)/dist(t-s);
    P onew=o-s;
    double dotov=dot(onew,vec);
    double delta=4*(dotov*dotov-dist(onew)*dist(onew)+r*r);
    delta=sqrt(delta);
    if(sgn(delta)<0) return ret;
    double t1=dotov+delta/2;
    ret.push_back(P(s.x+t1*vec.x,s.y+t1*vec.y,s.z+t1*vec.z));
    
    if(sgn(delta)>0){
        double t2=t1-delta;
        ret.push_back(P(s.x+t2*vec.x,s.y+t2*vec.y,s.z+t2*vec.z));
    }
    return ret;
}