{var%20f='http://v.t.sina.com.cn/share/share.php?appkey=1515056452',u=z||d.location,p=['&url=',e(u),'&title=',e(t||d.title),'&source=',e(r),'&sourceUrl=',e(l),'&content=',c||'gb2312','&pic=',e(p||'')].join('');function%20a(){if(!window.open([f,p].join(''),'mb',['toolbar=0,status=0,resizable=1,width=440,height=430,left=',(s.width-440)/2,',top=',(s.height-430)/2].join('')))u.href=[f,p].join('');};if(/Firefox/.test(navigator.userAgent))setTimeout(a,0);else%20a();})(screen,document,encodeURIComponent,'','','https://www.manongdao.com/data/attach/logo/logo.png', '推荐 地球回转人心会变 的问题《prime number summing still slow after using sieve》','https://www.manongdao.com/q-1268003.html','页面编码gb2312|utf-8默认gb2312'));){kind=link}
I had a go at a project euler coding challenge below, the answer given by the code is correct but I do not understand why its taking nearly a minute to run. It was finishing with similar times prior to using a sieve. Others users are reporting times as low as milliseconds.
I assume I am making a basic error somewhere...
// The sum of the primes below 10 is 2 + 3 + 5 + 7 = 17.
// Find the sum of all the primes below two million.
public static long Ex010()
{
var sum = 0L;
var sieve = new bool[2000000];
var primes = new List<int>(10000);
for (int i = 2; i < sieve.Length; i++)
{
if (sieve[i-1])
continue;
var isPrime = true;
foreach (var prime in primes)
{
if (i % prime == 0) {
isPrime = false;
break;
}
}
if (isPrime) {
primes.Add(i);
sum += i;
for (var x = i * 2; x < sieve.Length; x += i) {
sieve[x-1] = true;
}
}
}
return sum;
}
EDIT:
The only thing that seemed missing was this optimization :
if (prime > Math.Sqrt(i))
break;
It brings the time down to 160 ms.
EDIT 2:
Finally clicked, took out the foreach as was suggested many times. Its now 12 ms. Final solution :
public static long Ex010()
{
var sum = 0L;
var sieve = new bool[2000000];
for (int i = 2; i < sieve.Length; i++)
{
if (sieve[i-1])
continue;
sum += i;
for (var x = i * 2; x < sieve.Length; x += i) {
sieve[x-1] = true;
}
}
return sum;
}
(tl;dr: 2 million in 0.8 ms, 2 billion in 1.25 s; segmented odds-only SoE, presieving, wheeled striding)
As always, the limit of Euler task #10 seems designed to pose a mild challenge on a ZX81, Apple ][ or C64 but on modern hardware you generally have to multiply the limits by 1000 to make things even remotely interesting. Or set a time limit like 5 seconds and try to see by how many orders of magnitude the Euler limit can be exceeded...
Dennis_E's solution is simple and efficient, but I'd recommend applying two small improvements that give a marked performance boost without any effort at all.
Represent only odd numbers in the sieve
All even numbers except for the number 2 are composite. If you pull the number 2 out of thin air when needed then you can drop all even numbers from the sieve. This halves workload and memory footprint, for a doubling of performance at the marginal cost of writing
<<or>>in a few places (to convert between the realm of numbers and the realm of bit indices). This is usually known as an 'odds-only sieve' or 'wheeled representation mod 2'; it has the added advantage that it largely removes the need for guarding against index overflow.Skip a few extra small primes during decoding
Skipping a few small primes ('applying a wheel') is much easier when going through a range of numbers incrementally, compared to hopping wildly about with different strides as during sieving. This skipping only involves applying a cyclic sequence of differences between consecutive numbers that are not multiples of the primes in question, like 4,2,4,2... for skipping multiples of 2 and 3 (the mod 6 wheel) or 6,4,2,4,2,4,6,2... for skipping multiples of 2, 3 and 5.
The mod 6 wheel sequence alternates between only two numbers, which can easily be achieved by XORing with a suitable value. On top of the odds-only sieve the distances are halved, so that the sequence becomes 2,1,2,1... This skipping reduces the work during decoding by 1/3rd (for stepping mod 3), and the skipped primes can also be ignored during the sieving. The latter can have a marked effect on sieve times, since the smallest primes do the greatest number of crossings-off during the sieving.
Here's a simple Sieve of Eratosthenes, with both suggestions applied. Note: here and in the following I generally go with the flow of C#/.Net and use signed integers where I would normally use unsigned integers in any sane language. That's because I don't have the time to vet the code for the performance implications (penalties) resulting from the use of unsigned types, like that the compiler suddenly forgets how to replace division by a constant with multiplication of the inverse and so on.
The first
ifstatement handles the small fry (n in 0..6) by returning a suitable element of a precomputed list of numbers, and it serves to get all the special cases out of the way in one fell swoop. All other occurrences of shift operators are for converting between the realm of numbers and the realm of indices into the odds-only sieve.This is pretty much the same code that I normally use for sieving small primes, up to 64K or so (the potential least factors for numbers up to 32 bits). It does Euler's piddling 2 million in 4.5 milliseconds but throwing bigger numbers at it shows its Achilles heel: it does a lot of striding over large distances, which interacts badly with modern memory subsystems where decent access speed can only be got from caches. The performance drops markedly when the capacity of the level 1 cache (typically 32 KiByte) is exceeded significantly, and it goes down even further when exceeding the L2 and L3 capacities (typically several megabytes). How sharp the drop is depends on the quality (price tag) of the computer, of course...
Here are some timings taken on my laptop:
In the middlish ranges there are time increases that go well beyond the expected factor of about 11, and then things stabilise again.
And here comes how to speed up the beast an order of magnitude...
Process the range in cache-sized segments
The remedy is easy enough: instead of striding each prime all the way from one end of the range to the other - and hence all across the memory space - we sieve the range in cache-sized strips, memorising the final positions for each prime so that the next round can continue right where the previous round left off. If we don't need a big bad sieve full of bits at the end then we can process a strip (extract its primes) after it has been sieved and then discard its data, reusing the sieve buffer for the next strip. Both are variations on the theme of segmented sieving but the offsets are treated differently during the processing; when the distinction matters then the first approach (big bad sieve for the whole range) is usually called a segmented sieve and the latter an iterated sieve. The terms 'moving' or 'sliding' sieve might fit the latter somewhat but should be avoided because they normally refer to a totally different class of sieves (also known as deque sieves) that are deceptively simple but whose performance is worse by at least an order of magnitude.
Here's an example of an iterated sieve, a slightly modified version of a function that I normally use for sieving primes in given ranges [m, n], just like in SPOJ's PRIMES1 and PRINT. Here the parameter m is implicitly 0, so it doesn't need to be passed.
Normally the function takes an interface that is responsible for processing the raw sieve (and any lose primes it may get passed), and which can get queried for the number of primes that the processor skips ('decoder order') so that the sieve can ignore those during the sieving. For this exposition I've replaced this with a delegate for simplicity.
The factor primes get sieved by a stock function that may look somewhat familiar, and I've changed the logic of the sieve from 'is_composite' to 'not_composite' (and to a base type that can participate in arithmetic) for reasons that will be explained later.
decoder_orderis the number of additional primes skipped by the decoder (which would be1for the function shown earlier, because it skips multiples of the prime 3 during the prime extraction/summing, over and above the wheel prime 2).Here's sieve processor that is equivalent to the logic used in the function shown at the beginning, and a dummy function that can be used to measure the sieve time sans any decoding, for comparison:
This function needs to perform a bit of modulo magic to synchronise itself with the mod 3 sequence over the mod 2 sieve; the earlier function did not need to do this because its starting point was fixed, not a parameter. Here are the timings:
Quite a difference, n'est-ce pas? But we're not done yet.
Replace conditional branching with arithmetic
Modern processors like things to be simple and predictable; if branches are not predicted correctly then the CPU levies a heavy fine in extra cycles for flushing and refilling the instruction pipeline. Unfortunately, the decoding loop isn't very predictable because primes are fairly dense in the low number ranges we're talking about here:
If the average number of non-primes between primes times the cost of an addition is less than the penalty for a mispredicted branch then the following statement should be faster:
This explains why I inverted the logic of the sieve compared to normal, because with 'is_composite' semantics I'd have to do
And the rule is to pull everything out of inner loops that can be pulled out, so that I simply applied
1 ^ xto the whole array before even starting.However, Euler wants us to sum the primes instead of counting them. This can be done in a similar fashion, by turning the value
1into a mask of all 1 bits (which preserves everything when ANDing) and 0 zeroises any value. This is similar to the CMOV instruction, except that it works even on the oldest of CPUs and does not require a reasonably decent compiler:Result:
Unrolling and strength reduction
Now, a bit of overkill: a decoder with a fully unrolled wheel mod 15 (the unrolling can unlock some reserves of instruction-level parallelism).
As you can see, I performed a bit of strength reduction in order to make things easier for the compiler. Instead of summing
window_base + (i << 1), I sumiand1separately and perform the rest of the calculation only once, at the end of the function.Timings:
The first set of timings is for
decoder_order == 1(i.e. not telling the sieve about the extra skipped prime), for direct comparison to the other decoder versions. The second set is fordecoder_order == 2, which means the sieve could skip the crossings-off for the prime 5 as well. Here are the null timings (essentially the sieve time without the decode time), to put things a bit into perspective:This shows that the work on the decoder has decreased decode time for 2 billion from 1.21 s to 0.35 s, which is nothing to sneeze at. Similar speedups can be realised for the sieving as well, but that is nowhere near as easy as it was for the decoding.
Low-hanging fruit: presieving
Lastly, a technique that can sometimes offer dramatic speedups (especially for packed bitmaps and/or higher-order wheels) is blasting a canned bit pattern over the sieve before comencing a round of sieving, such that the sieve looks as if it had already been sieved by a handful of small primes. This is usually known as presieving. In the current case the speedup is marginal (not even 20%) but I'm showing it because it is a useful technique to have in one's toolchest.
Note: I've ripped the presieving logic from another Euler project, so it doesn't fit organically into the code I wrote for this article. But it should demonstrate the technique well enough.
For a quick test you can hack the presieving into the sieve function as follows:
plus
plus
and
in the static constructor or Main().
This way you can use m_presieve_level for turning presieving on and off. The BlockCopy also works correctly after calling
set_presieve_level(0), though, because then the modulus is 1.m_wheel_ordershould reflect the actual wheel order (= 1) plus the decoder order; it's currently set to 3, so it'll work only with the v5 decoder at level 2.Timings:
You are doing trial division in addition to a sieve.
The boolean array will already tell you if a number is prime, so you don't need the List of primes at all.
You can also speed it up by only sieving up to the square root of the limit.
If you want to save some memory also, you can use a BitArray instead of a boolean array.