iconEuler Examples

Distribution of Primes

This is a demonstration of the DLL feature in Euler. With DLLs, you can use C and even assembler code in Euler. Euler comes with the C compiler tinycc. You can compile files from the command line.

The following two lines are commented, since Euler cannot write to the installation directory. If you want to try the compilation, uncomment the lines and save the notebook to a directory with write access. For the source of the file, see

Source File

We copy the file to eulerhome(), which is usually "\Users\YourName\Euler". Compiling in the installation directory would cause an error.

>fn="Distribution of Primes.c"; filecopy(home()+fn,eulerhome()+fn); ...
 cd(eulerhome()); tccompile "Distribution of Primes";

The dll command loads the function dllprimes (which needs one parameter) from the DLL.

>dll("Distribution of Primes","dllprimes",1);

The function uses the sieve of Eratosthenes to compute all prime numbers up to a certain limit.

>dllprimes(100)
[2,  3,  5,  7,  11,  13,  17,  19,  23,  29,  31,  37,  41,  43,  47,
53,  59,  61,  67,  71,  73,  79,  83,  89,  97]

Let us try to find all prime numbers to 100 million.

>n=100000000;

tic and toc take the time. On my computer, it takes less than 3 seconds.

>tic; p=dllprimes(n); size(p), toc;
[1,  5.76146e+006]
Used 1.584 seconds

It is interesting to plot the distribution of the primes. It is known that

Distribution of Primes

where p[n] is the n-th prime number.

>k=1000:1000:cols(p); plot2d(k*log(p[k])/p[k],xl="*1000",>insimg);

Distribution of Primes

Usually this is stated in terms of the prime counting function

Distribution of Primes

Then

Distribution of Primes

But the above statement is equivalent, since n is the number of primes less than p[n]. For the function p[n] the following development is known

Distribution of Primes

We check the O-constant.

>plot2d((p[k]/(k*log(k))-1-(log(log(k))-1)/log(k))*log(k)^2,xl="*1000"):

Distribution of Primes

We can compute an approximation of the Meissel-Mertens constant.

>sum(1/p)-log(log(n))
0.261501242993

Or check Mertens second theorem. gamma$ is the Euler-Mascheroni constant.

>log(n)*prod(1-1/p), exp(-gamma$)
0.561457220953
0.561459483567

Comparison of Running Times

Euler has the sieve function. However, with the default stack size, only primes up to a million can be computed.

>tic; length(primes(100 000 000)), toc;
5761455
Used 1.966 seconds

The C code is faster, of course. The gain is a factor of 3.

>tic; length(dllprimes(100 000 000)), toc;
5761455
Used 1.53 seconds

Clearly, for other problems involving complicated data structures, Euler might be a completely wrong choice. Then the C interface is very nice. It can even keep data between subsequent calls.

The sieve function in Euler is straightforward coding in matrix languages, using an array for the odd numbers, and index arrays to avoid loops.

>type primes
function primes (n: integer)
    if n>=3
        len = floor((n-1)/2);
        sieve = ones ([1,len]);
        for i=1 to (sqrt(n)-1)/2;
            if (sieve[i]) then
                sieve[3*i+1:2*i+1:len] = 0; ..
            endif
        end
        return [2, 1+2*nonzeros(sieve)];
    elseif n>=2
        return 2;
    else
        return [];
    endif
endfunction

Examples