Occasionally, I like to dabble around with mathematical code just to keep myself honest. So, I wrote a little function a couple of months back to calculate the next prime given a number n (iff n < 2³²). I showed this to the illustrious Andy Grimm, and of course he suggested that it may be possible to get the upper bound of the maximum possible factor without actually performing an expensive square root in the is_prime function. I recall from reading an embedded systems book a couple of years ago that something like this should indeed be possible, but a solution isn’t really obvious to me right now. At any rate, here’s my implementation of next_prime(). If you know of a faster way to compute is_prime(), let me know.

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

short isprime(unsigned long n) {
  unsigned long increment=4, divisor=5;
  unsigned long max_divisor = (unsigned long)sqrt(n);
  if ((n==1)||((n&2ul)==1)||((n%3ul)==0)||(n%5ul==0)) return 0;
  for(divisor=5; divisor<max_divisor; divisor+=increment) {
    if (n%divisor==0) return 0;
    increment = 6 - increment;
  }
  return 1;
}

unsigned long next_prime(unsigned long n) {
  unsigned long increment, divisor=5;
  if (n<7) {
    if (n==1) return 2;
    if (n==2) return 3;
    if (n<5) return 5;
  }
  while(n%6!=1 && n%6!=5) n+=1;
  increment = (n%6==1) ? 4 : 2;
  while(!isprime(n)) {
    n+=increment;
    increment = 6 - increment;
  }
  return n;
}

int main(int argc, char **argv) {
  printf("next prime is %d\n", next_prime(strtoul(argv[1], NULL, 10)));
}