Untitled

/* homogeneous coordinates point */
namespace POINT2D
{
    template<class T> T _sqr(T x){ return x*x; }

    struct point
    {
        point(double _x = 0, double _y = 0, double _w = 1): x(_x), y(_y), w(_w) { normalize(); }
        double x, y, w;
        void normalize() { if (w != 0){ x /= w; y /= w; w = 1; } }
    };

    const point O(0, 0, 1);

    inline istream& operator >> (istream& istr, point& p){ p.w = 1.0; return istr >> p.x >> p.y; }

    inline ostream& operator << (ostream& ostr, const point& p){ return ostr << p.x << " " << p.y; }

    inline point operator + (const point& a, const point& b){ return point(a.x+b.x, a.y+b.y, a.w+b.w); }

    inline point operator - (const point& a, const point& b){ return point(a.x-b.x, a.y-b.y, a.w-b.w); }

    inline point operator * (double b, const point& a){ return point(b*a.x, b*a.y, b*a.w); }

    inline point operator / (const point& a, double b){ return point(a.x/b, a.y/b, a.w/b); }

    inline double operator * (const point& a, const point& b){ return a.x*b.x+a.y*b.y+a.w*b.w; }

    /* cross product operator:
       if a and b two points then c = a / b is a line crossing this points
       if a and b two lines then c = a / b is a point of intersection
       if a and b two 3d vectors then c = a / b is a vector a x b
    */
    inline point operator / (const point& a, const point& b){ return point(a.y*b.w-a.w*b.y, a.w*b.x-a.x*b.w, a.x*b.y-a.y*b.x); }

    /* square of length */
    inline double length2(const point& a){ return _sqr(a.x) + _sqr(a.y); }

    inline double length(const point& a){ return sqrt(length2(a)); }

    inline double cross_product_2d(const point& a, const point& b){ return (a/b).w; }
}