Problem 163: Cross-hatched Triangles

Project Euler Problem 163 Solution

EulerSolve provides an optimized solution for Project Euler Problem 163, Cross-hatched Triangles, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Start with an equilateral triangle of side length $n$, subdivided in the cross-hatched way described in the problem, and count every triangle whose three sides lie on drawn segments. For the Project Euler instance, $n=36$. The important point is that the implementations do not classify triangles by shape or size. Instead, they turn the whole picture into a finite arrangement of supporting lines and then count which triples of those lines really bound a triangle inside the outer hull. Mathematical Approach The derivation is geometric rather than recursive. Once the cross-hatching is written as explicit families of lines, the problem becomes: among all triples of generated lines, how many have three pairwise intersections that are inside the big triangle and are not collapsed into a single point? Coordinate Model of the Outer Triangle Set $$h=\frac{\sqrt{3}}{2},\qquad A=(0,0),\qquad B=(n,0),\qquad C=\left(\frac{n}{2},nh\right).$$ This places the outer figure as a standard equilateral triangle. A point $P$ lies inside or on the boundary precisely when it is on the left side of each directed edge $AB$, $BC$, and $CA$....

Detailed mathematical approach

Problem Summary

In an equilateral triangle a straight line is drawn from each vertex to the midpoint of the opposite side; this is the building block of size $1$. A triangle of size $n$ consists of $n^2$ such blocks: lines parallel to its sides cut it into $n(n+1)/2$ small triangles pointing up and $n(n-1)/2$ pointing down, and each of them carries the three lines of the block. $T(n)$ is the number of triangles in the drawing of size $n$ whose three sides run along drawn segments, two triangles being different when they differ in shape, size, orientation or location. The statement gives $T(1)=16$ and $T(2)=104$ and asks for $T(36)$. ^[1]

The whole solution in one sentence: all drawn segments lie on $9n-3$ straight lines that run right across the big triangle, and a triangle of the drawing is the same as three of these lines that cross pairwise inside the big triangle in three different points, so the programs test every set of three lines.

Mathematical Approach

In the building block, the three sides form a triangle, but the three medians do not: they meet in one point. So a triangle of the drawing should be a set of three drawn lines that cross pairwise in three different points. We prove this rule after showing that the drawing consists of whole lines, count the building block with it and then find all points where several lines cross.

1. Coordinates

Give the small triangles side $1$, so that their height is $h=\sqrt3/2$, ^[2] and put the big triangle at $A=(0,0)$, $B=(n,0)$ and $C=(n/2,\,nh)$. The corners of the small triangles, the lattice points, are the points $(m+j/2,\,jh)$ with whole numbers $j,m\ge0$ and $j+m\le n$; for $0\le j\le n-1$, the strip between the heights $jh$ and $(j+1)h$ holds $n-j$ small triangles pointing up and $n-j-1$ pointing down. This is a piece of the triangular tiling of the plane, in which six triangles meet at every corner. ^[3] A point $(x,y)$ lies in the closed big triangle exactly when $y\ge0$, $y\le\sqrt3\,x$ and $y\le\sqrt3\,(n-x)$.

2. Six directions

In an equilateral triangle, the line from a corner to the midpoint $M$ of the opposite side cuts it into two halves with equal sides, which are congruent, so the two angles at $M$ are equal and therefore right angles. The line is thus the perpendicular bisector of that side, and it halves the $60^\circ$ angle at its corner. These three lines are the medians, and they meet at the centroid, one third of the way from each midpoint to the opposite corner. ^[4] The sides of the small triangles have the directions $0^\circ$, $60^\circ$ and $120^\circ$, so the medians have the perpendicular directions $90^\circ$, $150^\circ$ and $30^\circ$, and every drawn segment has a direction that is a multiple of $30^\circ$.

3. The drawn lines run right across

Extend a drawn segment to a whole line; we show that its whole part inside the big triangle is drawn. A line through a side consists of sides of small triangles. A median runs through a small triangle $\Delta$ from a corner to the midpoint $M$ of the opposite side. The neighbour $\Delta'$ beyond that side shares it, so the line, the perpendicular bisector of the shared side, is a median of $\Delta'$ too and runs on to the far corner $V$ of $\Delta'$. Six angles of $60^\circ$ meet at $V$, so the triangle opposite $\Delta'$ at $V$ is $\Delta'$ turned by $180^\circ$ about $V$, and the line is a median of that triangle as well. Continuing in both directions, the line stays drawn until it leaves the big triangle. So the drawing is the union of the parts inside the big triangle of finitely many lines.

4. The $9n-3$ lines

Every drawn line contains a side or a median of a small triangle and therefore a lattice point. Conversely, at a lattice point every ray in one of the six directions runs along a side of one of the small triangles there or halves one of their $60^\circ$ angles, so by Sections 2 and 3 a line in one of these directions through a lattice point is drawn if it meets the big triangle in more than one point. At $(m+j/2,\,jh)$ the expressions $y/h$, $x-y/\sqrt3$, $x+y/\sqrt3$, $x-\sqrt3\,y$, $x+\sqrt3\,y$ and $2x$ take the whole-number values $j$, $m$, $m+j$, $m-j$, $m+2j$ and $2m+j$. On the big triangle each expression lies between its values at the corners, and a line meets the triangle in more than one point exactly when its value is in this range but not an end value taken at only one corner. So the drawn lines are

the horizontal lines $y=jh$ for $0\le j\le n-1$: $n$ lines;
the lines $x-y/\sqrt3=k$ at $60^\circ$ for $0\le k\le n-1$: $n$ lines;
the lines $x+y/\sqrt3=k$ at $120^\circ$ for $1\le k\le n$: $n$ lines;
the lines $x-\sqrt3\,y=s$ at $30^\circ$ for $-(n-1)\le s\le n-1$: $2n-1$ lines;
the lines $x+\sqrt3\,y=t$ at $150^\circ$ for $1\le t\le 2n-1$: $2n-1$ lines;
the vertical lines $x=u/2$ for $1\le u\le 2n-1$: $2n-1$ lines.

Each listed value is indeed taken at a lattice point of the big triangle: for the first three families on its border, and $s$, $t$, $u$ at $m=\max(s,0)$, $j=\max(-s,0)$, at $j=\lfloor t/2\rfloor$, $m=t-2j$ and at $m=\lfloor u/2\rfloor$, $j=u-2m$. For example, $x-\sqrt3\,y=n$ is left out because it meets the big triangle only at $B$. In all there are $3n+3(2n-1)=9n-3$ lines: $6$ for $n=1$, the sides and the medians, $15$ for $n=2$ and $321$ for $n=36$.

5. Triangles are sets of three lines

The triangles of the drawing correspond one to one to the sets of three lines $\ell_1$, $\ell_2$, $\ell_3$ that cross pairwise in three different points $P_{12}$, $P_{13}$, $P_{23}$ of the big triangle. Given such lines, the three points are not on one line: $P_{12}$ and $P_{13}$ are different points of $\ell_1$, and if $P_{23}$ were on $\ell_1$, it would be a common point of $\ell_1$ and $\ell_2$, that is, $P_{12}$. The side from $P_{12}$ to $P_{13}$ lies on $\ell_1$ and, as the big triangle is convex, inside it, so it is drawn by Section 3; so are the other two sides. Conversely, a side of a triangle of the drawing is covered by finitely many drawn lines, and a line that does not contain the side meets it in at most one point, so one of the lines contains it; the three side lines cross pairwise at the corners, which lie in the big triangle. A triangle and its side lines determine each other, so every triangle is counted exactly once.

A set of three lines therefore fails exactly when two of them are parallel, when a crossing lies outside the big triangle, or when all three pass through one point. Two equal crossings already mean the last case: if $P_{12}=P_{13}$, that point lies on all three lines.

6. The building block by hand

For $n=1$ the six lines are the three sides and the three medians. Their directions are all different, so all $15$ pairs of lines cross inside the triangle, in only $7$ different points: the corners, the midpoints of the sides and the centroid $G$. There are $\binom{6}{3}=20$ ways to choose three of the six lines, and four of them go through one point: at each corner the two sides and the median from it, and at $G$ the three medians. Hence $T(1)=20-4=16$, as in the statement. In the sketch the sixteen triangles are the six small triangles around $G$; the six halves into which a median cuts the big triangle, such as $ABM$ with $M$ the midpoint of $BC$; the three triangles $ABG$, $BCG$ and $CAG$, each made of two small ones; and the big triangle itself.

7. Only three shapes

At each corner of a triangle of the drawing, the angle between its two side lines is a multiple of $30^\circ$, because the directions of both lines are. So the three angles are positive multiples of $30^\circ$ with sum $180^\circ$, that is, $2+2+2$, $1+1+4$ or $1+2+3$ steps of $30^\circ$: every triangle is equilateral, isosceles with the angles $30^\circ$, $30^\circ$ and $120^\circ$, or a right triangle with the angles $30^\circ$, $60^\circ$ and $90^\circ$, half of an equilateral one. In the building block the big triangle is equilateral, $ABG$, $BCG$ and $CAG$ are isosceles, and the six small triangles and the six halves are right triangles.

8. Size 2 in numbers

For $n=2$ three directions have $2$ lines and three have $3$. Of the $\binom{15}{3}=455$ sets of three lines, $150$ contain a parallel pair: $13$ for each direction with two lines, whose pair can be completed by any of the other $13$ lines, and $3\cdot12+1=37$ for each direction with three lines. Of the other $305$ sets, $134$ have a crossing outside the big triangle, which leaves $171$. As Section 9 shows, $67$ of them go through one point, so $T(2)=171-67=104$, as in the statement. By shape, $11$ of them are equilateral, $21$ isosceles and $72$ right triangles.

9. Where the lines cross

By Section 3, a line that enters the inside of a small triangle runs through it as one of its medians. So inside a small triangle the drawing is its three sides and three medians, and two lines can cross only at a corner, at the midpoint of a side or at a centroid. Through such a point pass

$6$ lines at a lattice point other than $A$, $B$, $C$: in each direction, the line through it meets the big triangle in more than one point (Section 4);
$3$ lines at $A$, $B$ and $C$: two sides and a median;
$3$ lines at a centroid: the medians of its small triangle;
$2$ lines at the midpoint of a side: the line of the side and the median perpendicular to it.

A point with $r$ lines lies on $\binom{r}{3}$ concurrent sets of three lines: $20$ for $r=6$, $1$ for $r=3$ and none for $r=2$. The big triangle has $(n+1)(n+2)/2$ lattice points, $n-j+1$ at the height $jh$, and $n^2$ centroids, so the number of concurrent sets is

$$\begin{aligned}&K(n)\\&\quad=20\Bigl(\tfrac{(n+1)(n+2)}{2}-3\Bigr)+3+n^2\\&\quad=10(n+1)(n+2)+n^2-57.\end{aligned}$$

It gives $K(1)=4$ and $K(2)=67$, the numbers subtracted in Sections 6 and 8. The crossing points are the vertices of an arrangement of lines, and an arrangement is called simple when at most two lines cross at each vertex; ^[5] ours is far from simple, which is why concurrent sets must be recognised.

How the Code Works

The C++ program ^[6] reads the options --size= (default $36$) and --skip-checkpoints. Unless the checkpoints are skipped, run_checkpoints requires that solve(1) and solve(2) return $16$ and $104$, and then main prints solve(size).

build_lines makes every line from two points $p$, $q$, taken from the unit triangle with corners $(0,0)$, $(1,0)$, $(1/2,h)$ and the midpoints ac and bc of two of its sides and moved into place: the Line constructor stores $a=p_y-q_y$, $b=q_x-p_x$ and $c=q_xp_y-p_xq_y$, so that $ax+by=c$ holds at both points. Its six loops emit, in this order, the horizontal lines, the $30^\circ$ lines through $(\pm i,0)$, the $60^\circ$ lines through $(i,0)$, the $120^\circ$ and $150^\circ$ lines through $(i+1,0)$ and the vertical lines $x=i/2$: the $9n-3$ lines of Section 4.

solve fills an $M\times M$ table pair_data for the pairs $i\lt j$ of the $M=9n-3$ lines. intersect solves the two line equations by Cramer's rule, ^[7]

$$\begin{aligned}x&=\frac{c_1b_2-c_2b_1}{a_1b_2-a_2b_1},\\ y&=\frac{a_1c_2-a_2c_1}{a_1b_2-a_2b_1},\end{aligned}$$

and treats the lines as parallel when $|a_1b_2-a_2b_1|\lt10^{-10}$. inside_hull evaluates line_det, the cross product $(q-p)\times(r-p)$, which is twice the signed area of the triangle $pqr$ and positive when $r$ lies to the left of the walk from $p$ to $q$. ^[8] Since $A$, $B$, $C$ run anticlockwise, a crossing is kept unless one of its three values is below $-10^{-10}$, so points on the border count as inside. The main loop runs over all $i\lt j\lt k$, enters the loop over $k$ only when the pair $(i,j)$ is valid, needs valid pairs $(i,k)$ and $(j,k)$, and drops the triple when same_point finds two crossings whose coordinates differ by at most $10^{-10}$, which by Section 5 means a concurrent set. Every other triple is a triangle, counted in a 64-bit long long.

The coordinates contain $\sqrt3$ and are rounded, with a relative error of about $10^{-16}$ per operation, ^[9] and the tolerance $10^{-10}$ leaves a wide margin. Measured for $n=36$: parallel lines give a determinant below $10^{-16}$ in absolute value, as each family is built from one direction vector, other pairs at least $\sqrt3/4\approx0.43$; every computed crossing lies within $2.1\cdot10^{-14}$ of the exact one, while different crossings are at least $\sqrt3/6\approx0.29$ apart, the distance from a centroid to a side midpoint; ^[2] border points never give a value below $-3\cdot10^{-13}$, outside points always one below $-10$. So for $n=36$ every test decides as exact arithmetic would.

The Python program has the same steps in one function solve, for size $36$ only and without options or checkpoints; a first loop that builds the horizontal lines is discarded when the list is reset, invalid pairs are stored as None, and the count, an unbounded Python integer, is returned as a string. The Java program does the same in main, also only for size $36$, with makeLine, ld and sp as constructor, line_det and same_point, crossings in an array pts of $M^2$ references beside a boolean table valid, and a long count.

Complexity and Verification

With $M=9n-3$ lines the table costs $\binom{M}{2}$ intersections and $M^2$ entries, and the main loop at most $\binom{M}{3}$ tests, so the time is $O(n^3)$ and the memory $O(n^2)$. For $n=36$, $9\,345$ of the $51\,360$ pairs are parallel and $16\,395$ cross inside the big triangle; since the loop over $k$ starts only for these, it runs $1\,631\,449$ times instead of $\binom{321}{3}=5\,461\,280$. On the tested build an entry of the C++ table takes $24$ bytes, about $2.5$ MB in all. In tests on a desktop computer the C++ program needed about a hundredth of a second, Java a few hundredths and Python about a tenth.

Only the C++ program checks $T(1)$ and $T(2)$. For this article the count was repeated exactly in the coordinates $2x$ and $y/h$, where all line equations have integer coefficients and all crossings are fractions, so no tolerance is needed: it agrees with the C++ program for every $n\le40$, including $n=36$, and Python and Java print the same number as C++. A vertex-based count, testing all sets of three crossing points for being pairwise on a common line but not on one line, agrees for every $n\le14$. The exact program also confirmed the crossing points of Section 9 for $n\le30$ and $K(n)$ for $n\le40$. On the drawing with only the $3n$ grid lines the same count gives $1,5,13,27,48,\dots$, the known numbers $\lfloor n(n+2)(2n+1)/8\rfloor$ of triangles in a triangular grid, ^[10] for every $n\le40$. For testing other programs: $T(3)=303$, $T(4)=653$, $T(5)=1196$ and $T(6)=1978$.

Footnotes and References

The references below supply the background facts used above; the description of the drawing by whole lines, the correspondence between triangles and sets of three lines and the count $K(n)$ are derived in this article.

Project Euler 163 — Cross-hatched Triangles. The statement builds the size n drawing from n² cross-hatched size 1 triangles, gives T(1) = 16 and T(2) = 104 and asks for T(36).
Wikipedia — Equilateral triangle. Side t gives the height (√3/2)t, and the inradius √3t/6 is the distance from the centroid to the sides at their midpoints.
Wikipedia — Triangular tiling. In the tiling {3,6}, six triangles meet at every vertex and fill the full 360°.
Wikipedia — Median (geometry). The three medians meet at the centroid, which divides each of them in the ratio 2 : 1.
Wikipedia — Arrangement of lines. Vertices are the points on two or more lines; an arrangement is simple when at most two lines cross at each vertex.
C++ — Euler163.cpp. The immutable revision with build_lines, solve and run_checkpoints; the implementation record for this article.
Wikipedia — Cramer's rule. For a₁x + b₁y = c₁, a₂x + b₂y = c₂ with a₁b₂ − b₁a₂ ≠ 0, the unique solution is the quotient of determinants used by intersect.
cp-algorithms — Oriented area of a triangle. The cross product gives twice the signed area: negative for a clockwise turn, positive otherwise, zero for collinear points.
Wikipedia — Machine epsilon. For IEEE 754 double precision, 2⁻⁵³ ≈ 1.11 · 10⁻¹⁶ bounds the relative rounding error.
OEIS A002717. a(n) = ⌊n(n + 2)(2n + 1)/8⌋, the number of triangles in a triangular matchstick arrangement of side n.

Problem 163 source code

C++

#include <cmath>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

using i64 = long long;

constexpr double kEps = 1e-10;

struct Options {
    int size = 36;
    bool run_checkpoints = true;
};

struct Point {
    double x = 0.0;
    double y = 0.0;
};

struct Line {
    double a = 0.0;
    double b = 0.0;
    double c = 0.0;

    Line() = default;

    Line(const Point& p, const Point& q)
        : a(p.y - q.y), b(q.x - p.x), c(q.x * p.y - p.x * q.y) {}
};

struct Intersection {
    bool valid = false;
    Point p{};
};

bool almost_equal(double x, double y) {
    return std::fabs(x - y) <= kEps;
}

bool same_point(const Point& p, const Point& q) {
    return almost_equal(p.x, q.x) && almost_equal(p.y, q.y);
}

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--size=", options.size)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.size >= 1;
}

Intersection intersect(const Line& lhs, const Line& rhs) {
    const double det = lhs.a * rhs.b - rhs.a * lhs.b;
    if (std::fabs(det) < kEps) {
        return {false, {}};
    }

    const double x = (lhs.c * rhs.b - rhs.c * lhs.b) / det;
    const double y = (lhs.a * rhs.c - rhs.a * lhs.c) / det;
    return {true, {x, y}};
}

double line_det(const Point& p, const Point& q, const Point& r) {
    return (q.x - p.x) * (r.y - p.y) - (q.y - p.y) * (r.x - p.x);
}

bool inside_hull(const Point& p, const Point& a, const Point& b, const Point& c) {
    if (line_det(a, b, p) < -kEps) {
        return false;
    }
    if (line_det(b, c, p) < -kEps) {
        return false;
    }
    if (line_det(c, a, p) < -kEps) {
        return false;
    }
    return true;
}

std::vector<Line> build_lines(const int size) {
    const double h = std::sqrt(3.0) / 2.0;

    const Point a{0.0, 0.0};
    const Point b{1.0, 0.0};
    const Point c{0.5, h};

    const Point ab{(a.x + b.x) * 0.5, (a.y + b.y) * 0.5};
    const Point ac{(a.x + c.x) * 0.5, (a.y + c.y) * 0.5};
    const Point bc{(b.x + c.x) * 0.5, (b.y + c.y) * 0.5};

    std::vector<Line> lines;
    lines.reserve(static_cast<std::size_t>(9 * size - 3));

    // Family parallel to AB.
    for (int i = 0; i < size; ++i) {
        const double y = static_cast<double>(i) * h;
        lines.emplace_back(Point{a.x, y}, Point{b.x, y});
    }

    // Family through A and BC-like directions (both sides).
    for (int i = 0; i < size; ++i) {
        const double x = static_cast<double>(i);
        lines.emplace_back(Point{x, a.y}, Point{bc.x + x, bc.y});
        if (i > 0) {
            lines.emplace_back(Point{-x, a.y}, Point{bc.x - x, bc.y});
        }
    }

    // Family parallel to AC.
    for (int i = 0; i < size; ++i) {
        const double x = static_cast<double>(i);
        lines.emplace_back(Point{x, a.y}, Point{c.x + x, c.y});
    }

    // Family parallel to BC.
    for (int i = 0; i < size; ++i) {
        const double x = static_cast<double>(i);
        lines.emplace_back(Point{x + 1.0, b.y}, Point{c.x + x, c.y});
    }

    // Family through B and AC-like directions.
    for (int i = 0; i < 2 * size - 1; ++i) {
        const double x = static_cast<double>(i);
        lines.emplace_back(Point{x + 1.0, b.y}, Point{ac.x + x, ac.y});
    }

    // Family through C and AB-like directions.
    for (int i = 1; i < 2 * size; ++i) {
        const double x = static_cast<double>(i) * c.x;
        lines.emplace_back(Point{x, 0.0}, Point{x, h});
    }

    return lines;
}

i64 solve(const int size) {
    const double h = std::sqrt(3.0) / 2.0;
    const Point a{0.0, 0.0};
    const Point b{static_cast<double>(size), 0.0};
    const Point c{0.5 * static_cast<double>(size), h * static_cast<double>(size)};

    const std::vector<Line> lines = build_lines(size);
    const int n = static_cast<int>(lines.size());

    std::vector<Intersection> pair_data(static_cast<std::size_t>(n * n));
    auto at = [&](int i, int j) -> Intersection& {
        return pair_data[static_cast<std::size_t>(i) * static_cast<std::size_t>(n) +
                         static_cast<std::size_t>(j)];
    };

    for (int i = 0; i < n; ++i) {
        for (int j = i + 1; j < n; ++j) {
            Intersection inter = intersect(lines[static_cast<std::size_t>(i)],
                                           lines[static_cast<std::size_t>(j)]);
            if (inter.valid && !inside_hull(inter.p, a, b, c)) {
                inter.valid = false;
            }
            at(i, j) = inter;
        }
    }

    i64 count = 0;
    for (int i = 0; i < n; ++i) {
        for (int j = i + 1; j < n; ++j) {
            const Intersection& ij = at(i, j);
            if (!ij.valid) {
                continue;
            }
            for (int k = j + 1; k < n; ++k) {
                const Intersection& ik = at(i, k);
                const Intersection& jk = at(j, k);
                if (!ik.valid || !jk.valid) {
                    continue;
                }
                if (same_point(ij.p, ik.p) || same_point(ij.p, jk.p) || same_point(ik.p, jk.p)) {
                    continue;
                }
                ++count;
            }
        }
    }

    return count;
}

bool run_checkpoints() {
    if (solve(1) != 16) {
        std::cerr << "Checkpoint failed for size=1" << '\n';
        return false;
    }
    if (solve(2) != 104) {
        std::cerr << "Checkpoint failed for size=2" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }
    std::cout << solve(options.size) << '\n';
    return 0;
}

Python

import math

def solve():
    SIZE = 36
    EPS = 1e-10
    h = math.sqrt(3.0) / 2.0

    a = (0.0, 0.0)
    b = (1.0, 0.0)
    c = (0.5, h)
    ab = ((a[0]+b[0])*0.5, (a[1]+b[1])*0.5)
    ac = ((a[0]+c[0])*0.5, (a[1]+c[1])*0.5)
    bc = ((b[0]+c[0])*0.5, (b[1]+c[1])*0.5)

    lines = []
    for i in range(SIZE):
        y = i * h
        lines.append((a[1]-0.0, b[0]-a[0], b[0]*a[1]-a[0]*0.0 + y*(b[0]-a[0])))
    # Rebuild properly using line equation ax+by=c from two points
    def make_line(p, q):
        a_l = p[1] - q[1]
        b_l = q[0] - p[0]
        c_l = q[0]*p[1] - p[0]*q[1]
        return (a_l, b_l, c_l)

    lines = []
    for i in range(SIZE):
        y = i * h
        lines.append(make_line((0.0, y), (1.0, y)))
    for i in range(SIZE):
        x = float(i)
        lines.append(make_line((x, 0.0), (bc[0]+x, bc[1])))
        if i > 0:
            lines.append(make_line((-x, 0.0), (bc[0]-x, bc[1])))
    for i in range(SIZE):
        x = float(i)
        lines.append(make_line((x, 0.0), (c[0]+x, c[1])))
    for i in range(SIZE):
        x = float(i)
        lines.append(make_line((x+1.0, 0.0), (c[0]+x, c[1])))
    for i in range(2*SIZE - 1):
        x = float(i)
        lines.append(make_line((x+1.0, 0.0), (ac[0]+x, ac[1])))
    for i in range(1, 2*SIZE):
        x = float(i) * c[0]
        lines.append(make_line((x, 0.0), (x, h)))

    n = len(lines)
    A = (0.0, 0.0)
    B = (float(SIZE), 0.0)
    C = (0.5*SIZE, h*SIZE)

    def intersect(l1, l2):
        det = l1[0]*l2[1] - l2[0]*l1[1]
        if abs(det) < EPS:
            return None
        x = (l1[2]*l2[1] - l2[2]*l1[1]) / det
        y = (l1[0]*l2[2] - l2[0]*l1[2]) / det
        return (x, y)

    def line_det(p, q, r):
        return (q[0]-p[0])*(r[1]-p[1]) - (q[1]-p[1])*(r[0]-p[0])

    def inside(p):
        if line_det(A, B, p) < -EPS: return False
        if line_det(B, C, p) < -EPS: return False
        if line_det(C, A, p) < -EPS: return False
        return True

    pair_data = [[None]*n for _ in range(n)]
    for i in range(n):
        for j in range(i+1, n):
            pt = intersect(lines[i], lines[j])
            if pt and inside(pt):
                pair_data[i][j] = pt

    count = 0
    for i in range(n):
        for j in range(i+1, n):
            ij = pair_data[i][j]
            if ij is None:
                continue
            for k in range(j+1, n):
                ik = pair_data[i][k]
                jk = pair_data[j][k]
                if ik is None or jk is None:
                    continue
                if abs(ij[0]-ik[0]) < EPS and abs(ij[1]-ik[1]) < EPS:
                    continue
                if abs(ij[0]-jk[0]) < EPS and abs(ij[1]-jk[1]) < EPS:
                    continue
                if abs(ik[0]-jk[0]) < EPS and abs(ik[1]-jk[1]) < EPS:
                    continue
                count += 1

    return str(count)

if __name__ == '__main__':
    print(solve())

Java

public class Euler163 {
    static final double EPS = 1e-10;

    public static void main(String[] args) {
        int SIZE = 36;
        double h = Math.sqrt(3.0) / 2.0;
        double[] cx = { 0.5 }, cy = { h };
        double[] bcx = { 0.75 }, bcy = { h / 2 };
        double[] acx = { 0.25 }, acy = { h / 2 };

        double[][] lines = new double[9 * SIZE][3];
        int lc = 0;
        for (int i = 0; i < SIZE; i++) {
            double y = i * h;
            lines[lc++] = makeLine(0, y, 1, y);
        }
        for (int i = 0; i < SIZE; i++) {
            double x = i;
            lines[lc++] = makeLine(x, 0, 0.75 + x, h / 2);
            if (i > 0)
                lines[lc++] = makeLine(-x, 0, 0.75 - x, h / 2);
        }
        for (int i = 0; i < SIZE; i++) {
            double x = i;
            lines[lc++] = makeLine(x, 0, 0.5 + x, h);
        }
        for (int i = 0; i < SIZE; i++) {
            double x = i;
            lines[lc++] = makeLine(x + 1, 0, 0.5 + x, h);
        }
        for (int i = 0; i < 2 * SIZE - 1; i++) {
            double x = i;
            lines[lc++] = makeLine(x + 1, 0, 0.25 + x, h / 2);
        }
        for (int i = 1; i < 2 * SIZE; i++) {
            double x = i * 0.5;
            lines[lc++] = makeLine(x, 0, x, h);
        }
        int n = lc;
        double AX = 0, AY = 0, BX = SIZE, BY = 0, CX = 0.5 * SIZE, CY = h * SIZE;

        double[][] pts = new double[n * n][];
        boolean[][] valid = new boolean[n][n];
        for (int i = 0; i < n; i++)
            for (int j = i + 1; j < n; j++) {
                double det = lines[i][0] * lines[j][1] - lines[j][0] * lines[i][1];
                if (Math.abs(det) < EPS)
                    continue;
                double px = (lines[i][2] * lines[j][1] - lines[j][2] * lines[i][1]) / det;
                double py = (lines[i][0] * lines[j][2] - lines[j][0] * lines[i][2]) / det;
                if (ld(AX, AY, BX, BY, px, py) < -EPS || ld(BX, BY, CX, CY, px, py) < -EPS
                        || ld(CX, CY, AX, AY, px, py) < -EPS)
                    continue;
                pts[i * n + j] = new double[] { px, py };
                valid[i][j] = true;
            }
        long count = 0;
        for (int i = 0; i < n; i++)
            for (int j = i + 1; j < n; j++) {
                if (!valid[i][j])
                    continue;
                double[] ij = pts[i * n + j];
                for (int k = j + 1; k < n; k++) {
                    if (!valid[i][k] || !valid[j][k])
                        continue;
                    double[] ik = pts[i * n + k], jk = pts[j * n + k];
                    if (sp(ij, ik) || sp(ij, jk) || sp(ik, jk))
                        continue;
                    count++;
                }
            }
        System.out.println(count);
    }

    static double[] makeLine(double px, double py, double qx, double qy) {
        return new double[] { py - qy, qx - px, qx * py - px * qy };
    }

    static double ld(double px, double py, double qx, double qy, double rx, double ry) {
        return (qx - px) * (ry - py) - (qy - py) * (rx - px);
    }

    static boolean sp(double[] a, double[] b) {
        return Math.abs(a[0] - b[0]) < EPS && Math.abs(a[1] - b[1]) < EPS;
    }
}