Another Euler Brick in the Wall A Euler Brick is a cuboid where the length of all the edges are integers and all of the diagonals of the faces are integers as well. All sides must also be different.



Your program has to find as many different Euler Bricks where the sides are less than L within the time M (L and M will be defined later).

All of your calculated bricks have to be unique. This means that if your dimensions are the same, just in a different order, it must not be counted.

To calculate a brick means to find the lengths of the edges and the lengths of the diagonals of the faces of a valid Euler Brick.

Your bricks must also be indivisible by each other. In the same way that the Pythagorean triplet 6, 8, 10 is divisible by 3, 4, 5, your bricks must not be divisible by each other.

To test your program, I will need to the following:


The language your code is written in
A command to run your code
Your code


Your score will be the average number of bricks your program calculates (running the program five times). To do this, I will run your program and end it after M minutes have elapsed. Once you have submitted your answer, I will add a comment on your answer telling what your score is. If you make an edit to your code, please notify me directly so I can find your score again.

For that reason, I need your program to print the number of bricks either at throughout the running or at the end. This can be done like so (in Python it would help if people told me how to do this in other languages for this example):

try:
    # Calculation goes here

except KeyboardInterrupt: # I will use Ctrl-C
    print(numberofbricks)
    exit()


For this challenge L=10^30 and M=2 minutes. In the event of a tie, the actual time taken to calculate all of the bricks will be used. The lower time wins.

I will find your time using the script below. All programs will be tested on the same machine.

import time, subprocess
command = # How do you run your program?
s=time.time()
subprocess.Popen(command, shell=True)
e=time.time()
print(e-s)


Leaderboard

1. Peter Taylor  - 5731679685.5 bricks per minute
2. Steve Verrill - 680262303 bricks per minute
3. Beta Decay    - 0.5 bricks per minute
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.