#include <iostream>
#include <cmath>

#ifndef SIDE
#define SIDE 300
#endif

using T = int;

struct Point {
    double x, y, z;

    constexpr double mag2() const {
        return x*x + y*y + z*z;
    }

    constexpr double dot(Point rhs) const {
        return x*rhs.x + y*rhs.y + z*rhs.z;
    }

    constexpr Point cross(Point rhs) const {
        return Point{y*rhs.z - z*rhs.y, z*rhs.x - x*rhs.z, x*rhs.y - y*rhs.x};
    }

    friend constexpr Point operator+(Point p, Point q) {
        return Point{p.x + q.x, p.y + q.y, p.z + q.z};
    }

    friend constexpr Point operator-(Point p, Point q) {
        return Point{p.x - q.x, p.y - q.y, p.z - q.z};
    }

    friend constexpr bool operator==(Point a, Point b) {
        return a.x == b.x && a.y == b.y && a.z == b.z;
    }
};

static_assert(Point{1,0,0}.cross(Point{0,1,0}) == Point{0,0,1});

// These two points define the centerline of the hole.
#ifndef PETER
// Traditional solution, with the hole centered on the space diagonal.
static constexpr Point X1 = {0, 0, 0};
static constexpr Point X2 = {1, 1, 1};
#else
// Pieter Nieuwland's solution.
static constexpr Point X1 = {0.5, 0.5, 0.5};
static constexpr Point X2 = X1 + Point{6, 6, 3};
#endif

inline double sqr(double x) { return x * x; }

// http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html
inline double squared_distance_to_hole_centerline(Point X0)
{
    Point X1_X0 = (X1 - X0);
    static constexpr Point X2_X1 = (X2 - X1);
    double numerator = (X1_X0.mag2() * X2_X1.mag2()) - sqr(X1_X0.dot(X2_X1));
    static constexpr double denominator = X2_X1.mag2();
    return numerator / denominator;
}

inline double distance_to_vertical_plane(Point X0)
{
    // N is the direction of the vertical plane's normal vector.
    static constexpr Point N = {X2.y - X1.y, X1.x - X2.x, 0};
    static constexpr double a = N.x;
    static constexpr double b = N.y;
    static constexpr double c = N.z;
    static constexpr double d = -(a*X1.x + b*X1.y + c*X1.z);
    // The equation of the vertical plane is ax+by+cz+d = 0.
    static const double denominator = sqrt(a*a + b*b + c*c);
    double numerator = fabs(a*X0.x + b*X0.y + c*X0.z + d);
    return numerator / denominator;
}

inline double distance_to_horizontal_plane(Point X0)
{
    // dir1 is the direction of the vertical plane's normal vector.
    static constexpr Point dir1 = {X2.y - X1.y, X1.x - X2.x, 0};
    // dir2 is the direction of the hole's center line.
    static constexpr Point dir2 = {X2.x - X1.x, X2.y - X1.y, X2.z - X1.z};
    // N is the direction of the horizontal plane's normal vector.
    static constexpr Point N = dir1.cross(dir2);
    static constexpr double a = N.x;
    static constexpr double b = N.y;
    static constexpr double c = N.z;
    static constexpr double d = -(a*X1.x + b*X1.y + c*X1.z);
    // The equation of the vertical plane is ax+by+cz+d = 0.
    static const double denominator = sqrt(a*a + b*b + c*c);
    double numerator = fabs(a*X0.x + b*X0.y + c*X0.z + d);
    return numerator / denominator;
}

double max_d1 = -1;
double max_d2 = -1;

inline double is_in_hole(double x, double y, double z)
{
    Point p = {x, y, z};
    double d1 = distance_to_vertical_plane(p);
    double d2 = distance_to_horizontal_plane(p);
    assert(d1 >= 0);
    assert(d2 >= 0);
    max_d1 = std::max(d1, max_d1);
    max_d2 = std::max(d2, max_d2);
    static const double half_side_of_hole = 0.5;
    if (d1 < half_side_of_hole && d2 < half_side_of_hole) {
        return 1.0;
    } else if (d1 <= half_side_of_hole && d2 <= half_side_of_hole) {
        return 0.5;
    } else {
        return 0.0;
    }
}

int main()
{
    double hole_volume = 0.0;
    for (T x = 0; x < SIDE; ++x) {
      double xd = double(x) / SIDE;
      for (T y = 0; y < SIDE; ++y) {
        double yd = double(y) / SIDE;
        for (T z = 0; z < SIDE; ++z) {
          double zd = double(z) / SIDE;
          hole_volume += is_in_hole(xd, yd, zd);
        }
      }
    }
    std::cout << "max_d1 = " << max_d1 << "\n";
    std::cout << "max_d2 = " << max_d2 << "\n";
    T total = SIDE * SIDE * SIDE;
    std::cout << hole_volume << " of " << total << " points were inside the hole.\n";
    std::cout << (total - hole_volume) << " of " << total << " points remain.\n";
    std::cout << "The volume of the remaining cube is about " << (double(total - hole_volume) / double(total)) << ".\n";
}