git.ozlabs.org Git - ccan/blob - ccan/eratosthenes/eratosthenes.c

   1 /* Licensed under LGPLv2+ - see LICENSE file for details */
   2 #include <ccan/eratosthenes/eratosthenes.h>
   3 #include <ccan/bitmap/bitmap.h>
   4
   5 #include <assert.h>
   6 #include <stdlib.h>
   7
   8 #define VAL_TO_BIT(v)           (((v) - 3) / 2)
   9 #define LIMIT_TO_NBITS(l)       ((l > 2) ? ((l) - 2) / 2 : 0)
  10
  11 #define BIT_TO_VAL(b)           (((b) * 2) + 3)
  12
  13 void eratosthenes_init(struct eratosthenes *s)
  14 {
  15         s->limit = 0;
  16         s->b = NULL;
  17 }
  18
  19 void eratosthenes_reset(struct eratosthenes *s)
  20 {
  21         if (s->b)
  22                 free(s->b);
  23         eratosthenes_init(s);
  24 }
  25
  26 static void eratosthenes_once(struct eratosthenes *s, unsigned long limit, unsigned long p)
  27 {
  28         unsigned long n = VAL_TO_BIT(3*p);
  29         unsigned long obits = LIMIT_TO_NBITS(s->limit);
  30
  31         if (obits > n) {
  32                 n = obits + p - 1 - ((obits - n - 1) % p);
  33         }
  34
  35         assert((BIT_TO_VAL(n) % p) == 0);
  36         assert((BIT_TO_VAL(n) / p) > 1);
  37
  38         while (n < LIMIT_TO_NBITS(limit)) {
  39                 bitmap_clear_bit(s->b, n);
  40                 n += p;
  41         }
  42 }
  43
  44 static void eratosthenes_sieve_(struct eratosthenes *s, unsigned long limit)
  45 {
  46         unsigned long p = 3;
  47
  48         while ((p * p) < limit) {
  49                 unsigned long n;
  50
  51                 eratosthenes_once(s, limit, p);
  52
  53                 n = bitmap_ffs(s->b, VAL_TO_BIT(p) + 1, LIMIT_TO_NBITS(limit));
  54
  55                 /* We should never run out of primes */
  56                 assert(n < LIMIT_TO_NBITS(limit));
  57
  58                 p = BIT_TO_VAL(n);
  59         }
  60 }
  61
  62 void eratosthenes_sieve(struct eratosthenes *s, unsigned long limit)
  63 {
  64         if ((limit < 3) || (limit <= s->limit))
  65                 /* Nothing to do */
  66                 return;
  67
  68         if (s->limit < 3)
  69                 s->b = bitmap_alloc1(LIMIT_TO_NBITS(limit));
  70         else
  71                 s->b = bitmap_realloc1(s->b, LIMIT_TO_NBITS(s->limit),
  72                                        LIMIT_TO_NBITS(limit));
  73
  74         if (!s->b)
  75                 abort();
  76
  77         eratosthenes_sieve_(s, limit);
  78
  79         s->limit = limit;
  80 }
  81
  82 bool eratosthenes_isprime(const struct eratosthenes *s, unsigned long n)
  83 {
  84         assert(n < s->limit);
  85
  86         if ((n % 2) == 0)
  87                 return (n == 2);
  88
  89         if (n < 3) {
  90                 assert(n == 1);
  91                 return false;
  92         }
  93
  94         return bitmap_test_bit(s->b, VAL_TO_BIT(n));
  95 }
  96
  97 unsigned long eratosthenes_nextprime(const struct eratosthenes *s, unsigned long n)
  98 {
  99         unsigned long i;
 100
 101         if ((n + 1) >= s->limit)
 102                 return 0;
 103
 104         if (n < 2)
 105                 return 2;
 106
 107         if (n == 2)
 108                 return 3;
 109
 110         i = bitmap_ffs(s->b, VAL_TO_BIT(n) + 1, LIMIT_TO_NBITS(s->limit));
 111         if (i == LIMIT_TO_NBITS(s->limit))
 112                 /* Reached the end of the sieve */
 113                 return 0;
 114
 115         return BIT_TO_VAL(i);
 116 }