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.
>// closedll("Distribution of Primes");
Press F10 in this line. If you have copied the C file to the current
directory (the one of the notebook), the predefined editor will open
the C file.
>// 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  5761455 ]

Used 2.745 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 states in terms of the prime counting function

Distribution of Primes

stating

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"); insimg;

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(1000000)), toc;
78498

Used 0.062 seconds
The C code is faster, of course. The gain is a factor of 3.
>tic; length(dllprimes(1000000)), toc;
78498

Used 0.016 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, using an array
for the odd numbers.
>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
The following listing is the C code. The main function is a bit
cluttered by error checks.
>open("Distribution of Primes.c"); ...
>repeat; getline(), until eof(); end; close();
#include "dlldef.h"
#include <string.h>

#define BYTE unsigned char

void makeprimes (BYTE *b, int k)
// make primes using the sieve of Eratosthenes
{
    int i,j,step;
    for (i=0; i<k; i++) b[i]=1;
    for (i=0; i<k; i++)
    {
        if (b[i]==1)
        {
            step=2*i+3; // b contains the odd numbers
            for (j=i+step; j<k; j+=step) b[j]=0;
        }
        if ((2*i+1)*(2*i+1)>2*k+1) break;
    }
}

EXPORT char * dllprimes (header * p[], int np,
    char *rstart, char *rend)
// Get and check the Euler parameters, run algorithm, and return result.
// See the HTML documentation for more information.
{
    start(rstart,rend); // necessary start routine
    
    double x;
    x=getreal(p[0]); // get the real number and check param type
    IFERROR("Need a real for dllprimes.");
    
    int n=(int)x; // convert to integer
    if (n!=x || n<=2) ERROR("Need positive integer larger than 2");
    
    // get byte array from stack
    int k=(n-3)/2+1;    
    BYTE *b=(BYTE *)getram(k*sizeof(BYTE));
    CHECK(b,"Out of heap space.");
    
    makeprimes(b,k); // run the sieve
    
    // count prime numbers
    int i,count=0; 
    for (i=0; i<k; i++) if (b[i]) count++;
    
    // make output vector on stack
    header *out;
    out=new_matrix(1,count+1);
    IFERROR("Out of stack space in dllprimes.");
    double *m=matrixof(out);
    m[0]=2;
    
    // copy primes into it
    count=1;
    for (i=0; i<k; i++) if (b[i]) { m[count]=(2*i+3); count++; }
    
    moveresults(out); // necessary moving of result
    
    return newram; // always return newram, ram or 0
}

The sieve is in the makeprimes function. 

The dllprimes function scans the argument, takes place for the sieve
from the stack, runs the sieve, and collects the result.

Examples Homepage