Faster Primes

Do you remember my post about prime numbers? If you don't, you can catch up clicking on the link:
Prime Numbers: A Computational Approach

I've improved the function that checks if a number is a prime number. Here's the new code:

/*
 * Check if "number" is a prime number. If it is, the function will return 1.
 * If it isn't, it will return 0.
 *
 */
int isPrime(long long unsigned number)
{
    if(number == 2)
    {
        return 1;
    }
    else if(number % 2 == 0)
    {
        return 0;
    }
    long long unsigned aux;
    for(aux = 3; aux*aux <= number; aux+=2)
    {
        if(number % aux == 0)
        {
            return 0;
        }
    }
    return 1;
}

The key is the double increment in the for-loop (aux+=2). Apart from 2, all prime numbers are odd numbers. The first thing I do, is check if a prime number is 2. Secondly, I discard the possibility that the number is divisible by 2. Finally, I proceed with the algorithm the same way I did in the first program, but this time I start checking if the number can be divided by 3. And then, I can safely increment the number two units on each iteration (if 2 cannot divide a number, a number that can be divided by 2 won't be able to divide that number either).

This screencaptures (it's about the Unix time command again) show this improvement. I calculate primes up to 10000000 this time.

Old version

New version

You can download the source code with the improved function here.

The maths of the EuroMillions

My uncle is very engaged in this game, and I've found a little inspiration to write about it. Lottery... I hate it. However, lets see what science has to say about this.

The player must choose 5 numbers between 1 and 50 and 2 stars (stars are numbers as well) between 1 and 11. Numbers cannot be repeated.

Let's see how many combinations of 5 numbers and 2 stars can we make:

We can choose 5 numbers in 50*49*48*47*46 = 254251200 different ways, but the order doesn't matter. There are 5! = 120 ways of ordering 5 numbers, so there are 254251200/120 = 2118760 ways of selecting the numbers without taking into account the order of selection.

We can choose 2 stars in 11*10 = 110 different ways, but (again) the order doesn't matter. There are 2! = 2 says of ordering 2 stars (which are also numbers), so there are 110/2 = 55 ways of selecting the numbers without taking into account the order of selection.

Combining those numbers and stars, there are a total of 116531800 ways of betting in EuroMillions.

Taking this as a starting point, I've made a simple program (a script written in Python) that outputs the following:

If I've made it correctly, as the output says you need to play 80773689 times in a row to make your chances higher than 50%. In Spain, each ticket costs 2€, which means that if you want to avoid losses, you should only play when the jackpot is 161547378€ or higher (162 million in easy terms).

Notes:
- Actually, chances of winning are higher, since there are other prices, perhaps I should do the maths in another post, but here you can get a general idea.
- When I say that you should play for those jackpots, chance is still very little and you should buy quite a lot of tickets. You may need several million years to make profit playing lottery.
- Every lottery company expects to make money, so it's not a surprise if you end up with losses.
- Some people play because when they buy a ticket, they're buying a dream.
- Maths know nothing about dreams, just numbers.

Warehouse simulation

Today I'll share with you a warehouse simulation I had to do a month ago.

Let's consider a warehouse of a specific product whose sale price is 2€ per unit. Customer visits are distributed like a Poisson process (0.5 customers per hour) and the amount of products that each customer wants follows this distribution:

Request	1 unit	2 units	3 units	4 units
Probability	0.3	0.4	0.2	0.1

In order to satisfy the demand, the owner of the warehouse keeps a stock of products. When the inventory is low, he asks the distributor for more units. The orders policy is periodic. When the inventory is fewer than 30 and there aren't any pending orders, the owner asks for more units in a way that the inventory would end up having 100 units.

In each order there's a base price of 10€ (it doesn't matter the amount of units the owner asked). Additionally the cost per unit included in the order will depend on the amount of units asked. If the amount of units requested is less than 50, the price will be 1€/unit. However if the owner asks 50 or more units, the price would be 75 cents per unit.

The delivery time follows a normal distribution (with a mean of 48 hours and a standard deviation of 0.8 hours). The delivery will be paid as soon as it arrives. An agreement with the provider entitles us with a discount of 0.01% for every three hours of delay in the delivery. On the other hand, it will be 0.01% more expensive for every three hours it arrives sooner.

The owner of the warehouse spends 0.1€ per product each hour (due to physical storage, refrigeration, …).

If the customer asks for a greater amount than what is available, the warehouse sells all it has.

a) Simulate the behaviour of the warehouse during 5 months. Estimate the expected profits, the proportion of completely satisfied customers and the percentage of time the inventory remains zero. Let's assume the initial inventory is 70 units of product.

b) Represent graphically the evolution of the inventory level during the mentioned 5 months.

c) Try to identify what changes you could implement to improve the behaviour of the inventory model (get more profits and more satisfied customers).

For the first question, I've programmed a Python script. In order to run it, you need to have installed the NumPy library. I'll show you the output of a few simulations:

 

This is the graph that represents the evolution of the inventory. For this question, I commented the final part of the code (when the results of the simulation are printed) and put at the end of each iteration a line of code that printed in one line the stock available. Then, I used LibreOffice to plot the graphic. Maybe I should use better tools and techniques the next time.

For the last question. There's an obvious way to reduce the losses: sell all the stock and forget about new orders. That way, we don't have order expenses and we keep maintenance to a minimum. Of course, we'll have a lot of potential losses (a lot of unsatisfied clients), but it's the best way to solve this problem that I've found so far (maintenance is too expensive!). Nevertheless, I challenge the reader to get better results. Mine are these:

Don't change the product sale price, that's cheating.

Prime Numbers: A Computational Approach

I have two obsessions: prime numbers and brute force algorithms.

The first ones are almost useless in my daily life, but they're somehow... beautiful.

Brute force is a very useful way to solve simple problems, because you don't have to think about how to get a result when the computer can do it for you. Well, as I said, the problem has to be simple. Today's post isn't the case, but it will be fun.

There's that urban legend about the FBI buying prime numbers over 100 digits, which I don't know if it's true. However, it certainly meets my two obsessions. So, I made a program that prints all the prime numbers that exist.

Of course, this is just for fun. I don't have the methods neither the necessary machinery to actually get over 100 digits in a reasonable period of time. But I want to share one optimization to the method I discovered in a book while I was teaching C programming to a student of Civil Engineering (the former "Ingeniería de Caminos Canales y Puertos").

So, the first function I ever programmed that checked if a given number is a prime number was the following:

int isPrime(long long unsigned number)
{
    long long unsigned aux;
    for(aux = 2; aux < number; aux++)
    {
        if(number % aux == 0)
        {
            return 0;
        }
    }
    return 1;
}

As you can see, the function checks if any number divides the given number one by one until it gets to the number given in the first place. Let's say: pure brute force.

Now let's see the first optimization (founded in the previously mentioned book):

int isPrimeFancy(long long unsigned number)
{
    long long unsigned aux;
    for(aux = 2; aux*aux < number; aux++)
    {
        if(number % aux == 0)
        {
            return 0;
        }
    }
    return 1;
}

What has changed is the condition of iteration in the for loop. Why can the loop stop there? Let's give an intuitive approach.

Imagine that we want to check if 19 is a prime number. The first method (which we're sure that works because is pure brute force) will iterate exactly:

(19 - 2) times = 17 times

The second algorithm will do something like this:

does 2 divide 19? No
does 3 divide 19? No
does 4 divide 19? No

Wait. The next number is 5. If 5 divided 19, another number should have appeared first, but it hasn't. Why? Because 5 times 5 is the minimum value that would match 19 here, and that's not the case, because 5 times 5 is 25, which it's over 19. In next iterations the minimum possible match would be 6 times 6. We're getting even farther from 19. This will happen with 7, 8 and the following numbers. So, those iterations will be certainly unnecessary.

Let's compare the performance of these two algorithms in this specific case. The first one will iterate 17 times, while the second one will iterate only 3 times. So, using the second algorithm we'll save 14 iterations, more than the 400% of the execution time! Imagine how much time could be saved on bigger numbers.

So, in the general case, we have the first algorithm with a theoretical complexity of O(n) and the second one with a complexity of O(sqrt(n)). Now it's time to try this for real.

Let me introduce you to my laptop:
HP Pavilion dv6700
Memory: 3GiB
Processor: Intel® Core™2 Duo CPU T5750 @ 2.00GHz × 2
Running: Ubuntu 12.04 Precise Pangolin

In order to calculate the real performance of these two algorithms in my laptop, I've used the Linux command time in it's most simple way:

time ./prime

These are the results for the first algorithm checking numbers up to 999999.



And here they are for the second.



I can assure you the main program is the same. It checks number by number if it's a prime, and if it is, the program prints it on the screen. You can download the complete source code of this program by clicking here. By the way, the main program in that link will iterate forever, so you will need to add a little piece of code to make it stop when you want.

That's a difference, is it not?

That's all for today. In later posts I might do some research about how to get this to the next level...

As I said in my first post, any suggestions, mistakes, comments... congratulations... are welcome.