Distribution of Primes
>// closedll("Distribution of Primes");
>// tccompile "Distribution of Primes";
>dll("Distribution of Primes","dllprimes",1);
>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 ]
>n=100000000;
>tic; p=dllprimes(n); size(p), toc;
[ 1 5761455 ]
Used 2.745 seconds

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




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

>sum(1/p)-log(log(n))
0.261501242993
>log(n)*prod(1-1/p), exp(-gamma$)
0.561457220953
0.561459483567
>tic; length(primes(1000000)), toc;
78498
Used 0.062 seconds
>tic; length(dllprimes(1000000)), toc;
78498
Used 0.016 seconds
>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
>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
}
Examples Homepage