# Finding prime numbers

• Difficulty: 015/100
• Author: jeremy theler
• Keywords: LOAD_ROUTINE, CALL, FUNCTION, ROUTINE, PRINT_FUNCTION, MIN, MAX, STEP, FORMAT,

Sometimes, the features provided by wasora are not enough to perform a certain computation. This example shows how wasora can be extended by using arbitrary snippets of code by dynamically loading a shared object file. The user-provided routines can be programmed in whatever language you may feel is suitable for the particular computation, as long as it is able to produce position-independent objects which could be linked into a dynamic library.

The examples below use C functions to test whether a certain number $$n$$ is prime or not. First, a file called prime.c is prepared with the following contents:

#include <stdio.h>
#include <math.h>

int isprimeint(int);
double isprime(const double *);
double nprimes(const double *);
double checkprime(const double *);

// check if an integer n is prime or not
int isprimeint(int n) {

int i;
int sqrtn;
int nmod6;

if (n <= 1) {
return 0;  // negative numbers, zero and one are not primes
} else if (n == 2 || n == 3) {
return 1;  // two and tree are primes
}

// any other prime is of the form 6k+1 or 6k-1 (prove it)
nmod6 = n % 6;
if (nmod6 == 1 || nmod6 == 5) {
// check for canditate divisors from 2 up to sqrt(n)+1
sqrtn = sqrt(n) + 1;
for (i = 2; i < sqrtn; i++) {
if (n % i == 0) {
return 0;
}
}
// none found, so n is prime
return 1;
} else {
return 0;
}

}

// wasora-like routine to be called without checking the return value
double checkprime(const double *x) {
int n = (int)(round(*x));
printf("%d is%s prime\n", n, isprimeint(n)?"":" not");
return 0;
}

// wasora-like function that returns true or false
double isprime(const double *x) {
return isprimeint((int)round(*x));
}

// wasora-like routine that return the number of primes
// smaller or equal than n
double nprimes(const double *x) {

int n = (int)(round(*x));
int np = 0;
int i;

for (i = 0; i < n; i++) {
np += isprimeint(i);
}

return (double)np;

}

A dynamic library called libprime.so can be generated with

$gcc -c prime.c -fPIC$ gcc -shared -o libprime.so prime.o

And then the individual routines can be either called or used as functions from within wasora.

## checkprimes.was

This example dynamically loads the the routine checkprime() from the shared object libprime.so and uses the keyword CALL to execute it, using a single argument. Remember that wasora works only with double-precision arithmetich, and so the loaded (and called) routines need as a prototype:

double routine(double *x);

If the routine does not take any argument, NULL will be passed. It is the userâ€™s responsibility to make sure the routine reads the very same number of arguments as passed from the wasora input. In this case, the function checkprime() takes one argument, which wasora evaluates and then passes it to the routine as a pointer to an array of doubles. When calling a routine, the return value is discarded, so the result of the call is directly printed to the standard output by the routine itself (see the source prime.c above).

LOAD_ROUTINE ./libprime.so checkprime

CALL checkprime 2
CALL checkprime 147
CALL checkprime 2371
CALL checkprime 2^16-1
CALL checkprime 2^17-1
$gcc -O2 -c prime.c -fPIC$ gcc -shared -o libprime.so prime.o
$wasora checkprimes.was 2 is prime 147 is not prime 2371 is prime 65535 is not prime 131071 is prime$ 

## isprime.was

We now load the routine isprime() that again takes a single argument but now it returns either a one or a zero. Therefore, we can define a wasora function isprime(n) of a single variable and use it as a regular wasora function. Note that the argument $$n$$ is rounded to the nearest integer before checking wether it is prime or not.

The example prints the wasora function isprime(n) for $$n=1,2,\dots,25$$ with PRINT_FUNCTION, and after that it evaluates the function at the very same points the previous example did with checkprime. Again, see prime.c for details of the implementation.

LOAD_ROUTINE ./libprime.so isprime
FUNCTION isprime(n) ROUTINE isprime

PRINT_FUNCTION isprime MIN 1 MAX 25 STEP 1 FORMAT %g

PRINT %g 147    isprime(147)
PRINT %g 2371   isprime(2371)
PRINT %g 2^16-1 isprime(2^16-1)
PRINT %g 2^17-1 isprime(2^17-1)
$wasora isprime.was 1 0 2 1 3 1 4 0 5 1 6 0 7 1 8 0 9 0 10 0 11 1 12 0 13 1 14 0 15 0 16 0 17 1 18 0 19 1 20 0 21 0 22 0 23 1 24 0 25 0 147 0 2371 1 65535 0 131071 1$ 

## nprimes.was

The last example of the section loads the routine nprimes() which counts how many prime numbers less or equal than $$n$$ there are. As with the previous example, a wasora routine is defined which can be again used in any context where a function could be used. In this case, the value of nprimes(n) is compared to two approximations of the number of prime numbers less or equal than $$n$$, namely

$\pi_1(n) \approx \frac{n}{\ln n}$

$\pi_2(n) \approx \int_2^n \frac{dt}{\ln t}$

LOAD_ROUTINE ./libprime.so nprimes
FUNCTION nprimes(n) ROUTINE nprimes

pi1(n) := n/log(n)
pi2(n) := integral(1/log(t), t, 2, n)

PRINT_FUNCTION nprimes pi1 pi2 MIN 1e5 MAX 2e6 STEP 1e5 FORMAT %g
$wasora nprimes.was | tee nprimes.dat 100000 9592 8685.89 9628.76 200000 17984 16385.3 18035 300000 25997 23787.7 26085.6 400000 33860 31009.6 33921.6 500000 41538 38102.9 41605.2 600000 49098 45096.9 49171.7 700000 56543 52010.4 56643.5 800000 63951 58856.6 64035.9 900000 71274 65644.8 71360.6 1e+06 78498 72382.4 78626.1 1.1e+06 85714 79075.1 85839 1.2e+06 92938 85727.6 93004.9 1.3e+06 100021 92343.5 100128 1.4e+06 107126 98926.1 107213 1.5e+06 114155 105478 114261 1.6e+06 121127 112002 121277 1.7e+06 128141 118499 128262 1.8e+06 135072 124971 135219 1.9e+06 142029 131421 142148 2e+06 148933 137849 149053$ qdp nprimes.dat --ti "real \$n/\\ln(n)\$ \$\\int_2^n\\frac{dt}{\\ln(t)}\$" --xlabel "\$n\$" --key "bottom"
\$