# 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}],
#[] < #[] + #[] && #[]^2 == #[]^2 + #[]^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 E 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 E 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, {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, {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

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],
0 < # <= 1000 &]]


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

(* {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]//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

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], 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], previusInt[p/2]]
,True]
, {p, 1, 1000}] // Ordering[#, -1] & // Timing


{0.547, {840}}

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

In:= 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= {6.36, Null} *)

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

(* Out= {840, 8} *)


Here are the list of all "record setters", perimeters which achieve higher values than any previous perimeter, for n ≤ 10^32:

{12 60 120 240 420 720 840 1680 2520 4620 5040 9240 18480 27720 55440 110880 120120 166320 180180 240240 360360 720720 1081080 1441440 2042040 2162160 2882880 3603600 4084080 6126120 12252240 18378360 24504480 36756720 49008960 61261200 73513440 77597520 116396280 183783600 232792560 349188840 465585120 698377680 931170240 1163962800 1396755360 1784742960 2327925600 2677114440 3491888400 4655851200 5354228880 8031343320 10708457760 16062686640 21416915520 26771144400 32125373280 48188059920 53542288800 77636318760 80313433200 107084577600 155272637520 232908956280 310545275040 465817912560 621090550080 776363187600 931635825120 1397453737680 1552726375200 1863271650240 2329089562800 3105452750400 4658179125600 4813451763120 7220177644680 9316358251200 9626903526240 14440355289360 19253807052480 24067258815600 28880710578720 43321065868080 48134517631200 57761421157440 72201776446800 96269035262400 144403552893600 178097715235440 288807105787200 356195430470880 505412435127600 534293145706320 712390860941760 890488576177200 1068586291412640 1602879437118960 1780977152354400 2137172582825280 2671465728531600 3561954304708800 5342931457063200 7302006324653040 8014397185594800 10685862914126400 14604012649306080 16028794371189600 18700260099721200 21371725828252800 21906018973959120 29208025298612160 32057588742379200 36510031623265200 43812037947918240 65718056921877360 73020063246530400 87624075895836480 109530094869795600 146040126493060800 219060189739591200 313986271960080720 328590284609386800 438120379479182400 627972543920161440 657180569218773600 766710664088569200 876240758958364800 941958815880242160 1255945087840322880 1314361138437547200 1533421328177138400 1569931359800403600 1883917631760484320 2825876447640726480 3066842656354276800 3139862719600807200 3767835263520968640 4709794079401210800 6279725439201614400 9419588158802421600 18839176317604843200 28258764476407264800 37678352635209686400 43958078074411300800 56517528952814529600 65937117111616951200 88544128692742763040 131874234223233902400 147573547821237938400 177088257385485526080 263748468446467804800 295147095642475876800 442720643463713815200 885441286927427630400 1328161930391141445600 1770882573854855260800 2066029669497331137600 2656323860782282891200 3099044504245996706400 4692838820715366441120 5312647721564565782400 6198089008491993412800 7821398034525610735200 9385677641430732882240 12396178016983986825600 15642796069051221470400 23464194103576832205600 46928388207153664411200 70392582310730496616800 93856776414307328822400 109499572483358550292800 140785164621460993233600 164249358725037825439200 276877490422206620026080 281570329242921986467200 328498717450075650878400 461462484037011033376800 492748076175113476317600 553754980844413240052160 656997434900151301756800 922924968074022066753600 985496152350226952635200 1313994869800302603513600 1384387452111033100130400 2768774904222066200260800 4153162356333099300391200 5537549808444132400521600 6460474776518154467275200 8306324712666198600782400 9690712164777231700912800 16612649425332397201564800 19381424329554463401825600 28149211526257673035984800 29072136494331695102738400 33779053831509207643181760 38762848659108926803651200 56298423052515346071969600 58144272988663390205476800 77525697318217853607302400 84447634578773019107954400 168895269157546038215908800 253342903736319057323863200 337790538315092076431817600 394088961367607422503787200 506685807472638114647726400 591133442051411133755680800 1013371614945276229295452800 1182266884102822267511361600 1773400326154233401267042400 2263196606711116912093177920 2364533768205644535022723200 3546800652308466802534084800 4729067536411289070045446400 5657991516777792280232944800 10640401956925400407602254400 11315983033555584560465889600 16973974550333376840698834400 22631966067111169120931779200 26403960411629697307753742400 33947949100666753681397668800 39605940617444545961630613600 67895898201333507362795337600 79211881234889091923261227200 118817821852333637884891840800 158423762469778183846522454400 237635643704667275769783681600 316847524939556367693044908800 401717397691223251896539080800 475271287409334551539567363200 712906931114001827309351044800 792118812348890919232612272000 803434795382446503793078161600 1205152193073669755689617242400 1606869590764893007586156323200 1874681189225708508850515710400 2410304386147339511379234484800 2812021783838562763275773565600 4820608772294679022758468969600 5624043567677125526551547131200 8436065351515688289827320696800 11248087135354251053103094262400 16872130703031376579654641393600 22496174270708502106206188524800 29325370031459297388447352898400 33744261406062753159309282787200 50616392109094129738963924180800 56240435676771255265515471312000 58650740062918594776894705796800 87976110094377892165342058695200}


And here is the pseudo-code.

New Maximum (up to 6279725439201614400 = 77 sec in c#/64 bits)
This list should be put in 'The On-Line Encyclopedia of Integer Sequences (OEIS)'
The closest OEIS for this list is http://oeis.org/A099829

117,301,480,125,837,189,553,789,411,593,600 is high composite and contains :
25165824 divisors
2005111697 Pythagorean triplets

Algorithm: ∼ O(N^(1/3)) , a 2-fold time increase for 10-fold range

The perimeter P as a Pythagorean triplets formula:
P = 2M(M+N)
P/2 = M(M+N) , M and (M+N) are divisors of P/2

Then to find out the number of Pythagorean triplets of P
is to find the divisors of P/2 that follow these rules
M and N are co-primes
M > N , M > 1 , N > 0
M+N is odd , M+N > 1

Example: P = 420 has a total of 5 Pythagorean triplets
the divisors of P/2 = M(M+N) = 210 are:
1 2 3 5 6 7 10 14 15 21 30 35 42 70 105 210

Loop for each odd divisors (M+N) of 210 = 3 5 7 15
Lower bound: (M+N) > 1
Upper bound: M(M+M) > P/2 , 2*M*M > P/2 , 4*M*M > P , 2M > Sqrt(P) , (M+N) < Sqrt(P)
[maximum is 20 < Sqrt(420)]
Loop for each divisors M of 210
Lower bound: M > N , M+M > M+N , M > (M+N)/2
Upper bound: 0 < N , M < (M+N)
[(M+N) = 3 , M > 1 , M < 3 , M = 2]
[(M+N) = 5 , M > 2 , M < 5 , M = 3]
[(M+N) = 7 , M > 3 , M < 7 , M = 5]
[(M+N) = 7 , M > 3 , M < 7 , M = 6]
[(M+N) = 15 , M > 7 , M < 15 , M = 10]
[(M+N) = 15 , M > 7 , M < 15 , M = 14]
If [(M+N)-M] and M are co-primes Then 'triplet found'
[(M+N)-M] = N = 1 , M = 2 , co-primes (1)
[(M+N)-M] = N = 2 , M = 3 , co-primes (2)
[(M+N)-M] = N = 2 , M = 5 , co-primes (3)
[(M+N)-M] = N = 1 , M = 6 , co-primes (4)
[(M+N)-M] = N = 5 , M = 10 , not co-primes
[(M+N)-M] = N = 1 , M = 14 , co-primes (5)

As we found an efficient way to evaluate the number of Pythagorean triplets
of a specific perimeter P.
We also notice that P must be high composite to maximize
the number of divisors and the number of Pythagorean triplets of P

A pattern is found in the list of maximums numbers of Pythagorean triplets of P:
Step 1: For P >= 12 , for all new maximums , P are multiples of 12 (2*2*3)
Step 2: For P >= 60 , for all new maximums , P are multiples of 60 (2*2*3*5)
Step 3: For P >= 840 , for all new maximums , P are multiples of 420 (2*2*3*5*7)
Step 4: For P >= 9240 , for all new maximums , P are multiples of 4620 (2*2*3*5*7*11)
Step 5: For P >= 240240 , for all new maximums , P are multiples of 120120 (2*2*2*3*5*7*11*13)
Step x: For P >= ... and so on

To find new maximums (number of Pythagorean triplets of P):
At step 2 , evaluate new maximums , P in 60*[1..14] , without prime factors > 7
At step 3 , evaluate new maximums , P in 420*[2..22] , without prime factors > 11
At step 4 , evaluate new maximums , P in 4620*[2..52] , without prime factors > 13
At step 5 , evaluate new maximums , P in 120120*[2..34] , without prime factors > 17
...
Further optimizations has been implemented.

• Since your algorithm is not expressed in the Wolfram Language, it has no relevance to this site. – m_goldberg Jun 19 '15 at 4:50
• @m_goldberg I do not think I can agree with that. It certainly has some relevance, and unlike code in the wrong language it is pseudo-code as I requested here. While I would not encourage frequent posting of such answers in certain cases expertise from outside the Mathematica user community does benefit our users. If you feel strongly about this I encourage you to post (or point me to) a Meta question. – Mr.Wizard Jun 19 '15 at 6:09