Generating list of primes

I thought I'd share some basic ideas to make faster codes for some simple problems. One simple problem every post.

Generating list of all primes less than or equal to N

Simple enough problem, and I'm sure that Googling it will probably yield a host of solutions, but let's try coming up with a basic solution ourselves. So, to get started, a brute force solution would be to check divisibility of every number by numbers smaller than it.

This is Solution 0, which in C++ looks like this :-

vector <int> Primes;
for( int i = 2 ; i <= N ; i++) {
      bool isPrime = true;
      for( int j = 2 ; j < i ; j++)
            if( i % j == 0 ) {
                  isPrime = false;
                  break;
            }
      if( isPrime ) Primes.push_back( i );
}

On my poor decade-old computer (2.x GHz, 500MB RAM), when compiled with "gcc -O2", the above code takes 7s for N = 10^5. That is a crappy solution. We can improve it a lot.

To make it faster, we can ignore the even numbers and concern ourselves with just the odd numbers, but that only makes it slightly faster.
(Why only slightly faster? Aren't we removing half the numbers? Well yes, but it's not like an equal amount of time is spent on every number. Even numbers get discarded very quickly in the inner loop, the almost-primes take a lot more time to discard, and the primes take the most amount of time to confirm that there are no divisors.)
Also, since we are generating and saving a list of primes, we can use this list to speed things up - we don't need to divide by every single number but just smaller primes.

So we have Solution 1, which in C++ looks like this :-

vector <int> Primes;
Primes.push_back( 2 );
for( int i = 3 ; i <= N ; i += 2 ) {
      bool isPrime = true;
      for( vector <int> :: iterator j = Primes.begin() + 1 ; j != Primes.end() ; j++ )
            if( i % *j == 0 ) {
                  isPrime = false;
                  break;
            }
      if( isPrime ) Primes.push_back( i );
}

If we compare the run-times, we get:

	N = 10^5	N = 10^6
Solution 0	7s	~12 min
Solution 1	0.7s	~48s

This code is faster than the previous one, but this is still a terrible solution.

The above two solutions are very slow because their time complexity is essentially quadratic in N. (Actually, other than N^2, we also have some complicated terms of log(N) here, but the basic idea is that for every number, we are checking with all the primes smaller than it, which is (sort of) as bad as two loops of N).
This can be easily improved, since we only need to check with primes smaller than the square root of a number e.g. to check whether 101 is a prime, we only need to check with all primes up to 10 i.e. 2, 3, 5 and 7.

This gives us Solution 2, for which the C++ code is the same as above, except this line :

for( vector <int> :: iterator j = Primes.begin() + 1 ; j != Primes.end() && *j * *j <= i ; j++ )

Now we get:

	N = 10^5	N = 10^6	N = 10^7	N = 10^8
Solution 0	7s	~12 min
Solution 1	0.7s	~48s
Solution 2		0.2s	4.5s	~100s

The run-times show a sudden improvement as the time complexity comes down by a factor of N^1/2. (The asymptotic time complexity is an even more complicated expression now, and I'm not even gonna look into it now).

Now, how do we make this faster? Well, to reduce time from Solution 0 to Solution 2, instead of dividing by every smaller number, we were dividing only by primes smaller than the square root. But why go through this entire list every single time? I mean for every single number we first divide by 3, then 5, then 7, then 11, and so on until we find a divisor. For a number like 391 = 17*23, 5 iterations (3, 5, 7, 11, 13) are tried before the sixth iteration (17) marks this number as composite. Can't we save this effort? Isn't there any easier way to directly jump to its divisors 17 and 23?
Well, not really. It is a hard problem to find the factors of a given number. But we can approach this from another direction. It may be difficult to find divisors of a given number but it is easy to find multiples of a given divisor. For every prime number we find, we can find all its multiples in a single loop and discard these numbers as composite. Thus, now, to eliminate 391 = 17*23, we save the time of unnecessarily checking it with 3, 5, 7, 11 and 13 before it is eliminated by 17. As there are a lot of composite numbers to eliminate, this makes it considerably faster.
This solution is basically Sieve of Eratosthenes which was taught to us in middle school.

Our Solution 3 looks like this in C++ :-

vector <int> Primes;
Primes.push_back( 2 );
vector <bool> Sieve( N + 1 );
for( int i = 3 ; i <= N ; i += 2 ) Sieve[i] = 1;
for( int i = 3 ; i <= N ; i += 2 )
      if( Sieve[i] ) {
            Primes.push_back( i );
            for( int j = 3 * i ; j <= N ; j += 2 * i ) if( Sieve[j] ) Sieve[j] = 0;
      }

Note how we are still ignoring all the even numbers in our solution. Unlike our first solution, this actually makes a non-negligible difference right now. This is because unlike our first solution, the overhead involved with discarding each composite number is much lower, and getting rid of half the numbers makes it almost twice as fast.
Now we get:

	N = 10^5	N = 10^6	N = 10^7	N = 10^8	N = 10^9
Solution 0	7s	~12 min
Solution 1	0.7s	~48s
Solution 2		0.2s	4.5s	~100s
Solution 3			0.2s	4.2s	~60s

Now, how do we make this faster? Well, earlier in Solution 2, when we were checking divisibility by 2, 3, 5, 7 etc we noted that we only needed to check up to the square root of the number. Along similar lines of reasoning, we note that when we are discarding multiples of a prime p, we only need to consider the multiples greater than or equal to p^2, since the ones below have already been discarded by smaller primes.
This gives us our Solution 4, which in C++ looks the same as the previous one, the only change being inside this if block :

if( Sieve[i] ) {
      Primes.push_back( i );
      if( 1LL * i * i > N ) continue;
      for( int j = i * i ; j <= N ; j += 2 * i ) if( Sieve[j] ) Sieve[j] = 0;
}

Now we get:

	N = 10^5	N = 10^6	N = 10^7	N = 10^8	N = 10^9
Solution 0	7s	~12 min
Solution 1	0.7s	~48s
Solution 2		0.2s	4.5s	~100s
Solution 3			0.2s	4.2s	~60s
Solution 4			0.14s	2.4s	~28s

Now, the same question that we ask each time : how do we make this faster? Consider what happens when we discard multiples of 23. We start at 23x23, then 23x25, then 23x27, and so on. However, after 23x23, many of its multiples have already been discarded by smaller primes. 23x25 is discarded by 5, 23x27 by 3, 23x33 by 3, 23x35 by 5 and so on. So we could potentially save some time if we only focus on its multiples which have not already been discarded.These are numbers which are only divisible by 23 and larger primes.
This gives us our Solution 5, which in C++ looks like this :-

vector <int> Primes;
Primes.push_back( 2 );
vector <bool> Sieve( N + 1 );
for( int i = 3 ; i <= N ; i += 2 ) Sieve[i]=1;
for( int i = 3 ; i <= N ; i += 2 )
      if( Sieve[i] ) {
            Primes.push_back( i );
            if( 1LL * i * i > N ) continue;
            for( int j = N / i ; j >= i ; j-- ) if( Sieve[j] ) Sieve[i * j] = 0;
      }

Now we get:

	N = 10^5	N = 10^6	N = 10^7	N = 10^8	N = 10^9
Solution 0	7s	~12 min
Solution 1	0.7s	~48s
Solution 2		0.2s	4.5s	~100s
Solution 3			0.2s	4.2s	~60s
Solution 4			0.14s	2.4s	~28s
Solution 5				2s	~21s

Now, how do we make this faster? I'll continue this in Part 2.

Casablanca

Have you seen Casablanca ?

Classic Hollywood 1942 movie starring Humphrey Bogart and Ingrid Bergman. Definitely on any top 10 movie site. Simple plot, great actors, nice ending. Often been borrowed from, and its dialogues been referred to for decades after.

Bogart and Ingrid have a fabulous and rare chemistry, and Bogart's cold, cynical character has awesome witty dialogues - you can read nearly half the movie if you go to IMDB Memorable Quotes. Who hasn't heard the now-a-classic line "Play it Sam. Play As Time Goes By" ? And who can forget "Of all the gin joints ..." ?

Even the minor characters, played by actors like Peter Lorre and, especially, Claude Reins, are superb. ( "I think this is the beginning of a beautiful friendship ..." )

Well, all that is very well, but those are not the reasons I'm writing about Casablanca. My reason is Ingrid Bergman.

I swear, Ingrid in Casablanca is the most gorgeous woman I've ever seen on television.

Ingrid was a beautiful woman, no doubt, but in this movie she's just awesome. No other photograph, no other movie of her even comes close to her looks in this movie. And I don't mean that she wears a hot dress in this movie or anything like that; she just looks oh so beautiful. Every time I've seen the movie, every single time, just as soon as she entered the screen, I breathed "Wow!"

"Here's looking at you, kid !"


It's a good movie. If you haven't seen the movie, you might want to.
A Franc for your thoughts ?

p.s. It goes like
( Ilsa, to a pensive Rick: A Franc for your thoughts ?
Rick: In America they'd bring only a penny, and, huh, I guess that's about all they're worth.
Ilsa: Well, I'm willing to be overcharged. )

Arbit 2

In one of my internship interviews, I was asked a puzzle that I took over 10 minutes to answer:

"There are 100 prisoners about to be executed the next morning, in the following fashion - they are standed in a queue, and on every head a hat is placed, either white or black. (Number of white or black hats not known in advance). Now, starting from the last person, prisoner#100 to prisoner#1, each one says a colour of these two. (Every one can only see all the people in front but can hear what any person says.) Then only the prisoners who guessed their hat's colour correctly, survive. Rest die. They are told this the night before; come up with a strategy that saves as many of them as possible."

I started thinking with the easiest solution - even numbered person speaks the colour of the person in front, who repeats what they just heard. This way, 50 people definitely live, others have a 50-50 chance, so expected number of survivors is 75. But this was just the first thought, this can obviously be greatly improved.

( Solution follows this, so if you wanna think on your own first, do so before reading further. )

The solution I finally gave was that the prisoner#100 counts the white hats she can see, and says white if the number's even and black if it's odd. Prisoner#99 now counts the white hats she can see, and if her count is same (odd/even) as prisoner#100's count (odd/even), her hat must be black, otherwise white. So she says the colour of her hat. Prisoner#98 now counts the white hats (counts prisoner#99's hat too if it was white), compares this value with prisoner#100's count, and by same reasoning gets her colour. And so on to the first prisoner, who will only consider the colours of prisoners 99 to 2 when counting hats, and by same logic gets her colour.

So, the first 99 people definitely survive, and the last person has a 50-50 chance.

So far, so good. But later, after the interview, I started thinking, what if there were not 2 types of hats, but three - red, blue and green? Can the solution be modified to still save 99 people? What if there were n types of hats? Can 99 people be saved then? I found that a simple generalization serves the purpose.

( Again, if you want to think on your own, do so before reading further. )

Lets start with 3 colours of hats - Red, Blue and Green.
Prisoner#100 counts number of red hats (r) and the number of blue hats (b), and calculates m100 = (2*r + b) mod 3 . She then says R, B or G depending on whether m100 is 0, 1 or 2 respectively. Note that depending on this colour, everyone now knows m100. Prisoner#99 now counts hats and finds m99 = (2*r + b) mod 3. Now her colour is R, B or G if (m100-m99) mod 3 is 2, 1 or 0 respectively. Again, prisoner#98 counts hats, also considers prisoner#99's colour, then finds (m100-m98) mod 3, and by same reasoning finds her colour. And so on until the first person.

Just as before, first 99 prisoners live, and the last person has a 1 in 3 chance.

Similarly, for n colours - C1, C2... Cn. We just generalise this process.
Prisoner#100 counts numbers of hats of colour C1(r1), number of hats of colour C2(r2) and so on till Cn-1 (rn-1), and finds m100 = { (n-1)*r1+(n-2)*r2+... 1*rn-1 } mod n. She then says C1, C2... or Cn-1 if m100 is 0, 1, 2 .. or n-1 respectively. Prisoner#99 now counts hats and finds m99, and her colour is C1, C2... or Cn if (m100-m99) mod n is n-1, n-2 ... or 1 respectively. Again, prisoner#98 counts hats, also considers prisoner#99's colour, then finds (m100-m98) mod n, and by same reasoning finds her colour. And so on until the first person.

Once again, first 99 prisoners live, and the last person has a 1 in n chance.

The most interesting aspect of this generalization comes when the number of colours, n, tends to infinite. Just imagine the odds against the 100 men - they all have to guess the colour of their hats when the number of options to choose from is unimaginably huge. Despite the odds, their strategy lets the first 99 prisoners live. But the last person has a 1 in n chance of survival!

So, here's the funny part - infinite colours, yet 99 prisoners guess correctly thanks to the last person, who certainly dies. Funny, huh?

Arbit

For my first log, I'll tell ya about something arbitrary that I find cool.

As a child, you must've known the age-old children's puzzle - Draw the following without lifting pen from paper and without repeating any line.
The age old solution is - you fold the paper and use the backside as a bridge to get from one point to another. It's a very general trick, and can be used to draw anything without lifting pen from paper.

Let's look deeper. Using graph theory, it's easy to show that the above figure has no path that uses each edge exactly once. (Just read euler paths on wiki.) The basic idea is that if there's such a path, every time it visits a vertex, it must leave the vertex using a different edge, and so the degree of every vertex must be even, except for the starting and ending vertex (if they are not same), in which case, they both must have odd degrees. In the above figure, the four ends of the square all have odd degrees, which makes it impossible for such a path to exist.

Actually, the above condition on degrees can be shown to be both necessary and sufficient for an Euler path to exist. So, consider the following figure:
(Same as figure 1, except each line is drawn twice.) Clearly, every vertex has even degree, as I just doubled the degree of every vertex. So, this figure has an Euler cycle.

This gives us a new solution to the age-old problem of drawing figure 1. You draw figure 2 instead, and then use the axiom - "What I can do twice, I can do once" ! Neat huh? My new solution too is a general trick; it can be used to draw anything connected. Ta-da!