# 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:

1. The language your code is written in
2. A command to run your code

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)


1. Peter Taylor  - 5731679685.5 bricks per minute
2. Steve Verrill - 680262303 bricks per minute
3. Beta Decay    - 0.5 bricks per minute

• Why is there an [underhanded] category if they all get put on hold. "There are either too many possible answers..." which is kinda to key to obfuscation. – AndoDaan Aug 26 '14 at 17:32
• According to mathworld.wolfram.com/EulerBrick.html no "perfect" Euler brick with integral internal diagonal is known. If anyone comes up with a viable method, the whole world maths community will congratulate them. Voting to reopen because I like to dream, but I would prefer a straight Euler Brick code challenge: (bricks found in a minute) / (bytes). – Level River St Aug 26 '14 at 19:06
• @AndoDaan Good underhanded questions ask you to do one (specific) thing and pretend to do another. Questions like this here are just saying "pretend you're doing this thing, but in fact you can do whatever you want". That's what's too broad. In addition, as Steve and Eric point out, it is unlikely that it's even possible to write a program which does correctly compute the brick which means there are even less invalid programs. – Martin Ender Aug 26 '14 at 22:02
• Okay, much better. One more thing now: isn't your challenge trivially solved by hardcoding the smallest Euler brick and spitting out 1,000,000 multiples of it? – Martin Ender Aug 27 '14 at 16:15
• You could probably fix that by having a scoring like "Find the most Euler bricks with side-length less than L within m minutes." Where you pick L and m such that just producing multiples of a few hardcoded Bricks will hit L way before hitting m, but such that it's unlikely that anyone can find all bricks less than L in the given time (and if people do, the tie-breaker is the actual time they took). – Martin Ender Aug 27 '14 at 16:23

## 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 <stdlib.h>
#include <inttypes.h>
#include <stdio.h>
#include <signal.h>
#include <string.h>

// 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 < 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(&act, 0, sizeof(act));
act.sa_sigaction = &handler;
act.sa_flags = SA_SIGINFO;
if (sigaction(SIGINT, &act, NULL) < 0) return 1;

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

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

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

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

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

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

// Saunderson brick
scalar x = t.b * ( ((t.a << 1) >= t.c) ? ((t.a << 1) - t.c) : (t.c - (t.a << 1)) ) * ((t.a << 1) + t.c);
scalar y = t.a * ((t.b << 1) - t.c) * ((t.b << 1) + t.c);
scalar z = t.a * t.b * t.c << 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 >> 64), (uint64_t)x);
printf(", 0x%016"PRIx64"%016"PRIx64, (uint64_t)(y >> 64), (uint64_t)y);
printf(", 0x%016"PRIx64"%016"PRIx64, (uint64_t)(z >> 64), (uint64_t)z);
printf("; diags 0x%016"PRIx64"%016"PRIx64, (uint64_t)(u >> 64), (uint64_t)u);
printf(", 0x%016"PRIx64"%016"PRIx64, (uint64_t)(v >> 64), (uint64_t)v);
printf(", 0x%016"PRIx64"%016"PRIx64"\n", (uint64_t)(w >> 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.

• Man you did your research. I had a 'similar' approach, but it was about 8 orders of magnitude more naive. I was using a similar process to generate co-prime integers which can be used to generate pythagorean triples, and then comparing linear combinations of the triples to find cuboids. I guess I didn't hit google hard enough ;) – stokastic Aug 29 '14 at 0:16
• Which version of GCC are you using because I'm getting a load of errors when I try to compile your code? – Beta Decay Aug 29 '14 at 12:49
• Here is what I'm getting. – Beta Decay Aug 29 '14 at 12:54
• @BetaDecay, GCC 4.6, 64-bit Linux. There seem to be two important differences between my environment and yours: firstly, your target apparently doesn't support __int128 (are you using a 32-bit OS?); secondly, you don't apparently have a POSIX-compliant signal.h. – Peter Taylor Aug 29 '14 at 13:24
• This is amazing. A couple of comments: 1.Your non-Saunderson brick isn't primitive, it's divisible by ab. 2.With recursion you can avoid the memory issue (though the brick sizes would be all over the place. Considering both of these points would increase the number of (sufficiently small) bricks found. I'm still debating whether to do my own program. – Level River St Aug 30 '14 at 23:34

# C using Saunderson's parameterisation and a different tree generation method.

This is heavily influenced by Peter Taylor's answer, and by the following page: https://sites.google.com/site/tpiezas/0021. For consistency, my variable naming convention largely follows Peter's.

Instead of using a buffer, I use recursion to generate the three children of each node of the tree. The tree I generate is equivalent to Berggren's tree, but I use a two-parameter method (mainly for interest and variation.) The parameters p and q for all integers generate a pythagorean triple a=p^2-q^2, b=2pq, c=p^2+q^2. If we want to limit this to primitive triples only, we start with the vector (2,1) and multiply by the 2x2 matrices below. This is equivalent to Berggren's / Barning's method, in which we start with the vector (3,4,5) and multiply by the 3x3 matrices below. See http://en.wikipedia.org/wiki/Tree_of_primitive_Pythagorean_triples .

Barning's matrices                       2x2 matrices
1 -2  2      1  2  2     -1  2  2        2 -1    2  1    1  2
2 -1  2      2  1  2     -2  1  2        1  0    1  0    0  1
2 -2  3      2  2  3     -2  2  3


It can be seen from Barning's matrices that the hypotenuse c of the second child of each node will be between 3 and 7 times larger than that of its parent (more advanced treatment would narrow this range.) the growth rates of the other branches are rather slower.

In practice it was necessary to limit the depth of recursion in order to avoid stack overflow (segmentation fault.) However, it seems likely that only a few bricks were lost. Peter's program runs on my machine in about 6 minutes and generates 843035517 Euler bricks. My program takes just over twice as long and generates 1688135295 bricks with recursion depth limited to 5000, or 1671621046 bricks with depth limited to 1000. Therefore the rate of generation is about the same. With depth limited to 40, 1134616080 bricks are generated, of which 1131433861 are Saunderson and only 3182219 are conjugate. And the program runs about twice as fast with the code for generating conjugate bricks commented out!

Instead of generating the conjugate brick of (x,y,z) as (yz,xz,xy) I generate a conjugate brick ab times smaller. To do this intermediate values x',y',u',v' are stored during the calculation of x,y,z,u,v,w.

#include <inttypes.h>

int writeStdOut,count=0;

typedef __int128 scalar;

scalar L=(scalar)1000000*1000000*1000000*1000000*1000000,T=(scalar)1<<124;

f(scalar p, scalar q, int depth){
scalar p2,q2,
a,b,c,a2,b2,c2,
uprime,vprime,u,v,w,
xprime,yprime,x,y,z;

p2=p*p;q2=q*q;
a=p2-q2;b=p*q<<1;c=p2+q2;

z=a*b*c<<2;
if(z<L){
a2=a*a;b2=b*b;c2=c*c;

xprime=(b2<<2)-c2;yprime=(a2<<2)-c2;
if(xprime<0)xprime=-xprime;if(yprime<0)yprime=-yprime;
x=a*xprime;y=b*yprime;

uprime=b2*5+a2;vprime=a2*5+b2;
u=a*uprime;v=b*vprime; w=c2*c;

count++;
if (writeStdOut) {
printf("Sides: 0x%016"PRIx64"%016"PRIx64, (uint64_t)(x >> 64), (uint64_t)x);
printf(", 0x%016"PRIx64"%016"PRIx64, (uint64_t)(y >> 64), (uint64_t)y);
printf(", 0x%016"PRIx64"%016"PRIx64, (uint64_t)(z >> 64), (uint64_t)z);
printf("; diags 0x%016"PRIx64"%016"PRIx64, (uint64_t)(u >> 64), (uint64_t)u);
printf(", 0x%016"PRIx64"%016"PRIx64, (uint64_t)(v >> 64), (uint64_t)v);
printf(", 0x%016"PRIx64"%016"PRIx64"\n", (uint64_t)(w >> 64), (uint64_t)w);
}
if(T/a/c>xprime && T/b/c>yprime){
x=a*c*xprime<<2;y=b*c*yprime<<2; if(x<0 | y<0)puts("overflow");
if(x<L && y<L){

z=xprime*yprime;
u=uprime*yprime;v=vprime*xprime;w=c2*c2<<2;
count++;
if (writeStdOut) {
printf("SIDES: 0x%016"PRIx64"%016"PRIx64, (uint64_t)(x >> 64), (uint64_t)x);
printf(", 0x%016"PRIx64"%016"PRIx64, (uint64_t)(y >> 64), (uint64_t)y);
printf(", 0x%016"PRIx64"%016"PRIx64, (uint64_t)(z >> 64), (uint64_t)z);
printf("; diags 0x%016"PRIx64"%016"PRIx64, (uint64_t)(u >> 64), (uint64_t)u);
printf(", 0x%016"PRIx64"%016"PRIx64, (uint64_t)(v >> 64), (uint64_t)v);
printf(", 0x%016"PRIx64"%016"PRIx64"\n", (uint64_t)(w >> 64), (uint64_t)w);
}

}
}
if(depth<40){f((p<<1)-q,p,depth+1);f((p<<1)+q,p,depth+1);f(p+(q<<1),q,depth+1);}
}
}

main(int argc, char **argv){
writeStdOut=argc>1;
f(2,1,0);
printf(" %d bricks",count);
}

• Your answer deserves to get you at least as much rep as mine got me. Thanks for pointing out the bugs and for suggesting a massive improvement. – Peter Taylor Sep 7 '14 at 23:17
• I had to modify your code to include a signal interrupt, but after 2 minutes your code calculated 680262303 bricks. – Beta Decay Sep 9 '14 at 13:04
• Thanks for the bonus @PeterTaylor! I had started thinking about some rather more naive ideas when I saw your answer. All I did was copy your maths. I've finally got round to running your code today. I will comment soon. – Level River St Sep 14 '14 at 21:20

# Python 3 - Brute forcing it

## Edit: This produces a measly 0.5 bricks per minute...

To get the ball rolling again, I decided to write a simple entry. As always with brute force methods, it's extremely slow.

def mult(x,y):
c=0
for i in range(len(x)):
if y[i]%x[i]==0:c+=1
if c==3:return True
return False

pyth=lambda x,y:True if (x**2+y**2)**0.5==int((x**2+y**2)**0.5) else False
prev=[]
ymax=245
zmax=269
try:
for x in range(44,10**30):
for y in range(117,ymax):
ymax+=1
for z in range(240,zmax):
zmax+=1
curr=[x,y,z]
eq=0
ml=0
for i in range(3):
for j in range(3):
if i!=j:
if curr[i]==curr[j]:
eq+=1
for k in prev:
if mult(curr, k):
ml+=1
if eq==0 and ml==0:
f=pyth(curr,curr)
s=pyth(curr,curr)
b=pyth(curr,curr)
if f and s and b:
prev.append(curr)
except KeyboardInterrupt:
print('\n\n'+str(len(prev)))