Plan 9 from Bell Labs’s /usr/web/sources/plan9/sys/src/libgeometry/arith3.c

 ``` #include #include #include #include /* * Routines whose names end in 3 work on points in Affine 3-space. * They ignore w in all arguments and produce w=1 in all results. * Routines whose names end in 4 work on points in Projective 3-space. */ Point3 add3(Point3 a, Point3 b){ a.x+=b.x; a.y+=b.y; a.z+=b.z; a.w=1.; return a; } Point3 sub3(Point3 a, Point3 b){ a.x-=b.x; a.y-=b.y; a.z-=b.z; a.w=1.; return a; } Point3 neg3(Point3 a){ a.x=-a.x; a.y=-a.y; a.z=-a.z; a.w=1.; return a; } Point3 div3(Point3 a, double b){ a.x/=b; a.y/=b; a.z/=b; a.w=1.; return a; } Point3 mul3(Point3 a, double b){ a.x*=b; a.y*=b; a.z*=b; a.w=1.; return a; } int eqpt3(Point3 p, Point3 q){ return p.x==q.x && p.y==q.y && p.z==q.z; } /* * Are these points closer than eps, in a relative sense */ int closept3(Point3 p, Point3 q, double eps){ return 2.*dist3(p, q)=den) return p1; return add3(p0, mul3(q, num/den)); } /* * distance from point p to segment [p0,p1] */ #define SMALL 1e-8 /* what should this value be? */ double pldist3(Point3 p, Point3 p0, Point3 p1){ Point3 d, e; double dd, de, dsq; d=sub3(p1, p0); e=sub3(p, p0); dd=dot3(d, d); de=dot3(d, e); if(dd