linux/lib/math/prime_numbers.c
<<
>>
Prefs
   1// SPDX-License-Identifier: GPL-2.0-only
   2#define pr_fmt(fmt) "prime numbers: " fmt
   3
   4#include <linux/module.h>
   5#include <linux/mutex.h>
   6#include <linux/prime_numbers.h>
   7#include <linux/slab.h>
   8
   9#define bitmap_size(nbits) (BITS_TO_LONGS(nbits) * sizeof(unsigned long))
  10
  11struct primes {
  12        struct rcu_head rcu;
  13        unsigned long last, sz;
  14        unsigned long primes[];
  15};
  16
  17#if BITS_PER_LONG == 64
  18static const struct primes small_primes = {
  19        .last = 61,
  20        .sz = 64,
  21        .primes = {
  22                BIT(2) |
  23                BIT(3) |
  24                BIT(5) |
  25                BIT(7) |
  26                BIT(11) |
  27                BIT(13) |
  28                BIT(17) |
  29                BIT(19) |
  30                BIT(23) |
  31                BIT(29) |
  32                BIT(31) |
  33                BIT(37) |
  34                BIT(41) |
  35                BIT(43) |
  36                BIT(47) |
  37                BIT(53) |
  38                BIT(59) |
  39                BIT(61)
  40        }
  41};
  42#elif BITS_PER_LONG == 32
  43static const struct primes small_primes = {
  44        .last = 31,
  45        .sz = 32,
  46        .primes = {
  47                BIT(2) |
  48                BIT(3) |
  49                BIT(5) |
  50                BIT(7) |
  51                BIT(11) |
  52                BIT(13) |
  53                BIT(17) |
  54                BIT(19) |
  55                BIT(23) |
  56                BIT(29) |
  57                BIT(31)
  58        }
  59};
  60#else
  61#error "unhandled BITS_PER_LONG"
  62#endif
  63
  64static DEFINE_MUTEX(lock);
  65static const struct primes __rcu *primes = RCU_INITIALIZER(&small_primes);
  66
  67static unsigned long selftest_max;
  68
  69static bool slow_is_prime_number(unsigned long x)
  70{
  71        unsigned long y = int_sqrt(x);
  72
  73        while (y > 1) {
  74                if ((x % y) == 0)
  75                        break;
  76                y--;
  77        }
  78
  79        return y == 1;
  80}
  81
  82static unsigned long slow_next_prime_number(unsigned long x)
  83{
  84        while (x < ULONG_MAX && !slow_is_prime_number(++x))
  85                ;
  86
  87        return x;
  88}
  89
  90static unsigned long clear_multiples(unsigned long x,
  91                                     unsigned long *p,
  92                                     unsigned long start,
  93                                     unsigned long end)
  94{
  95        unsigned long m;
  96
  97        m = 2 * x;
  98        if (m < start)
  99                m = roundup(start, x);
 100
 101        while (m < end) {
 102                __clear_bit(m, p);
 103                m += x;
 104        }
 105
 106        return x;
 107}
 108
 109static bool expand_to_next_prime(unsigned long x)
 110{
 111        const struct primes *p;
 112        struct primes *new;
 113        unsigned long sz, y;
 114
 115        /* Betrand's Postulate (or Chebyshev's theorem) states that if n > 3,
 116         * there is always at least one prime p between n and 2n - 2.
 117         * Equivalently, if n > 1, then there is always at least one prime p
 118         * such that n < p < 2n.
 119         *
 120         * http://mathworld.wolfram.com/BertrandsPostulate.html
 121         * https://en.wikipedia.org/wiki/Bertrand's_postulate
 122         */
 123        sz = 2 * x;
 124        if (sz < x)
 125                return false;
 126
 127        sz = round_up(sz, BITS_PER_LONG);
 128        new = kmalloc(sizeof(*new) + bitmap_size(sz),
 129                      GFP_KERNEL | __GFP_NOWARN);
 130        if (!new)
 131                return false;
 132
 133        mutex_lock(&lock);
 134        p = rcu_dereference_protected(primes, lockdep_is_held(&lock));
 135        if (x < p->last) {
 136                kfree(new);
 137                goto unlock;
 138        }
 139
 140        /* Where memory permits, track the primes using the
 141         * Sieve of Eratosthenes. The sieve is to remove all multiples of known
 142         * primes from the set, what remains in the set is therefore prime.
 143         */
 144        bitmap_fill(new->primes, sz);
 145        bitmap_copy(new->primes, p->primes, p->sz);
 146        for (y = 2UL; y < sz; y = find_next_bit(new->primes, sz, y + 1))
 147                new->last = clear_multiples(y, new->primes, p->sz, sz);
 148        new->sz = sz;
 149
 150        BUG_ON(new->last <= x);
 151
 152        rcu_assign_pointer(primes, new);
 153        if (p != &small_primes)
 154                kfree_rcu((struct primes *)p, rcu);
 155
 156unlock:
 157        mutex_unlock(&lock);
 158        return true;
 159}
 160
 161static void free_primes(void)
 162{
 163        const struct primes *p;
 164
 165        mutex_lock(&lock);
 166        p = rcu_dereference_protected(primes, lockdep_is_held(&lock));
 167        if (p != &small_primes) {
 168                rcu_assign_pointer(primes, &small_primes);
 169                kfree_rcu((struct primes *)p, rcu);
 170        }
 171        mutex_unlock(&lock);
 172}
 173
 174/**
 175 * next_prime_number - return the next prime number
 176 * @x: the starting point for searching to test
 177 *
 178 * A prime number is an integer greater than 1 that is only divisible by
 179 * itself and 1.  The set of prime numbers is computed using the Sieve of
 180 * Eratoshenes (on finding a prime, all multiples of that prime are removed
 181 * from the set) enabling a fast lookup of the next prime number larger than
 182 * @x. If the sieve fails (memory limitation), the search falls back to using
 183 * slow trial-divison, up to the value of ULONG_MAX (which is reported as the
 184 * final prime as a sentinel).
 185 *
 186 * Returns: the next prime number larger than @x
 187 */
 188unsigned long next_prime_number(unsigned long x)
 189{
 190        const struct primes *p;
 191
 192        rcu_read_lock();
 193        p = rcu_dereference(primes);
 194        while (x >= p->last) {
 195                rcu_read_unlock();
 196
 197                if (!expand_to_next_prime(x))
 198                        return slow_next_prime_number(x);
 199
 200                rcu_read_lock();
 201                p = rcu_dereference(primes);
 202        }
 203        x = find_next_bit(p->primes, p->last, x + 1);
 204        rcu_read_unlock();
 205
 206        return x;
 207}
 208EXPORT_SYMBOL(next_prime_number);
 209
 210/**
 211 * is_prime_number - test whether the given number is prime
 212 * @x: the number to test
 213 *
 214 * A prime number is an integer greater than 1 that is only divisible by
 215 * itself and 1. Internally a cache of prime numbers is kept (to speed up
 216 * searching for sequential primes, see next_prime_number()), but if the number
 217 * falls outside of that cache, its primality is tested using trial-divison.
 218 *
 219 * Returns: true if @x is prime, false for composite numbers.
 220 */
 221bool is_prime_number(unsigned long x)
 222{
 223        const struct primes *p;
 224        bool result;
 225
 226        rcu_read_lock();
 227        p = rcu_dereference(primes);
 228        while (x >= p->sz) {
 229                rcu_read_unlock();
 230
 231                if (!expand_to_next_prime(x))
 232                        return slow_is_prime_number(x);
 233
 234                rcu_read_lock();
 235                p = rcu_dereference(primes);
 236        }
 237        result = test_bit(x, p->primes);
 238        rcu_read_unlock();
 239
 240        return result;
 241}
 242EXPORT_SYMBOL(is_prime_number);
 243
 244static void dump_primes(void)
 245{
 246        const struct primes *p;
 247        char *buf;
 248
 249        buf = kmalloc(PAGE_SIZE, GFP_KERNEL);
 250
 251        rcu_read_lock();
 252        p = rcu_dereference(primes);
 253
 254        if (buf)
 255                bitmap_print_to_pagebuf(true, buf, p->primes, p->sz);
 256        pr_info("primes.{last=%lu, .sz=%lu, .primes[]=...x%lx} = %s\n",
 257                p->last, p->sz, p->primes[BITS_TO_LONGS(p->sz) - 1], buf);
 258
 259        rcu_read_unlock();
 260
 261        kfree(buf);
 262}
 263
 264static int selftest(unsigned long max)
 265{
 266        unsigned long x, last;
 267
 268        if (!max)
 269                return 0;
 270
 271        for (last = 0, x = 2; x < max; x++) {
 272                bool slow = slow_is_prime_number(x);
 273                bool fast = is_prime_number(x);
 274
 275                if (slow != fast) {
 276                        pr_err("inconsistent result for is-prime(%lu): slow=%s, fast=%s!\n",
 277                               x, slow ? "yes" : "no", fast ? "yes" : "no");
 278                        goto err;
 279                }
 280
 281                if (!slow)
 282                        continue;
 283
 284                if (next_prime_number(last) != x) {
 285                        pr_err("incorrect result for next-prime(%lu): expected %lu, got %lu\n",
 286                               last, x, next_prime_number(last));
 287                        goto err;
 288                }
 289                last = x;
 290        }
 291
 292        pr_info("%s(%lu) passed, last prime was %lu\n", __func__, x, last);
 293        return 0;
 294
 295err:
 296        dump_primes();
 297        return -EINVAL;
 298}
 299
 300static int __init primes_init(void)
 301{
 302        return selftest(selftest_max);
 303}
 304
 305static void __exit primes_exit(void)
 306{
 307        free_primes();
 308}
 309
 310module_init(primes_init);
 311module_exit(primes_exit);
 312
 313module_param_named(selftest, selftest_max, ulong, 0400);
 314
 315MODULE_AUTHOR("Intel Corporation");
 316MODULE_LICENSE("GPL");
 317