# How to improve the performance of solutions to Project Euler (#39)?

This is the problem 39 of Project Euler, which I asked in the chat room two days ago. My original code runs as slowly as snails, and finally I got two answers from JM and Rojo. Unfortunately, both of them are faster than mine but still not enough.

Here is the problem:

If p is the perimeter of a right angle triangle with integral length sides, {a,b,c}, there are exactly three solutions for p = 120.

{20,48,52}, {24,45,51}, {30,40,50}

For which value of p <= 1000, is the number of solutions maximised?

My terrible answer, runs more than one minute:

Last@Sort@ Table[{Select[
IntegerPartitions[i, {3}],
#[[1]] < #[[2]] + #[[3]] && #[[1]]^2 == #[[2]]^2 + #[[3]]^2 &]
// Length, i}, {i, 1, 1000}]


JM's suggestion (not finished):

Count[IntegerPartitions[900, {3}], tri_ /; First[tri] < Total[Rest[tri]]
&& Norm[Rest[tri]]^2 == First[tri]^2]//Timing


{1.248, 4}

This is from Rojo (not finished):

Transpose[IntegerPartitions[900, {3}]]^2 /. {a_, b_, c_} :>
Length@a - Total@Unitize[a - b - c]//Timing


{0.124, 4}

I have noticed that many answers based on other languages (e.g., C, C++,...) in the forum of Project Euler can cost just several seconds, so how fast can MMa get? Since there are so many programming styles here, are there any rules for us to choose different programming styles for different problems? For this question, which style is the best?

-
+1 on the question because it generated very interesting answers :) –  Rojo Aug 29 '12 at 15:58
For $p \leq 10^6$, the answer is always of the form $2^a \times 3^b \times 5 \times 7 \times 11 \times 13 \times \cdots \times p_r$ for some $(a,b,r)$. If we could figure out the pattern in the exponents $(a,b)$, there might be a much faster method. –  David Speyer Apr 25 '13 at 12:45
Here are the list of all "record setters", perimeters which achieve higher values than any previous perimeter, for $n \leq 10^6$: {12, 60, 120, 240, 420, 720, 840, 1680, 2520, 4620, 5040, 9240, 18480, 27720, 55440, 110880, 120120, 166320, 180180, 240240, 360360, 720720} –  David Speyer Apr 25 '13 at 12:45

We are challenged to determine "how fast MMa can get" and, in so doing, to suggest rules "to choose different programming styles." The original solution takes 116 seconds (on my machine). At the time the question was posted, the solution time had been reduced by a factor of 1000 (10 doublings of speed) to 0.124 seconds by suggestions from users in chat.

This solution takes 1300 microseconds (0.0013 seconds) on the same machine, for a further 100-fold speedup (another 7 doublings):

euler39[p_] := Commonest @ Flatten[Table[ l,
{m, 2, Sqrt[p/2]},
{n, Select[Range[m - 1], GCD[m #, m + #] == 1 && OddQ[m - #] &]},
{l, 2 m (m + n), p, 2 m (m + n)}]];
Timing[Table[euler39[1000], {i, 1, 1000}];]


{1.311, Null}

It scales nicely (changing Table to ParallelTable to double the speed on larger problems):

AbsoluteTiming[euler39p[10^8]]


{120.8409117, {77597520}}

That is almost linear performance.

Note the simplicity of the basic operations: this program could be ported to any machine that can loop, add, multiply, and tally (the square root can be eliminated). I estimate that an efficient compiled implementation could perform the same calculation in just a few microseconds, using just 500 bytes of RAM, for up to another nine doublings in speed.

This solution was obtained through a process that, in my experience, generalizes to almost all forms of scientific computation:

1. The problem was analyzed theoretically to identify an efficient algorithm. Resulting speed: 0.062 seconds.

2. A timing analysis identified an MMa post-processing bottleneck. Some tweaking of this improved the timing. Speed: 0.0036 seconds (3600 microseconds).

3. In comments, J.M. and Simon Woods each suggested better MMa constructs, together reducing the execution time to 2400 microseconds.

4. The MMa bottleneck was removed altogether by a re-examination of the algorithm and the data structure, achieving a final reduction to 1300 microseconds (and considerably less RAM usage).

Ultimately a speedup factor of 90,000 was achieved, and this was done solely by means of algorithmic improvements: none of it can be attributed to programming style. Better MMa programmers than me will doubtlessly be able to squeeze most of the next nine speed doublings by compiling the code and making other optimizations, but--short of obtaining a direct $O(1)$ formula for the answer (which more or less would circumvent the whole point of the exercise, which is to use the computer for investigating a problem rather than for mere implementation of a theory-derived solution)--no more real speedup is possible. Note that compilation would also take us out of the MMa way of doing computation and bring us down to the procedural level of C and other compiled code.

The important lesson of this experience is that algorithm design is paramount. Don't worry about programming style or tweaking code: use your mathematical and computer science knowledge to find a better algorithm; implement a prototype; profile it; and--always focusing on the algorithm--see what can be done to eliminate bottlenecks. In my experience, one rarely has to go beyond this stage.

Detail of the story, as amended several times during development of this solution, follow.

This problem invites us to learn a tiny bit of elementary number theory, in the expectation it can result in a substantial change in the algorithm: that's how to really speed up a computation. With its help we learn that these Pythagorean triples can be parameterized by integers $\lambda \gt 0$ and $m \gt n \gt 0$ with $m$ relatively prime to $n$. We may take $x = \lambda(2 m n)$, $y = \lambda(m^2-n^2)$, and $z = \lambda(m^2+n^2)$, whence the perimeter is $p = 2 \lambda m (m+n)$.

The restrictions imposed by $p\le 1000$ and the obvious fact that $p$ is even give the limits for a triple loop over the parameters, implemented in a Table command below. The rest can be done without much thought--inelegantly and slowly--with brute force post-processing to avoid double counting $(x,y,z)$ and $(y,x,z)$ as solutions, to gather and count the solutions for each $p$, and select the commonest one.

(Although a triple loop sounds awful--one's instinctive reaction is to recoil at what looks like a $O(p^3)$ algorithm--notice that $m$ cannot exceed $\sqrt{p/2}$ and $n$ must be smaller yet, leaving few options for $\lambda$ in general. This gives us something like a $O(p f(p))$ algorithm with $f$ slowly growing, which scales very well. This limitation in the loop lengths is the key to the speed of this approach.)

euler39[p_] := Module[{candidates, scores, best},
candidates =
Flatten[Table[{Through[{Min, Max}[2 l m n, l (m^2 - n^2)]], 2 l m (m + n)},
{m, 3, Floor[Sqrt[p/2]]},
{n, 1, m - 1},
{l, 1, If[GCD[m, n] > 1, 0,  Floor[p / (2 m (m + n))]]}], 2];
scores = {Last[Last[#]], Length[#]} & /@
DeleteDuplicates /@
Gather[candidates[[Ordering[candidates[[;; , 2]]]]], Last[#1] == Last[#2] &];
best = Max[scores[[;; , 2]]];
Select[scores, Last[#] >= best &]
];


The amount of speedup is surprising. Accurate timing requires repetition because the calculation is so fast:

Timing[Table[euler39[1000], {i, 1, 1000}]]


{3.619, {{{840, 8}}, {{840, 8}}, ...

I.e., the time to solve the problem is $0.0036$ seconds or $1/17000$ minutes. This makes larger versions of the problem accessible (using ParallelTable instead of Table to exploit some extra cores in part of the algorithm):

euler39[5 10^6] // AbsoluteTiming


{55.1441541, {{4084080, 168}}}

Even accounting for the parallelization, the timing is scaling nicely: it appears to be acting like $O(p\log(p))$. The limiting factor in MMA is RAM: the program needed about 4 GB for this last calculation and attempted to claim almost 20 GB for euler39[10^7] (but failed due to lack of RAM on this machine). This, too, could be streamlined if necessary using a more compact data structure, and perhaps could allow arguments up to $10^8$ or so.

Perhaps a solution that is faster yet (for smaller values of $p$, anyway) can be devised by factoring $p$, looping over the factors $\lambda$, and memoizing values for smaller $p$. But, at $1/300$ of a second, we have already achieved a four order of magnitude speedup, so it doesn't seem worth the bother.

Remarkably, this is much faster than the built-in PowersRepresentations solution found by Simon Woods.

### Edit

At this point, J.M. and Simon Woods weighed in with better MMa code, together speeding up the solution by 50% (see the comments). In pondering their achievement, and wondering how much further one could go, it became apparent that the bottleneck lay in the post processing to remove duplicates. What if we could generate each solution exactly once? There would no longer be any need for a complicated data structure--we could just tally the number of times each perimeter was obtained--and no post-processing at all.

To assure no duplication, we need to check that when generating a triple $\{x,y,p\}$ with $x^2 + y^2 = z^2$ and $x+y+z=p$, we do not later also generate $\{y,x,p\}$: that's how the duplicates arise. The initial effort tracked possible duplicates by forcing $x\le y$. The improved idea is to look at parity.

The parameter $\lambda$ is intended to be the greatest common divisor of $\{x,y,p\}$. When it is, $x=2 m n$ and $y = m^2-n^2$ must be relatively prime. Because $x$ is obviously even, $y$ must be odd: that uniquely determines which of these two numbers is $x$ and which is $y$. Therefore, we do not need to check for duplicates if, in the looping, we guarantee that $2 m n$ and $m^2-n^2$ are relatively prime. A quick way to check is that (a) $m n$ and $m+n$ are relatively prime and (b) $m$ and $n$ have opposite parity. Making this check is essentially all the work performed by the algorithm: the rest is just looping and counting.

By eliminating the check for duplicates, the new solution doubled the speed once more, from 2400 microseconds to 1300 microseconds. Where does it spend its time? For an argument $p$ (such as $1000$),

• Approximately $p/2$ calculations of a GCD (for the second loop over n).

• A loop of length $p/2 / (m(m+n))$ for each combination of $(m,n)$.

An easy upper bound for the total number of iterations is $\frac{p}{8}\log{p}$, demonstrating the $O(p\log{p})$ scaling. If we assume the GCD calculations take an average of $\log{p}$ arithmetic operations each, the total number of operations is less than $p\log{p}$ plus comparable loop overhead together with incrementing an array of counts. The post-processing would merely scan that array for the location of its maximum. At $3 \times 10^9$ operations per second and $p=10^3$, the timing for good compiled code would be 0.3 microseconds. Problems up to $p \approx 10^{10}$ could be handled in reasonable time (under a minute) and without extraordinary amounts of RAM.

-
Outstanding speed up! It appears to be slightly faster to replace Min and Max with Sort, and run DeleteDuplicates over the whole list without gathering. With[{mmax=1000},cand=DeleteDuplicates@Flatten[Table[{Sort[{2 l m n,l (m^2-n^2)}],2 l m (m+n)},{m,3,Floor[Sqrt[mmax/2]]},{n,1,m-1},{l,1,If[GCD[m,n]>1,0,Quotient[mmax,2 m (m+n)]]}],2]; Commonest[cand[[All,2]]]] –  Simon Woods Aug 29 '12 at 19:17
@whuber, excellent edit and the new code is even faster! 5 orders of magnitude speed increase must be a site record mustn't it? A minor point - in your first paragraph, the 0.124 seconds for Rojo's code is for a single value of p, not the whole problem. –  Simon Woods Aug 30 '12 at 16:07
@Simon Oops, I overlooked that. But it's ok; there's no need to press that point. (My personal credo is that in general only speedups of $200$x or greater are worth working hard on, and it looks like that's been achieved by many of the replies here.) –  whuber Aug 30 '12 at 16:18
NB Compiling the code gets another factor of 3 for large problems but is not parallelizable (for some reason). (E.g., the $10^8$ problem takes 38 seconds; the $5\times 10^8$ problem is now feasible and takes 196 seconds and allocates 2 GB RAM. Its solution, by the way, is 465,585,120 = $(2^5)(3^2)(5)(7)(11)(13)(17)(19)$. The MMa kernel crashes on the $10^9$ problem :-(. ). –  whuber Aug 31 '12 at 2:16
@J.M. Doesn't GCD[m #,m + #]==1 imply GCD[m, #]==1 ? Then why is the former preferred ? –  The-Ever-Kid Feb 8 '13 at 6:16

All primitive Pythagorean triples ${a,b,c}$ are given by the parametric equations $a=m^2-n^2$, $b=2 m n$, $c=m^2+n^2$, where $m$ and $n$ are relatively prime, positive integers of opposite parity with $m>n$. The (even) perimeter of the corresponding primitive triangle is $2 m (m+n)$.

RelativePrimesOppositeParity[m] finds such $n$ corresponding to the input integer $m$. The implementation with Complement is much faster than with CoprimeQ or GCD.

RelativePrimesOppositeParity[m_Integer] :=
If[OddQ[m],
Complement[Range[2,m-1,2],Apply[Sequence,
Map[Range[2 #, m - 1, 2 #] &, FactorInteger[m][[All, 1]]]]],
Complement[Range[1,m-1,2], Apply[Sequence,
Map[Range[#, m - 1, 2 #] &, FactorInteger[m][[2 ;;, 1]]]]]
]


If the primitive Pythagorean triangle has perimeter $p$, then perimeters of any similar Pythagorean triangles are $p$ multiplied by factors $k$, $1\le k\le\left\lfloor p_{max}/p\right\rfloor$.

Make a table $t$ of allowed pairs $m$ and $n$, with maximum $m_{max}=\left\lfloor\sqrt{p_{max}/2}\right\rfloor$. Reset variable $t$ to the Tally of primitive perimeters multiplied by factors $k$. Find the maximum number of duplicate perimeters. Pick the corresponding {perimeter,max}.

ProjectEuler0039[pmax_] := Block[{max,t},
t = Table[2m(m+n),{m,Floor[Sqrt[pmax/2]]},{n,RelativePrimesOppositeParity[m]}];
t = Tally[Flatten[t/.{p_:>p*Range[Floor[pmax/p]]}]];
max = Max[t[[All,2]]];
Pick[t, t[[All, 2]], max]]


On my machine with $p_{max}=5\times 10^6$ or $p_{max}=10^8$, ProjectEuler0039 was twice as fast as euler39 in whuber's answer (the quality of which was unsurpassable!).

I am learning a great deal from the superb questions and answers on this site. My sincere thanks to all who participate!

-

Not fastest, but gives an idea of how Solve might be used with tolerable efficiency.

In[181]:= Timing[
soln = Solve[{x^2 + y^2 - z^2 == 0, z >= y - 1 >= x - 1 >= 0,
0 <= x + y + z <= 1000}, {x, y, z}, Integers];]

(* Out[181]= {6.36, Null} *)

In[182]:= mostCommonPerimeter =
Last[SortBy[Tally[Map[Total, {x, y, z} /. soln]], Last]]

(* Out[182]= {840, 8} *)

-

Here's an approach based on finding all right angled triangles with a hypotenuse <=500 and measuring the perimeters. The answer is the Commonest perimeter which is less than 1000. This runs in about 1 second.

rats[n_] := DeleteDuplicates[
Cases[Divisors[n^2, GaussianIntegers -> True],
z_Complex /; Abs[z] == n :> Sort[{Re@z, Im@z, n}]]];

Commonest[Select[Flatten[Table[Total[rats[h], {2}], {h, 5, 500}]], 0 < # <= 1000 &]]

(* 840 *)


Update

It is a bit faster to use PowersRepresentations to get the right angled triangles. Here is the solution as a one-liner:

Commonest[Select[
Flatten[# + Total[Rest@PowersRepresentations[#^2, 2, 2], {2}] & /@ Range[500]],
0 < # <= 1000 &]]

-

If you run the problem by Solve, it helps you a little along the way of looking differently on the question:

Solve[a + b + c == p && a^2 + b^2 == c^2 && a > 0 && b > 0 && c > 0 &&
a > b, {a, b, c}, Integers]


So what we find is that we know a range within which a must lie, and for each potential a we can calculate b and c, but really we are only interested in whether b and c are integers above zero, so we look at what the requirements for this are. For b we know it will be positive, so we just need the fraction to give an integer, and for c we know it will be an integer we just need it to be above zero, so we have:

nextInt[n_]:=Round[n+(1-Mod[n,1])]
previusInt[n_]:=Round[n-If[#==0,1,#]&@Mod[n,1]]

Table[
Count[
IntegerQ[p ((# - p/2)/(# - p))] &&
p > -p ((# - p/2)/(# - p)) + #
&/@ Range[nextInt[p - Sqrt[p^2]/Sqrt[2]], previusInt[p/2]]
,True]
,{p, 1, 1000}]
//Ordering[#, -1]& // Timing


{0.891, {840}}

I'm sure that this can be improved upon perhaps with a little more mathematical reasoning. For one thing, you could only iterate over even p, but when comparing solutions for speed that seems like cheating.

Slightly more reasoning

Actually, we don't need to ensure c is positive, the previousInt function correctly ensures that a<p/2 and since a>b we know already that c>0. So we just need to verify that p(a-p/2) is divisible by a-p, the check was an error on my part while using a different rounding function.

Table[Count[
Divisible[p  (# - p/2), # - p] & /@
Range[nextInt[p - Sqrt[p^2]/Sqrt[2]], previusInt[p/2]]
,True]
, {p, 1, 1000}] // Ordering[#, -1] & // Timing


{0.547, {840}}

-

The following is a very straightforward algorithm, implemented using inner loop vectorization technique, which is an adaptation of my reply to this blog post:

fn =
Compile[{{p, _Integer}},
Module[{al = Range[p], bl = 1, c = Range[p], d = Range[p],
zeros = 0*Range[p],result = 0*Range[p], rctr = 0},
For[bl = 1, bl <= p/2 + 1, bl++,
c = p - al - bl;
d = UnitStep[-Abs[c*c - bl*bl - al*al]];
If[d != zeros,
result[[++rctr]] = First@c[[Pick[Range[Length[d]], d, 1]]];
]
];
Union@Take[result, rctr]]
];


This function returns a list of hypotenuse lengths for a given perimeter, for example:

fn[120]

(* {50,51,52} *)


What is remarkable about this function is that, due to the vectorized inner loop, compilation to C does not bring any speed enhancements.

Using the above function on my 6 cores, I get:

Position[#,Max[#]]&@ParallelMap[Length[fn[#]]&,Range[1000]]//AbsoluteTiming

(*
{0.7812500,{{840}}}
*)


while on a single core I get it in about 3.2 seconds.

-
What's the purpose of setting c = Range[p] as it is reset later? –  rcollyer Aug 30 '12 at 2:20
@rcollyer Type information. A bit of waste, I agree - {0} would have been enough. –  Leonid Shifrin Aug 30 '12 at 9:35
I figured it was type information, but why is it needed? That is what has me confused. –  rcollyer Aug 30 '12 at 11:28
@rcollyer Because otherwise the compiler would assign some default type (integer or real number, forgot which), and then will report type conflict when first seeing it being a list rather than a number. –  Leonid Shifrin Aug 30 '12 at 19:05