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"
$ 
nprimes-dat.png
nprimes-dat.png