Another Euler Brick in the Wall
C using Saunderson's parameterisation and Berggren's tree

In 1740, Saunderson found the parameterisation that if (a, b, c) is a Pythagorean triple (i.e. a^2 + b^2 = c^2) then (a(2b-c)(2b+c), b(2a-c)(2a+c), 4abc) is a rational cuboid (this was before Euler studied them) with diagonals (c^3, a(5b^2+a^2), b(5a^2+b^2)). This doesn't generate all Euler bricks. Note that if the brick is primitive iff the Pythagorean triple is primitive.

If (x, y, z) is an Euler brick then so is (xy, xz, yz). Applying this once to a Saunderson brick creates a new brick which isn't Saunderson and which can be made primitive with division by ab. Applying the same transformation a second time gives the non-primitive xyz(x, y, z).

In 1934, Berggren discovered that the primitive Pythagorean triples form a ternary tree rooted at (3, 4, 5) which can be generated by a simple matrix multiplication and which contains each triple precisely once.

This implementation generates primitive Pythagorean triples and stores them in a stack, limiting the memory usage to 1.5GB. An earlier version used a circular buffer, but the stack gives much better results. This is probably a combination of fewer values being discarded and better cache locality.

Thanks to Steve Verrill for pointing out some bugs and for suggesting the use of a stack instead of a circular buffer. I've also followed him in removing the conjugate bricks for speed: previously the search space (minus the discarded branches) was exhausted within the two minutes, and there was value in spending time to get a few conjugates. But the new version can easily run for 5 minutes (as I discovered in the process of learning the value of the volatile keyword!), and most of the conjugates overflow, so on average the tests are wasted time.

#include &lt;stdlib.h&gt;
#include &lt;inttypes.h&gt;
#include &lt;stdio.h&gt;
#include &lt;signal.h&gt;
#include &lt;string.h&gt;

// Stack size. This is 1.5GB of memory. 
#define N 33554432

typedef unsigned __int128 scalar;

typedef struct {
    scalar a, b, c;
} triple;

// Stack operations
#define take(st) st[--stp];
#define offer(A,B,C) if (stp &lt; N) { st[stp].a=(A); st[stp].b=(B); st[stp].c=(C); stp++; }

volatile sig_atomic_t stop;

void handler(int sig, siginfo_t *siginfo, void *context) {
    stop = 1;
}

int main(int argc, char **argv) {
    // Register interrupt handler
    stop = 0;
    struct sigaction act;
    memset(&amp;act, 0, sizeof(act));
    act.sa_sigaction = &amp;handler;
    act.sa_flags = SA_SIGINFO;
    if (sigaction(SIGINT, &amp;act, NULL) &lt; 0) return 1;

    // L = 10^30
    const scalar L = ((scalar)0xc9f2c9cd0UL &lt;&lt; 64) + 0x4674edea40000000UL;
    const scalar T = (scalar)1 &lt;&lt; 64;   // used in an anti-overflow test

    // To avoid the optimiser removing the calculation of the bricks.
    int writeToStdOut = argc &gt; 1;

    triple *st = (triple *)malloc(N * sizeof(triple));
    if (!st) return 2;

    uint32_t stp = 1; // Index of first empty space in stack
    st[0].a = 3; st[0].b = 4; st[0].c = 5;
    uint64_t count = 0;

    while (stp &gt; 0) {
        if (stop) break;
        triple t = take(st);

        // Extend queue of Pythagorean triples
        offer(t.a + ((t.c - t.b) &lt;&lt; 1), ((t.a + t.c) &lt;&lt; 1) - t.b, ((t.a - t.b + t.c) &lt;&lt; 1) + t.c)
        offer(t.a + ((t.c + t.b) &lt;&lt; 1), ((t.a + t.c) &lt;&lt; 1) + t.b, ((t.a + t.b + t.c) &lt;&lt; 1) + t.c)
        offer(-t.a + ((t.c + t.b) &lt;&lt; 1), ((-t.a + t.c) &lt;&lt; 1) + t.b, ((-t.a + t.b + t.c) &lt;&lt; 1) + t.c)

        // Saunderson brick
        scalar x = t.b * ( ((t.a &lt;&lt; 1) &gt;= t.c) ? ((t.a &lt;&lt; 1) - t.c) : (t.c - (t.a &lt;&lt; 1)) ) * ((t.a &lt;&lt; 1) + t.c);
        scalar y = t.a * ((t.b &lt;&lt; 1) - t.c) * ((t.b &lt;&lt; 1) + t.c);
        scalar z = t.a * t.b * t.c &lt;&lt; 2;
        scalar u = t.c * t.c * t.c;
        scalar v = t.b * (t.b * t.b + 5 * t.a * t.a);
        scalar w = t.a * (t.a * t.a + 5 * t.b * t.b);
        count++;
        if (writeToStdOut) {
            printf("Sides: 0x%016"PRIx64"%016"PRIx64, (uint64_t)(x &gt;&gt; 64), (uint64_t)x);
            printf(", 0x%016"PRIx64"%016"PRIx64, (uint64_t)(y &gt;&gt; 64), (uint64_t)y);
            printf(", 0x%016"PRIx64"%016"PRIx64, (uint64_t)(z &gt;&gt; 64), (uint64_t)z);
            printf("; diags 0x%016"PRIx64"%016"PRIx64, (uint64_t)(u &gt;&gt; 64), (uint64_t)u);
            printf(", 0x%016"PRIx64"%016"PRIx64, (uint64_t)(v &gt;&gt; 64), (uint64_t)v);
            printf(", 0x%016"PRIx64"%016"PRIx64"\n", (uint64_t)(w &gt;&gt; 64), (uint64_t)w);
        }
    }

    printf("%ld\n", count);

    return 0;
}


Compiles with gcc -O3 euler_brick.c -o euler_brick and handles SIGINT (Ctrl-C).  On my computer the earlier circular buffer version ran out of usably small Pythagorean triples after 85 seconds having found 843035517 (0.843 thousand million) Euler bricks. The new version runs to the two minutes (tested with timeout -sINT 2m ./euler_brick) and finds about 9.1 thousand million bricks. Exact numbers vary from run to run; the highest I've observed is 9113459216.

Note that this is single-threaded code. The process is embarrassingly parallel - just start each thread going down a different branch of the tree - and so the limitation to scaling it to multicore would be the shared memory bus.