Stacking Pythagorean Triangles

Question

Background

A Pythagorean triangle is a right triangle where each side length is an integer (that is, the side lengths form a Pythagorean triple):

Using the sides of this triangle, we can attach two more non-congruent Pythagorean triangles as follows:

We can continue with this pattern as we see fit, so long as no two triangles overlap, and connecting sides have equal length:

The question is, how many non-congruent Pythagorean triangles can we fit in a given space?

The Input

You will receive two integers as input, W and H, through function arguments, STDIN, strings, or whatever you like. The integers may be received as decimal, hexadecimal, binary, unary (good luck, Retina), or any other integer base. You may assume that max(W, H) <= 2^15 - 1.

The Output

Your program or function should calculate a list of non-overlapping connected non-congruent Pythagorean triangles and output a list of sets of three coordinates each, where the coordinates in a set form one of the Pythagorean triangles when connected by lines. Coordinates must be real numbers in our space (x must be in the interval [0, W] and y must be in the interval [0, H]), and distance should be accurate to machine precision. The order of the triangles and the exact format of each coordinate is not important.

It must be possible to "walk" from one triangle to any other only stepping over connected boundaries.

Using the above diagram as an example, let our input be W = 60, H = 60.

Our output could then be the following list of coordinates:

(0, 15), (0, 21), (8, 15)
(0, 21), (14.4, 40.2), (8, 15)
(0, 15), (8, 0), (8, 15)
(8, 0), (8, 15), (28, 15)
(8, 15), (28, 15), (28, 36)
(28, 15), (28, 36), (56, 36)

Now, 6 triangles is most certainly not the best we can do given our space. Can you do better?

Rules and Scoring

• Your score for this challenge is the number of triangles your program generates given the input of W = 1000, H = 1000. I reserve the right to change this input if I suspect someone hardcoding this case.

• You may not use builtins that calculate Pythagorean triples, and you may not hardcode a list of more than 2 Pythagorean triples (if you hardcode your program to always start with a (3, 4, 5) triangle, or some similar starting circumstance, that is okay).

• You may write your submission in any language. Readability and commenting are encouraged.

• You may find a list of Pythagorean triples here.

• Standard Loopholes are disallowed.

Answer 1

Python 3, 109

This was certainly a deceptively hard challenge. There were many times writing the code that I found myself questioning my basic geometry abilities. That being said, I'm pretty happy with the result. I put no effort in to coming up with a complex algorithm for placing the triangles, and instead my code just blunders through always placing the smallest it can find. I hope I can improve this later on, or my answer will spurn others to find a better algorithm! All in all, a very fun problem, and it produces some interesting pictures.

Here's the code:

import time
import math

W = int(input("Enter W: "))
H = int(input("Enter H: "))

middle_x = math.floor(W/2)
middle_y = math.floor(H/2)

sides = [ # each side is in the format [length, [x0, y0], [x1, y1]]
    [3,[middle_x,middle_y],[middle_x+3,middle_y]],
    [4,[middle_x,middle_y],[middle_x,middle_y+4]],
    [5,[middle_x+3,middle_y],[middle_x,middle_y+4]]
    ]

triangles = [[0,1,2]] # each triangle is in the format [a, b, c] where a, b and c are the indexes of sides

used_triangles = [[3,4,5]] # a list of used Pythagorean triples, where lengths are ordered (a < b < c)

max_bounds_length = math.sqrt(W**2 + H**2)

def check_if_pyth_triple(a,b): # accepts two lists of the form [l, [x0,y0], [x1,y1]] defining two line segments
    # returns 0 if there are no triples, 1 if there is a triple with a right angle on a,
    # and 2 if there is a triple with the right angle opposite a
    c = math.sqrt(a[0]**2 + b[0]**2)
    if c.is_integer():
        if not sorted([a[0], b[0], c]) in used_triangles:
            return 1
        return 0
    else:
        if a[0] > b[0]:
            c = math.sqrt(a[0]**2 - b[0]**2)
            if c.is_integer() and not sorted([a[0], b[0], c]) in used_triangles:
                return 2
        return 0

def check_if_out_of_bounds(p):
    out = False
    if p[0] < 0 or p[0] > W:
        out = True
    if p[1] < 0 or p[1] > H:
        out = True
    return out

def in_between(a,b,c):
    maxi = max(a,c)
    mini = min(a,c)
    return mini < b < maxi

def sides_intersect(AB,CD): # accepts two lists of the form [l, [x0,y0], [x1,y1]] defining two line segments
    # doesn't count overlapping lines
    A = AB[1]
    B = AB[2]
    C = CD[1]
    D = CD[2]

    if A[0] == B[0]: # AB is vertical
        if C[0] == D[0]: # CD is vertical
            return False
        else:
            m1 = (C[1] - D[1])/(C[0] - D[0]) # slope of CD
            y = m1*(A[0] - C[0]) + C[1] # the y value of CD at AB's x value
            return in_between(A[1], y, B[1]) and in_between(C[0], A[0], D[0])
    else:
        m0 = (A[1] - B[1])/(A[0] - B[0]) # slope of AB
        if C[0] == D[0]: # CD is vertical
            y = m0*(C[0] - A[0]) + A[1] # the y value of CD at AB's x value
            return in_between(C[1], y, D[1]) and in_between(A[0],C[0],B[0])
        else:
            m1 = (C[1] - D[1])/(C[0] - D[0]) # slope of CD
            if m0 == m1:
                return False
            else:
                x = (m0*A[0] - m1*C[0] - A[1] + C[1])/(m0 - m1)
                return in_between(A[0], x, B[0]) and in_between(C[0], x, D[0])

def check_all_sides(b,triangle):
    no_intersections = True
    for side in sides:
        if sides_intersect(side, b):
            no_intersections = False
            break

    return no_intersections

def check_point_still_has_room(A): # This function is needed for the weird case when all 2pi degrees
    # around a point are filled by triangles, but you could fit in a small triangle into another one
    # already built around the point. Doing this won't cause sides_intersect() to detect it because
    # the sides will all be parallel. Crazy stuff.
    connecting_sides = []
    for side in sides:
        if A in side:
            connecting_sides.append(side)

    match_count = 0
    slopes = []
    for side in connecting_sides:
        B = side[1]
        if A == B:
            B = side[2]
        if not A[0] == B[0]:
            slope = round((A[1]-B[1])/(A[0]-B[0]),4)
        else:
            if A[1] < B[1]:
                slope = "infinity"
            else:
                slope = "neg_infinity"
        if slope in slopes:
            match_count -= 1
        else:
            slopes.append(slope)
            match_count += 1

    return match_count != 0

def construct_b(a,b,pyth_triple_info,straight_b_direction,bent_b_direction):
    # this function finds the correct third point of the triangle given a and the length of b
    # pyth_triple_info determines if a is a leg or the hypotenuse
    # the b_directions determine on which side of a the triangle should be formed
    a_p = 2 # this is the index of the point in a that is not the shared point with b
    if a[1] != b[1]:
        a_p = 1

    vx = a[a_p][0] - b[1][0] # v is our vector, and these are the coordinates, adjusted so that
    vy = a[a_p][1] - b[1][1] # the shared point is the origin

    if pyth_triple_info == 1:
        # because the dot product of orthogonal vectors is zero, we can use that and the Pythagorean formula
        # to get this simple formula for generating the coordinates of b's second point
        if vy == 0:
            x = 0
            y = b[0]
        else:
            x = b[0]/math.sqrt(1+((-vx/vy)**2)) # b[0] is the desired length
            y = -vx*x/vy

        x = x*straight_b_direction # since the vector is orthagonal, if we want to reverse the direction,
        y = y*straight_b_direction # it just means finding the mirror point

    elif pyth_triple_info == 2: # this finds the intersection of the two circles of radii b[0] and c 
        # around a's endpoints, which is the third point of the triangle if a is the hypotenuse
        c = math.sqrt(a[0]**2 - b[0]**2)
        D = a[0]
        A = (b[0]**2 - c**2 + D**2 ) / (2*D)
        h = math.sqrt(b[0]**2 - A**2)
        x2 = vx*(A/D)
        y2 = vy*(A/D)        
        x = x2 + h*vy/D
        y = y2 - h*vx/D

        if bent_b_direction == -1: # this constitutes reflection of the vector (-x,-y) around the normal vector n,
            # which accounts for finding the triangle on the opposite side of a
            dx = -x
            dy = -y
            v_length = math.sqrt(vx**2 + vy**2)
            nx = vx/v_length
            ny = vy/v_length

            d_dot_n = dx*nx + dy*ny

            x = dx - 2*d_dot_n*nx
            y = dy - 2*d_dot_n*ny

    x = x + b[1][0] # adjust back to the original frame
    y = y + b[1][1]

    return [x,y]

def construct_triangle(side_index):
    a = sides[side_index] # a is the base of the triangle
    a_p = 1
    b = [1, a[a_p], []] # side b, c is hypotenuse

    for index, triangle in enumerate(triangles):
        if side_index in triangle:
            triangle_index = index
            break

    triangle = list(triangles[triangle_index])
    triangle.remove(side_index)

    add_tri = False

    straight_b = construct_b(a,b,1,1,1)

    bent_b = construct_b(a,b,2,1,1)

    A = sides[triangle[0]][1]
    if A in a:
        A = sides[triangle[0]][2]

    Ax = A[0] - b[1][0] # adjusting A so that it's a vector
    Ay = A[1] - b[1][1]

    # these are for determining if construct_b() is going to the correct side
    triangle_on_side = (a[2][0]-a[1][0])*(A[1]-a[1][1]) - (a[2][1]-a[1][1])*(A[0]-a[1][0])
    straight_b_on_side = (a[2][0]-a[1][0])*(straight_b[1]-a[1][1]) - (a[2][1]-a[1][1])*(straight_b[0]-a[1][0])
    bent_b_on_side = (a[2][0]-a[1][0])*(bent_b[1]-a[1][1]) - (a[2][1]-a[1][1])*(bent_b[0]-a[1][0])

    straight_b_direction = 1
    if (triangle_on_side > 0 and straight_b_on_side > 0) or (triangle_on_side < 0 and straight_b_on_side < 0):
        straight_b_direction = -1

    bent_b_direction = 1
    if (triangle_on_side > 0 and bent_b_on_side > 0) or (triangle_on_side < 0 and bent_b_on_side < 0):
        bent_b_direction = -1


    a_ps = []
    for x in [1,2]:
        if check_point_still_has_room(a[x]): # here we check for that weird exception
            a_ps.append(x)

    while True:
        out_of_bounds = False
        if b[0] > max_bounds_length:
            break

        pyth_triple_info = check_if_pyth_triple(a,b)

        for a_p in a_ps:
            if a_p == 1: # this accounts for the change in direction when switching a's points
                new_bent_b_direction = bent_b_direction
            else:
                new_bent_b_direction = -bent_b_direction

            b[1] = a[a_p]
            if pyth_triple_info > 0:
                b[2] = construct_b(a,b,pyth_triple_info,straight_b_direction,new_bent_b_direction)

                if check_if_out_of_bounds(b[2]): # here is the check to make sure we don't go out of bounds
                    out_of_bounds = True
                    break

                if check_all_sides(b,triangle):
                    if pyth_triple_info == 1:
                        c = [math.sqrt(a[0]**2 + b[0]**2), a[3-a_p], b[2]]
                    else:
                        c = [math.sqrt(a[0]**2 - b[0]**2), a[3-a_p], b[2]]

                    if check_all_sides(c,triangle):
                        add_tri = True
                        break

        if out_of_bounds or add_tri:
            break

        b[0] += 1 # increment the length of b every time the loop goes through

    if add_tri: # this adds a new triangle
        sides.append(b)
        sides.append(c)
        sides_len = len(sides)
        triangles.append([side_index, sides_len - 2, sides_len - 1])
        used_triangles.append(sorted([a[0], b[0], c[0]])) # so we don't use the same triangle again

def build_all_triangles(): # this iterates through every side to see if a new triangle can be constructed
    # this is probably where real optimization would take place so more optimal triangles are placed first
    t0 = time.clock()

    index = 0
    while index < len(sides):
        construct_triangle(index)
        index += 1

    t1 = time.clock()

    triangles_points = [] # this is all for printing points
    for triangle in triangles:
        point_list = []
        for x in [1,2]:
            for side_index in triangle:
                point = sides[side_index][x]
                if not point in point_list:
                    point_list.append(point)
        triangles_points.append(point_list)

    for triangle in triangles_points:
        print(triangle)

    print(len(triangles), "triangles placed in", round(t1-t0,3), "seconds.")

def matplotlib_graph(): # this displays the triangles with matplotlib
    import pylab as pl
    import matplotlib.pyplot as plt
    from matplotlib import collections as mc

    lines = []
    for side in sides:
        lines.append([side[1],side[2]])

    lc = mc.LineCollection(lines)
    fig, ax = pl.subplots()
    ax.add_collection(lc)
    ax.autoscale()
    ax.margins(0.1)
    plt.show()

build_all_triangles()

Here's a graph of the output for W = 1000 and H = 1000 with 109 triangles:

Here is W = 10000 and H = 10000 with 724 triangles:

Call the matplotlib_graph() function after build_all_triangles() to graph the triangles.

I think the code runs reasonably fast: at W = 1000 and H = 1000 it takes 0.66 seconds, and at W = 10000 and H = 10000 it takes 45 seconds using my crappy laptop.

Answer 2

C++, 146 triangles (part 1/2)

Result as Image

Algorithm Description

This uses a breadth-first search of the solution space. In each step, it starts with all unique configurations of k triangles that fit in the box, and builds all unique configurations of k + 1 triangles by enumerating all options of adding an unused triangle to any of the configurations.

The algorithm is basically set up to find the absolute maximum with an exhaustive BFS. And it does that successfully for smaller sizes. For example, for a box of 50x50, it finds the maximum in about 1 minute. But for 1000x1000, the solution space is way too big. To allow it to terminate, I trim the list of solutions after each step. The number of solutions that is kept is given by a command line argument. For the solution above, a value of 50 was used. This resulted in a runtime of about 10 minutes.

The outline of the main steps looks like this:

1. Generate all Pythagorean triangles that could potentially fit inside the box.
2. Generate the initial solution set consisting of solutions with 1 triangle each.
3. Loop over generations (triangle count).
   1. Eliminate invalid solutions from the solution set. These are solutions that either do not fit inside the box, or have overlap.
   2. If solution set is empty, we are done. The solution set from the previous generation contains the maxima.
   3. Trim solution set to given size if trim option was enabled.
   4. Loop over all solutions in current generation.
      1. Loop over all sides in perimeter of solution.
         1. Find all triangles that have a side length matching the perimeter side, and that are not in the solution yet.
         2. Generate the new solutions resulting from adding the triangles, and add the solutions to the solution set of the new generation.
4. Print solutions.

One critical aspect in the whole scheme is that configurations will generally be generated multiple times, and we are only interested in unique configurations. So we need a unique key that defines a solution, which must be independent of the order of the triangles used while generating the solution. For example, using coordinates for the key would not work at all, since they can be completely different if we arrived at the same solution in multiple orders. What I used is the set of triangle indices in the global list, plus a set of "connector" objects that define how the triangles are connected. So the key only encodes the topology, independent of the construction order and position in 2D space.

While more an implementation aspect, another part that is not entirely trivial is deciding if and how the whole thing fits into the given box. If you really want to push the boundaries, it is obviously necessary to allow for rotation to fit inside the box.

I'll try and add some comments to the code in part 2 later, in case somebody wants to dive into the details of how this all works.

Result in Official Text Format

(322.085, 641.587) (318.105, 641.979) (321.791, 638.602)
(318.105, 641.979) (309.998, 633.131) (321.791, 638.602)
(318.105, 641.979) (303.362, 639.211) (309.998, 633.131)
(318.105, 641.979) (301.886, 647.073) (303.362, 639.211)
(301.886, 647.073) (297.465, 638.103) (303.362, 639.211)
(301.886, 647.073) (280.358, 657.682) (297.465, 638.103)
(301.886, 647.073) (283.452, 663.961) (280.358, 657.682)
(301.886, 647.073) (298.195, 666.730) (283.452, 663.961)
(301.886, 647.073) (308.959, 661.425) (298.195, 666.730)
(301.886, 647.073) (335.868, 648.164) (308.959, 661.425)
(335.868, 648.164) (325.012, 669.568) (308.959, 661.425)
(308.959, 661.425) (313.666, 698.124) (298.195, 666.730)
(313.666, 698.124) (293.027, 694.249) (298.195, 666.730)
(313.666, 698.124) (289.336, 713.905) (293.027, 694.249)
(298.195, 666.730) (276.808, 699.343) (283.452, 663.961)
(335.868, 648.164) (353.550, 684.043) (325.012, 669.568)
(303.362, 639.211) (276.341, 609.717) (309.998, 633.131)
(276.808, 699.343) (250.272, 694.360) (283.452, 663.961)
(335.868, 648.164) (362.778, 634.902) (353.550, 684.043)
(362.778, 634.902) (367.483, 682.671) (353.550, 684.043)
(250.272, 694.360) (234.060, 676.664) (283.452, 663.961)
(362.778, 634.902) (382.682, 632.942) (367.483, 682.671)
(382.682, 632.942) (419.979, 644.341) (367.483, 682.671)
(419.979, 644.341) (379.809, 692.873) (367.483, 682.671)
(353.550, 684.043) (326.409, 737.553) (325.012, 669.568)
(353.550, 684.043) (361.864, 731.318) (326.409, 737.553)
(353.550, 684.043) (416.033, 721.791) (361.864, 731.318)
(416.033, 721.791) (385.938, 753.889) (361.864, 731.318)
(385.938, 753.889) (323.561, 772.170) (361.864, 731.318)
(385.938, 753.889) (383.201, 778.739) (323.561, 772.170)
(383.201, 778.739) (381.996, 789.673) (323.561, 772.170)
(323.561, 772.170) (292.922, 743.443) (361.864, 731.318)
(323.561, 772.170) (296.202, 801.350) (292.922, 743.443)
(250.272, 694.360) (182.446, 723.951) (234.060, 676.664)
(335.868, 648.164) (330.951, 570.319) (362.778, 634.902)
(330.951, 570.319) (381.615, 625.619) (362.778, 634.902)
(330.951, 570.319) (375.734, 565.908) (381.615, 625.619)
(330.951, 570.319) (372.989, 538.043) (375.734, 565.908)
(323.561, 772.170) (350.914, 852.648) (296.202, 801.350)
(323.561, 772.170) (362.438, 846.632) (350.914, 852.648)
(234.060, 676.664) (217.123, 610.807) (283.452, 663.961)
(217.123, 610.807) (249.415, 594.893) (283.452, 663.961)
(375.734, 565.908) (438.431, 559.733) (381.615, 625.619)
(382.682, 632.942) (443.362, 567.835) (419.979, 644.341)
(443.362, 567.835) (471.667, 606.601) (419.979, 644.341)
(323.561, 772.170) (393.464, 830.433) (362.438, 846.632)
(372.989, 538.043) (471.272, 556.499) (375.734, 565.908)
(372.989, 538.043) (444.749, 502.679) (471.272, 556.499)
(372.989, 538.043) (365.033, 521.897) (444.749, 502.679)
(443.362, 567.835) (544.353, 553.528) (471.667, 606.601)
(544.353, 553.528) (523.309, 622.384) (471.667, 606.601)
(544.353, 553.528) (606.515, 572.527) (523.309, 622.384)
(419.979, 644.341) (484.688, 697.901) (379.809, 692.873)
(444.749, 502.679) (552.898, 516.272) (471.272, 556.499)
(217.123, 610.807) (170.708, 516.623) (249.415, 594.893)
(484.688, 697.901) (482.006, 753.837) (379.809, 692.873)
(484.688, 697.901) (571.903, 758.147) (482.006, 753.837)
(419.979, 644.341) (535.698, 636.273) (484.688, 697.901)
(276.808, 699.343) (228.126, 812.299) (250.272, 694.360)
(228.126, 812.299) (185.689, 726.188) (250.272, 694.360)
(228.126, 812.299) (192.246, 829.981) (185.689, 726.188)
(393.464, 830.433) (449.003, 936.807) (362.438, 846.632)
(393.464, 830.433) (468.505, 926.625) (449.003, 936.807)
(416.033, 721.791) (471.289, 833.915) (385.938, 753.889)
(471.289, 833.915) (430.252, 852.379) (385.938, 753.889)
(350.914, 852.648) (227.804, 874.300) (296.202, 801.350)
(192.246, 829.981) (114.401, 834.898) (185.689, 726.188)
(114.401, 834.898) (155.433, 715.767) (185.689, 726.188)
(217.123, 610.807) (91.773, 555.523) (170.708, 516.623)
(91.773, 555.523) (141.533, 457.421) (170.708, 516.623)
(141.533, 457.421) (241.996, 407.912) (170.708, 516.623)
(141.533, 457.421) (235.365, 394.457) (241.996, 407.912)
(241.996, 407.912) (219.849, 525.851) (170.708, 516.623)
(241.996, 407.912) (304.896, 419.724) (219.849, 525.851)
(91.773, 555.523) (55.917, 413.995) (141.533, 457.421)
(571.903, 758.147) (476.260, 873.699) (482.006, 753.837)
(571.903, 758.147) (514.819, 890.349) (476.260, 873.699)
(571.903, 758.147) (587.510, 764.886) (514.819, 890.349)
(587.510, 764.886) (537.290, 898.778) (514.819, 890.349)
(587.510, 764.886) (592.254, 896.801) (537.290, 898.778)
(587.510, 764.886) (672.455, 761.831) (592.254, 896.801)
(55.917, 413.995) (113.819, 299.840) (141.533, 457.421)
(113.819, 299.840) (149.275, 293.604) (141.533, 457.421)
(544.353, 553.528) (652.112, 423.339) (606.515, 572.527)
(652.112, 423.339) (698.333, 461.597) (606.515, 572.527)
(535.698, 636.273) (651.250, 731.917) (484.688, 697.901)
(651.250, 731.917) (642.213, 756.296) (484.688, 697.901)
(304.896, 419.724) (299.444, 589.636) (219.849, 525.851)
(304.896, 419.724) (369.108, 452.294) (299.444, 589.636)
(304.896, 419.724) (365.965, 299.326) (369.108, 452.294)
(304.896, 419.724) (269.090, 347.067) (365.965, 299.326)
(114.401, 834.898) (0.942, 795.820) (155.433, 715.767)
(114.401, 834.898) (75.649, 947.412) (0.942, 795.820)
(192.246, 829.981) (124.489, 994.580) (114.401, 834.898)
(269.090, 347.067) (205.435, 217.901) (365.965, 299.326)
(205.435, 217.901) (214.030, 200.956) (365.965, 299.326)
(182.446, 723.951) (68.958, 600.078) (234.060, 676.664)
(182.446, 723.951) (32.828, 633.179) (68.958, 600.078)
(652.112, 423.339) (763.695, 288.528) (698.333, 461.597)
(763.695, 288.528) (808.220, 324.117) (698.333, 461.597)
(763.695, 288.528) (811.147, 229.162) (808.220, 324.117)
(652.112, 423.339) (627.572, 321.247) (763.695, 288.528)
(627.572, 321.247) (660.872, 244.129) (763.695, 288.528)
(652.112, 423.339) (530.342, 344.618) (627.572, 321.247)
(652.112, 423.339) (570.488, 453.449) (530.342, 344.618)
(627.572, 321.247) (503.633, 267.730) (660.872, 244.129)
(365.965, 299.326) (473.086, 450.157) (369.108, 452.294)
(365.965, 299.326) (506.922, 344.440) (473.086, 450.157)
(365.965, 299.326) (394.633, 260.827) (506.922, 344.440)
(394.633, 260.827) (537.381, 303.535) (506.922, 344.440)
(811.147, 229.162) (979.067, 234.338) (808.220, 324.117)
(698.333, 461.597) (706.660, 655.418) (606.515, 572.527)
(811.147, 229.162) (982.117, 135.385) (979.067, 234.338)
(982.117, 135.385) (999.058, 234.954) (979.067, 234.338)
(365.965, 299.326) (214.375, 186.448) (394.633, 260.827)
(811.147, 229.162) (803.145, 154.590) (982.117, 135.385)
(803.145, 154.590) (978.596, 102.573) (982.117, 135.385)
(214.375, 186.448) (314.969, 126.701) (394.633, 260.827)
(314.969, 126.701) (508.984, 192.909) (394.633, 260.827)
(314.969, 126.701) (338.497, 88.341) (508.984, 192.909)
(338.497, 88.341) (523.725, 138.884) (508.984, 192.909)
(338.497, 88.341) (359.556, 11.163) (523.725, 138.884)
(808.220, 324.117) (801.442, 544.012) (698.333, 461.597)
(801.442, 544.012) (739.631, 621.345) (698.333, 461.597)
(660.872, 244.129) (732.227, 78.877) (763.695, 288.528)
(660.872, 244.129) (644.092, 40.821) (732.227, 78.877)
(808.220, 324.117) (822.432, 544.659) (801.442, 544.012)
(660.872, 244.129) (559.380, 47.812) (644.092, 40.821)
(660.872, 244.129) (556.880, 242.796) (559.380, 47.812)
(556.880, 242.796) (528.882, 242.437) (559.380, 47.812)
(808.220, 324.117) (924.831, 449.189) (822.432, 544.659)
(924.831, 449.189) (922.677, 652.177) (822.432, 544.659)
(922.677, 652.177) (779.319, 785.836) (822.432, 544.659)
(779.319, 785.836) (696.630, 771.054) (822.432, 544.659)
(779.319, 785.836) (746.412, 969.918) (696.630, 771.054)
(779.319, 785.836) (848.467, 840.265) (746.412, 969.918)
(848.467, 840.265) (889.327, 872.428) (746.412, 969.918)
(746.412, 969.918) (619.097, 866.541) (696.630, 771.054)
(779.319, 785.836) (993.200, 656.395) (848.467, 840.265)
(993.200, 656.395) (935.157, 864.450) (848.467, 840.265)
(993.200, 656.395) (995.840, 881.379) (935.157, 864.450)
(338.497, 88.341) (34.607, 5.420) (359.556, 11.163)
(338.497, 88.341) (189.294, 204.357) (34.607, 5.420)
(189.294, 204.357) (158.507, 228.296) (34.607, 5.420)
(158.507, 228.296) (38.525, 230.386) (34.607, 5.420)
(158.507, 228.296) (41.694, 412.358) (38.525, 230.386)

Code

See part 2 for the code. This was broken into 2 parts to work around post size limits.

The code is also available on PasteBin.

Answer 3

C++, 146 triangles (part 2/2)

Continued from part 1. This was broken into 2 parts to work around post size limits.

Code

Comments to be added.

#include <cmath>
#include <vector>
#include <set>
#include <map>
#include <sstream>
#include <iostream>

class Vec2 {
public:
    Vec2()
      : m_x(0.0f), m_y(0.0f) {
    }

    Vec2(float x, float y)
      : m_x(x), m_y(y) {
    }

    float x() const {
        return m_x;
    }

    float y() const {
        return m_y;
    }

    void normalize() {
        float s = 1.0f / sqrt(m_x * m_x + m_y * m_y);
        m_x *= s;
        m_y *= s;
    }

    Vec2 operator+(const Vec2& rhs) const {
        return Vec2(m_x + rhs.m_x, m_y + rhs.m_y);
    }

    Vec2 operator-(const Vec2& rhs) const {
        return Vec2(m_x - rhs.m_x, m_y - rhs.m_y);
    }

    Vec2 operator*(float s) const {
        return Vec2(m_x * s, m_y * s);
    }

private:
    float m_x, m_y;
};

static float cross(const Vec2& v1, const Vec2& v2) {
    return v1.x() * v2.y() - v1.y() * v2.x();
}

class Triangle {
public:
    Triangle()
      : m_sideLenA(0), m_sideLenB(0), m_sideLenC(0) {
    }

    Triangle(int sideLenA, int sideLenB, int sideLenC)
      : m_sideLenA(sideLenA),
        m_sideLenB(sideLenB),
        m_sideLenC(sideLenC) {
    }

    int getSideLenA() const {
        return m_sideLenA;
    }

    int getSideLenB() const {
        return m_sideLenB;
    }

    int getSideLenC() const {
        return m_sideLenC;
    }

private:
    int m_sideLenA, m_sideLenB, m_sideLenC;
};

class Connector {
public:
    Connector(int sideLen, int triIdx1, int triIdx2, bool flipped);

    bool operator<(const Connector& rhs) const;

    void print() const {
        std::cout << m_sideLen << "/" << m_triIdx1 << "/"
                  << m_triIdx2 << "/" << m_flipped << " ";
    }

private:
    int m_sideLen;
    int m_triIdx1, m_triIdx2;
    bool m_flipped;
};

typedef std::vector<Triangle> TriangleVec;
typedef std::multimap<int, int> SideMap;

typedef std::set<int> TriangleSet;
typedef std::set<Connector> ConnectorSet;

class SolutionKey {
public:
    SolutionKey() {
    }

    void init(int triIdx);
    void add(int triIdx, const Connector& conn);

    bool containsTriangle(int triIdx) const;
    int minTriangle() const;

    bool operator<(const SolutionKey& rhs) const;

    void print() const;

private:
    TriangleSet m_tris;
    ConnectorSet m_conns;
};

typedef std::map<SolutionKey, class SolutionData> SolutionMap;

class SolutionData {
public:
    SolutionData()
      : m_lastPeriIdx(0),
        m_rotAng(0.0f),
        m_xShift(0.0f), m_yShift(0.0f) {
    }

    void init(int triIdx);

    bool fitsInBox();
    bool selfOverlaps() const;

    void nextGeneration(
        const SolutionKey& key, bool useTrim, SolutionMap& rNewSols) const;

    void print() const;

private:
    void addTriangle(
        const SolutionKey& key, int periIdx, int newTriIdx,
        SolutionMap& rNewSols) const;

    std::vector<int> m_periTris;
    std::vector<int> m_periLens;
    std::vector<bool> m_periFlipped;
    std::vector<Vec2> m_periPoints;

    int m_lastPeriIdx;

    std::vector<Vec2> m_triPoints;

    float m_rotAng;
    float m_xShift, m_yShift;
};

static int BoxW  = 0;
static int BoxH  = 0;
static int BoxD2 = 0;

static TriangleVec AllTriangles;
static SideMap AllSides;

Connector::Connector(
    int sideLen, int triIdx1, int triIdx2, bool flipped)
  : m_sideLen(sideLen),
    m_flipped(flipped) {
    if (triIdx1 < triIdx2) {
        m_triIdx1 = triIdx1;
        m_triIdx2 = triIdx2;
    } else {
        m_triIdx1 = triIdx2;
        m_triIdx2 = triIdx1;
    }
}

bool Connector::operator<(const Connector& rhs) const {
    if (m_sideLen < rhs.m_sideLen) {
        return true;
    } else if (m_sideLen > rhs.m_sideLen) {
        return false;
    }

    if (m_triIdx1 < rhs.m_triIdx1) {
        return true;
    } else if (m_triIdx1 > rhs.m_triIdx1) {
        return false;
    }

    if (m_triIdx2 < rhs.m_triIdx2) {
        return true;
    } else if (m_triIdx2 > rhs.m_triIdx2) {
        return false;
    }

    return m_flipped < rhs.m_flipped;
}

void SolutionKey::init(int triIdx) {
    m_tris.insert(triIdx);
}

void SolutionKey::add(int triIdx, const Connector& conn) {
    m_tris.insert(triIdx);
    m_conns.insert(conn);
}

bool SolutionKey::containsTriangle(int triIdx) const {
    return m_tris.count(triIdx);
}

int SolutionKey::minTriangle() const {
    return *m_tris.begin();
}

bool SolutionKey::operator<(const SolutionKey& rhs) const {
    if (m_tris.size() < rhs.m_tris.size()) {
        return true;
    } else if (m_tris.size() > rhs.m_tris.size()) {
        return false;
    }

    TriangleSet::const_iterator triIt1 = m_tris.begin();
    TriangleSet::const_iterator triIt2 = rhs.m_tris.begin();
    while (triIt1 != m_tris.end()) {
        if (*triIt1 < *triIt2) {
           return true;
        } else if (*triIt2 < *triIt1) {
           return false;
        }
        ++triIt1;
        ++triIt2;
    }

    if (m_conns.size() < rhs.m_conns.size()) {
        return true;
    } else if (m_conns.size() > rhs.m_conns.size()) {
        return false;
    }

    ConnectorSet::const_iterator connIt1 = m_conns.begin();
    ConnectorSet::const_iterator connIt2 = rhs.m_conns.begin();
    while (connIt1 != m_conns.end()) {
        if (*connIt1 < *connIt2) {
           return true;
        } else if (*connIt2 < *connIt1) {
           return false;
        }
        ++connIt1;
        ++connIt2;
    }

    return false;
}

void SolutionKey::print() const {
    TriangleSet::const_iterator triIt = m_tris.begin();
    while (triIt != m_tris.end()) {
        std::cout << *triIt << " ";
        ++triIt;
    }
    std::cout << "\n";

    ConnectorSet::const_iterator connIt = m_conns.begin();
    while (connIt != m_conns.end()) {
        connIt->print();
        ++connIt;
    }
    std::cout << "\n";
}

void SolutionData::init(int triIdx) {
    const Triangle& tri = AllTriangles[triIdx];

    m_periTris.push_back(triIdx);
    m_periTris.push_back(triIdx);
    m_periTris.push_back(triIdx);

    m_periLens.push_back(tri.getSideLenB());
    m_periLens.push_back(tri.getSideLenC());
    m_periLens.push_back(tri.getSideLenA());

    m_periFlipped.push_back(false);
    m_periFlipped.push_back(false);
    m_periFlipped.push_back(false);

    m_periPoints.push_back(Vec2(0.0f, 0.0f));
    m_periPoints.push_back(Vec2(tri.getSideLenB(), 0.0f));
    m_periPoints.push_back(Vec2(0.0f, tri.getSideLenA()));

    m_triPoints = m_periPoints;

    m_periPoints.push_back(Vec2(0.0f, 0.0f));
}

bool SolutionData::fitsInBox() {
    int nStep = 8;
    float angInc = 0.5f * M_PI / nStep;

    for (;;) {
        bool mayFit = false;
        float ang = 0.0f;

        for (int iStep = 0; iStep <= nStep; ++iStep) {
            float cosAng = cos(ang);
            float sinAng = sin(ang);

            float xMin = 0.0f;
            float xMax = 0.0f;
            float yMin = 0.0f;
            float yMax = 0.0f;
            bool isFirst = true;

            for (int iPeri = 0; iPeri < m_periLens.size(); ++iPeri) {
                const Vec2& pt = m_periPoints[iPeri];
                float x = cosAng * pt.x() - sinAng * pt.y();
                float y = sinAng * pt.x() + cosAng * pt.y();

                if (isFirst) {
                    xMin = x;
                    xMax = x;
                    yMin = y;
                    yMax = y;
                    isFirst = false;
                } else {
                    if (x < xMin) {
                        xMin = x;
                    } else if (x > xMax) {
                        xMax = x;
                    }
                    if (y < yMin) {
                        yMin = y;
                    } else if (y > yMax) {
                        yMax = y;
                    }
                }
            }

            float w = xMax - xMin;
            float h = yMax - yMin;

            bool fits = false;
            if ((BoxW >= BoxH) == (w >= h)) {
                if (w <= BoxW && h <= BoxH) {
                    m_rotAng = ang;
                    m_xShift = 0.5f * BoxW - 0.5f * (xMax + xMin);
                    m_yShift = 0.5f * BoxH - 0.5f * (yMax + yMin);
                    return true;
                }
            } else {
                if (h <= BoxW && w <= BoxH) {
                    m_rotAng = ang + 0.5f * M_PI;
                    m_xShift = 0.5f * BoxW + 0.5f * (yMax + yMin);
                    m_yShift = 0.5f * BoxH - 0.5f * (xMax + xMin);
                    return true;
                }
            }

            w -= 0.125f * w * angInc * angInc + 0.5f * h * angInc;
            h -= 0.125f * h * angInc * angInc + 0.5f * w * angInc;

            if ((BoxW < BoxH) == (w < h)) {
                if (w <= BoxW && h <= BoxH) {
                    mayFit = true;
                }
            } else {
                if (h <= BoxW && w <= BoxH) {
                    mayFit = true;
                }
            }

            ang += angInc;
        }

        if (!mayFit) {
            break;
        }

        nStep *= 4;
        angInc *= 0.25f;
    }

    return false;
}

static bool intersects(
    const Vec2& p1, const Vec2& p2,
    const Vec2& q1, const Vec2& q2) {

    if (cross(p2 - p1, q1 - p1) * cross(p2 - p1, q2 - p1) > 0.0f) {
        return false;
    }

    if (cross(q2 - q1, p1 - q1) * cross(q2 - q1, p2 - q1) > 0.0f) {
        return false;
    }

    return true;
}

bool SolutionData::selfOverlaps() const {
    int periSize = m_periPoints.size();

    int triIdx = m_periTris[m_lastPeriIdx];
    const Triangle& tri = AllTriangles[triIdx];
    float offsScale = 0.0001f / tri.getSideLenC();

    const Vec2& pt1 = m_periPoints[m_lastPeriIdx];
    const Vec2& pt3 = m_periPoints[m_lastPeriIdx + 1];
    const Vec2& pt2 = m_periPoints[m_lastPeriIdx + 2];

    Vec2 pt1o = pt1 + ((pt2 - pt1) + (pt3 - pt1)) * offsScale;
    Vec2 pt2o = pt2 + ((pt1 - pt2) + (pt3 - pt2)) * offsScale;
    Vec2 pt3o = pt3 + ((pt1 - pt3) + (pt2 - pt3)) * offsScale;

    float xMax = m_periPoints[0].x();
    float yMax = m_periPoints[0].y();
    for (int iPeri = 1; iPeri < m_periLens.size(); ++iPeri) {
        if (m_periPoints[iPeri].x() > xMax) {
            xMax = m_periPoints[iPeri].x();
        }
        if (m_periPoints[iPeri].y() > yMax) {
            yMax = m_periPoints[iPeri].y();
        }
    }

    Vec2 ptOut(xMax + 0.3f, yMax + 0.7f);
    int nOutInter = 0;

    for (int iPeri = 0; iPeri < m_periLens.size(); ++iPeri) {
        int iNextPeri = iPeri + 1;
        if (iPeri == m_lastPeriIdx) {
            ++iNextPeri;
        } else if (iPeri == m_lastPeriIdx + 1) {
            continue;
        }

        if (intersects(
            m_periPoints[iPeri], m_periPoints[iNextPeri], pt1o, pt3o)) {
            return true;
        }

        if (intersects(
            m_periPoints[iPeri], m_periPoints[iNextPeri], pt2o, pt3o)) {
            return true;
        }

        if (intersects(
            m_periPoints[iPeri], m_periPoints[iNextPeri], pt3o, ptOut)) {
            ++nOutInter;
        }
    }

    return nOutInter % 2;
}

void SolutionData::nextGeneration(
    const SolutionKey& key, bool useTrim, SolutionMap& rNewSols) const
{
    int nPeri = m_periLens.size();
    for (int iPeri = (useTrim ? 0 : m_lastPeriIdx); iPeri < nPeri; ++iPeri) {
        int len = m_periLens[iPeri];
        SideMap::const_iterator itCand = AllSides.lower_bound(len);
        SideMap::const_iterator itCandEnd = AllSides.upper_bound(len);
        while (itCand != itCandEnd) {
            int candTriIdx = itCand->second;
            if (!key.containsTriangle(candTriIdx) &&
                candTriIdx > key.minTriangle()) {
                addTriangle(key, iPeri, candTriIdx, rNewSols);
            }
            ++itCand;
        }
    }
}

void SolutionData::print() const {
    float cosAng = cos(m_rotAng);
    float sinAng = sin(m_rotAng);

    int nPoint = m_triPoints.size();

    for (int iPoint = 0; iPoint < nPoint; ++iPoint) {
        const Vec2& pt = m_triPoints[iPoint];
        float x = cosAng * pt.x() - sinAng * pt.y() + m_xShift;
        float y = sinAng * pt.x() + cosAng * pt.y() + m_yShift;
        std::cout << "(" << x << ", " << y << ")";

        if (iPoint % 3 == 2) {
            std::cout << std::endl;
        } else {
            std::cout << " ";
        }
    }
}

void SolutionData::addTriangle(
    const SolutionKey& key, int periIdx, int newTriIdx,
    SolutionMap& rNewSols) const {

    int triIdx = m_periTris[periIdx];
    bool flipped = m_periFlipped[periIdx];
    int len = m_periLens[periIdx];

    Connector conn1(len, triIdx, newTriIdx, flipped);
    SolutionKey newKey1(key);
    newKey1.add(newTriIdx, conn1);
    bool isNew1 = (rNewSols.find(newKey1) == rNewSols.end());

    Connector conn2(len, triIdx, newTriIdx, !flipped);
    SolutionKey newKey2(key);
    newKey2.add(newTriIdx, conn2);
    bool isNew2 = (rNewSols.find(newKey2) == rNewSols.end());

    if (!(isNew1 || isNew2)) {
        return;
    }

    SolutionData data;

    int periSize = m_periLens.size();
    data.m_periTris.resize(periSize + 1);
    data.m_periLens.resize(periSize + 1);
    data.m_periFlipped.resize(periSize + 1);
    data.m_periPoints.resize(periSize + 2);
    for (int k = 0; k <= periIdx; ++k) {
        data.m_periTris[k] = m_periTris[k];
        data.m_periLens[k] = m_periLens[k];
        data.m_periFlipped[k] = m_periFlipped[k];
        data.m_periPoints[k] = m_periPoints[k];
    }
    for (int k = periIdx + 1; k < periSize; ++k) {
        data.m_periTris[k + 1] = m_periTris[k];
        data.m_periLens[k + 1] = m_periLens[k];
        data.m_periFlipped[k + 1] = m_periFlipped[k];
        data.m_periPoints[k + 1] = m_periPoints[k];
    }
    data.m_periPoints[periSize + 1] = m_periPoints[periSize];

    data.m_lastPeriIdx = periIdx;

    data.m_periTris[periIdx] = newTriIdx;
    data.m_periTris[periIdx + 1] = newTriIdx;

    int triSize = m_triPoints.size();
    data.m_triPoints.resize(triSize + 3);
    for (int k = 0; k < triSize; ++k) {
        data.m_triPoints[k] = m_triPoints[k];
    }

    const Triangle& tri = AllTriangles[newTriIdx];
    int lenA = tri.getSideLenA();
    int lenB = tri.getSideLenB();
    int lenC = tri.getSideLenC();

    const Vec2& pt1 = m_periPoints[periIdx];
    const Vec2& pt2 = m_periPoints[periIdx + 1];

    Vec2 v = pt2 - pt1;
    v.normalize();
    Vec2 vn(v.y(), -v.x());

    float dA = lenA;
    float dB = lenB;
    float dC = lenC;

    int len1 = 0, len2 = 0;
    Vec2 pt31, pt32;

    if (len == lenA) {
        len1 = lenB;
        len2 = lenC;
        pt31 = pt1 + vn * dB;
        pt32 = pt2 + vn * dB;
    } else if (len == lenB) {
        len1 = lenC;
        len2 = lenA;
        pt31 = pt2 + vn * dA;
        pt32 = pt1 + vn * dA;
    } else {
        len1 = lenA;
        len2 = lenB;
        pt31 = pt1 + v * (dA * dA / dC) + vn * (dA * dB / dC);
        pt32 = pt1 + v * (dB * dB / dC) + vn * (dA * dB / dC);
    }

    if (isNew1) {
        data.m_periLens[periIdx] = len1;
        data.m_periLens[periIdx + 1] = len2;
        data.m_periFlipped[periIdx] = false;
        data.m_periFlipped[periIdx + 1] = false;
        data.m_periPoints[periIdx + 1] = pt31;

        data.m_triPoints[triSize] = pt1;
        data.m_triPoints[triSize + 1] = pt31;
        data.m_triPoints[triSize + 2] = pt2;

        rNewSols.insert(std::make_pair(newKey1, data));
    }

    if (isNew2) {
        data.m_periLens[periIdx] = len2;
        data.m_periLens[periIdx + 1] = len1;
        data.m_periFlipped[periIdx] = true;
        data.m_periFlipped[periIdx + 1] = true;
        data.m_periPoints[periIdx + 1] = pt32;

        data.m_triPoints[triSize] = pt1;
        data.m_triPoints[triSize + 1] = pt32;
        data.m_triPoints[triSize + 2] = pt2;

        rNewSols.insert(std::make_pair(newKey2, data));
    }
}

static void enumerateTriangles() {
    for (int c = 2; c * c <= BoxD2; ++c) {
        for (int a = 1; 2 * a * a < c * c; ++a) {
            int b = static_cast<int>(sqrt(c * c - a * a) + 0.5f);
            if (a * a + b * b == c * c) {
                Triangle tri(a, b, c);

                int triIdx = AllTriangles.size();
                AllTriangles.push_back(Triangle(a, b, c));

                AllSides.insert(std::make_pair(a, triIdx));
                AllSides.insert(std::make_pair(b, triIdx));
                AllSides.insert(std::make_pair(c, triIdx));
            }
        }
    }
}

static void eliminateInvalid(SolutionMap& rSols) {
    SolutionMap::iterator it = rSols.begin();
    while (it != rSols.end()) {
        SolutionMap::iterator itNext = it;
        ++itNext;

        SolutionData& rSolData = it->second;

        if (!rSolData.fitsInBox()) {
            rSols.erase(it);
        } else if (rSolData.selfOverlaps()) {
            rSols.erase(it);
        }

        it = itNext;
    }
}

static void trimSolutions(SolutionMap& rSols, int trimCount) {
    if (trimCount >= rSols.size()) {
        return;
    }

    SolutionMap::iterator it = rSols.begin();
    for (int iTrim = 0; iTrim < trimCount; ++iTrim) {
        ++it;
    }

    rSols.erase(it, rSols.end());
}

static void nextGeneration(
    const SolutionMap& srcSols, bool useTrim, SolutionMap& rNewSols) {
    SolutionMap::const_iterator it = srcSols.begin();
    while (it != srcSols.end()) {
        const SolutionKey& solKey = it->first;
        const SolutionData& solData = it->second;
        solData.nextGeneration(solKey, useTrim, rNewSols);
        ++it;
    }
}

static void printSolutions(const SolutionMap& sols) {
    std::cout << std::fixed;
    std::cout.precision(3);

    SolutionMap::const_iterator it = sols.begin();
    while (it != sols.end()) {
        const SolutionKey& solKey = it->first;
        solKey.print();
        const SolutionData& solData = it->second;
        solData.print();
        std::cout << std::endl;
        ++it;
    }
}

int main(int argc, char* argv[]) {
    if (argc < 3) {
        std::cerr << "usage: " << argv[0] << " width height [trimCount]"
                  << std::endl;
        return 1;
    }

    std::istringstream streamW(argv[1]);
    streamW >> BoxW;
    std::istringstream streamH(argv[2]);
    streamH >> BoxH;

    int trimCount = 0;
    if (argc > 3) {
        std::istringstream streamTrim(argv[3]);
        streamTrim >> trimCount;
    }

    BoxD2 = BoxW * BoxW + BoxH * BoxH;

    enumerateTriangles();
    int nTri = AllTriangles.size();

    SolutionMap solGen[2];
    int srcGen = 0;

    for (int iTri = 0; iTri < nTri; ++iTri) {
        const Triangle& tri = AllTriangles[iTri];

        SolutionKey solKey;
        solKey.init(iTri);

        SolutionData solData;
        solData.init(iTri);

        solGen[srcGen].insert(std::make_pair(solKey, solData));
    }

    int level = 1;

    for (;;) {
        eliminateInvalid(solGen[srcGen]);
        std::cout << "level: " << level
                  << " solutions: " << solGen[srcGen].size() << std::endl;
        if (solGen[srcGen].empty()) {
            break;
        }

        if (trimCount > 0) {
            trimSolutions(solGen[srcGen], trimCount);
        }

        solGen[1 - srcGen].clear();
        nextGeneration(solGen[srcGen], trimCount > 0, solGen[1 - srcGen]);

        srcGen = 1 - srcGen;
        ++level;
    }

    printSolutions(solGen[1 - srcGen]);

    return 0;
}