!Стереома Рубинчика

#include <cmath>
#include <iostream>
#include <cstdio>


using namespace std;


double sqr(double a) // квадрат числа
{
        return a * a;
}


bool doubleEqual(double a, double b) // сравниваем на равенство с eps
{
        return fabs(a - b) < 1e-9;
}


bool doubleLessOrEqual(double a, double b) // <= с eps
{
        return a < b || doubleEqual(a, b);
}


bool doubleLess(double a, double b) // < с eps
{
        return a < b && !doubleEqual(a, b);
}


bool doubleGreaterOrEqual(double a, double b) // >= с eps
{
        return a > b || doubleEqual(a, b);
}


bool doubleGreater(double a, double b) // > с eps
{
        return a > b && !doubleEqual(a, b);
}


double mySqrt(double a) // mySqrt с проверкой, что корректен аргумент
{
        if(doubleLess(a, 0) )
        {
                throw "sqrt(-1)";
        }
        if(a < 0) return 0;
        return sqrt(a);
}


struct Point{ // класс точки или вектора, далее мы эти понятия разделять не будем
                                // но условимся (для удобного чтения) большими буквами обозначать точки (A, B, C, D, ...)
                                // маленькими - вектора(v, u, w,...)
private: double x, y, z; // 3 поля, других не будет
public:
        Point(): x(0), y(0), z(0) {} // конструктор по умолчанию

        Point(double x, double y, double z): x(x), y(y), z(z) {} // намеренно сделаем 2 конструктора вместо Point(x = 0...)
                                                                                        //Это изавит нас от ошибок вида Point A = 2;

        void scan() // читаем координаты точки
        {
                scanf("%lf %lf %lf", &x, &y, &z);
        }

        void print() const // выводим координаты точки
        {
                printf("%.10lf %.10lf %.10lf\n", x, y, z);
        }

        Point operator+(const Point & p) const // сложение 2х векторов
        {
                return Point(x + p.x, y + p.y, z + p.z);
        }

        Point operator-(const Point & p) const // вычитание 2х векторов
        {
                return Point(x - p.x, y - p.y, z - p.z);
        }

        Point operator-() const // вычитание 2х векторов
        {
                return Point(-x, -y, -z);
        }

        Point operator*(double k) const // умножение вектора на скаляр
        {
                return Point(x * k, y * k, z * k);
        }

        Point operator/(double k) const // деление вектора на скаляр
        {
                return Point(x / k, y / k, z / k);
        }

        double operator%(const Point & p) const // скалярное произведение
        {
                return x * p.x + y * p.y + z * p.z;
        }

        Point operator*(const Point & p) const // векторное произведение
        {
                return
                        Point
                        (
                                 y * p.z - z * p.y,
                                 z * p.x - x * p.z,
                                 x * p.y - y * p.x
                        );
        }

        double length() const // длина вектора по определению из ан.гема (корень из скалярного квадрата)
        {
                return mySqrt(*this % *this);
        }

        double distTo(const Point & p) const //расстояние между 2мя точками - модуль вектора между ними
        {
                return (*this - p).length();
        }

        double distTo(const Point & A, const Point & B) const // расстояние от точки до прямой (всегда >= 0)
        {
                double d = A.distTo(B);
                if(doubleEqual(d, 0) ) // прямая должна задаваться 2мя точками
                {
                        throw "A = B";
                }
                double s = ( (*this - A) * (*this - B) ).length(); // площадь треугольника
                return s / d; // метод площадей
        }

        Point normalize(double k = 1) const // выставляет длину вектора в k
        {
                double len = length(); // Текущая длина
                if(doubleEqual(len, 0) )
                {
                        if(doubleEqual(k, 0) )
                        {
                                return Point();
                        }
                        throw "zero-size vector";
                }
                return *this * (k / len);
        }

        Point getH(const Point & A, const Point & B) const // опускаем высоту из точки на прямую (A, B)
        {
                Point C = *this;
                Point v = B - A; // направляющий вектор прямой
                Point u = C - A; // вектор, проекция которого нам нужна
                double k = v % u / v.length(); // нашли длину проекции
                v = v.normalize(k); // нашли проекцию u на v
                Point H = A + v; // подвинули точку A в проекцию конца вектора u, отложенного из A
                return H;
        }

        Point rotate(Point normal) const // поворот на 90 градусов против часовой стрелки (положительное направление)
        {
                Point v = *this;
                if(!doubleEqual(v % normal, 0) || !doubleEqual( normal.length(), 1) )
                {
                        throw "no axe";
                }
                return v * normal;
        }

        Point rotate(double alpha, const Point & normal) const // поворот на угол alpha против часовой стрелки
                                                                        // (по часовой стрелке, если alpha < 0)
        {
                return rotate(cos(alpha), sin(alpha), normal ); // делегируем задачу другому экземпляру функции
        }

        Point rotate(double cosa, double sina, const Point & normal) const // поворот с заданными косинусом и синусом
        {
                Point v = *this;
                Point u = v.rotate(normal); // вектор, перпендикулярный v, теперь (v, u) - базис, в котором мы знаем ответ
                Point w = v * cosa + u * sina; // зная координаты в базисе (v, u), нашли w
                return w;
        }

        bool isZero() const // проверка на то, что точка нулевая без сложных операций
        {
                return doubleEqual(x, 0) && doubleEqual(y, 0) && doubleEqual(z, 0);
        }

        bool isOnLine(const Point & A, const Point & B) const // точка на прямой?
        {
                return ( (A - *this) * (B - *this) ).isZero();
        }

        bool isInSegment(const Point & A, const Point & B) const // точка внутри отрезка?
        {
                return isOnLine(A, B) && doubleLessOrEqual( (A - *this) % (B - *this), 0 );
        }

        bool isInSegmentStrictly(const Point & A, const Point & B) const // точка внутри отрезка строго?
        {
                return isOnLine(A, B) && doubleLess( (A - *this) % (B - *this), 0 );
        }

        double getAngle() const
        {
                return atan2(y, x); // угол между вектором и осью ОХ
        }

        double getAngle(Point u) const
        {
                Point v = *this;
                return atan2((v * u).length() , v % u); // угол между векторами (Не ориентированный)
        }

        bool isOnPlane(const Point & A, const Point & B, const Point & C) const // проверка принадлежности точки к плоскости
        {
                return doubleEqual( (A - *this) * (B - *this) % (C - *this), 0);
        }

};


int getIntersection // ищем пересечение прямой (A, B) и прямой (C, D)
        (
                const Point & A,
                const Point & B,
                const Point & C,
                const Point & D,
                Point & O
        )
{
        if( !doubleEqual( (B - A) * (C - A) % (D - A), 0)  )
        {
                throw "It's not plane";
        }
        if( doubleEqual( ( (A - B) * (C - D) ).length(), 0) )
        {
                if(A.isOnLine(C, D) ) return 2;
                return 0;
        }
        Point normal = ( (A - B) * (C - B) ).normalize();
        Point v = B - A; // направляющий вектор прямой (A, B)
        double s1 = (C - A) * (D - A) % normal; // площадь треугольника A, C, D
        double s2 = (D - B) * (C - B) % normal; // площадь треугольника B, D, C
        double s = s1 + s2; // площадь четурёхугольника
        v = v / s;
        v = v * s1; // вектора AO и AB пропорциональны площадям s1 и s
        O = A + v; // нашли точку пересечения
        return 1; // 1 точка в пересечении
}



int getIntersection // ищем пересечение прямой (A, B) и плоскости (C, D, E)
        (
                const Point & A,
                const Point & B,
                const Point & C,
                const Point & D,
                const Point & E,
                Point & O
        )
{
        Point v = B - A; // направляющий вектор прямой (A, B)
        double V1 = (C - A) * (D - A) % (E - A); // объём тетраэдра A, C, D, E
        double V2 = (D - B) * (C - B) % (E - B); // объём тетраэдра B, D, C, E
        double V = V1 + V2; // площадь четурёхугольника
        v = v / V;
        if(doubleEqual(V, 0) )
        {
                if(A.isOnPlane(C, D, E) ) return 2;
                return 0;
        }
        v = v * V1; // вектора AO и AB пропорциональны площадям s1 и s
        O = A + v; // нашли точку пересечения
        return 1; // 1 точка в пересечении
}


bool getIntersection // ищем пересечение плоскости (A, nA) и плоскости (B, nB). Плоскость даём точкой и нормалью
        (
                const Point & A,
                const Point & nA,
                const Point & B,
                const Point & nB,
                Point & P,
                Point & Q
        )
{
        Point n = nA * nB;
        if(n.isZero() ) return false;
        Point v = n * nA;
        double k = (B - A) % nB / (v % nB);

        v = v * k;
        P = A + v;
        Q = P + n;
        return true;
}