1 /*
2 * Copyright © 2016 Intel Corporation
3 *
4 * Permission is hereby granted, free of charge, to any person obtaining a
5 * copy of this software and associated documentation files (the "Software"),
6 * to deal in the Software without restriction, including without limitation
7 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
8 * and/or sell copies of the Software, and to permit persons to whom the
9 * Software is furnished to do so, subject to the following conditions:
10 *
11 * The above copyright notice and this permission notice (including the next
12 * paragraph) shall be included in all copies or substantial portions of the
13 * Software.
14 *
15 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
16 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
17 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
18 * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
19 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
20 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
21 * IN THE SOFTWARE.
22 */
23
24 #include "igt_primes.h"
25
26 #include <stdlib.h>
27 #include <stdbool.h>
28 #include <string.h>
29 #include <math.h>
30
31 /**
32 * SECTION:igt_primes
33 * @short_description: Prime numbers helper library
34 * @title: Primes
35 * @include: igt_primes.h
36 */
37
38 #define BITS_PER_CHAR 8
39 #define BITS_PER_LONG (sizeof(long)*BITS_PER_CHAR)
40
41 #define BITMAP_FIRST_WORD_MASK(start) (~0UL << ((start) & (BITS_PER_LONG - 1)))
42 #define BITMAP_LAST_WORD_MASK(nbits) (~0UL >> (-(nbits) & (BITS_PER_LONG - 1)))
43
44 #define __round_mask(x, y) ((__typeof__(x))((y)-1))
45 #define round_up(x, y) ((((x)-1) | __round_mask(x, y))+1)
46 #define round_down(x, y) ((x) & ~__round_mask(x, y))
47
48 #define min(x, y) ({ \
49 typeof(x) _min1 = (x); \
50 typeof(y) _min2 = (y); \
51 (void) (&_min1 == &_min2); \
52 _min1 < _min2 ? _min1 : _min2; \
53 })
54
55 #define max(x, y) ({ \
56 typeof(x) _max1 = (x); \
57 typeof(y) _max2 = (y); \
58 (void) (&_max1 == &_max2); \
59 _max1 > _max2 ? _max1 : _max2; \
60 })
61
__bit__(unsigned long nr)62 static inline unsigned long __bit__(unsigned long nr)
63 {
64 return 1UL << (nr % BITS_PER_LONG);
65 }
66
set_bit(unsigned long nr,unsigned long * addr)67 static inline void set_bit(unsigned long nr, unsigned long *addr)
68 {
69 addr[nr / BITS_PER_LONG] |= __bit__(nr);
70 }
71
clear_bit(unsigned long nr,unsigned long * addr)72 static inline void clear_bit(unsigned long nr, unsigned long *addr)
73 {
74 addr[nr / BITS_PER_LONG] &= ~__bit__(nr);
75 }
76
test_bit(unsigned long nr,const unsigned long * addr)77 static inline bool test_bit(unsigned long nr, const unsigned long *addr)
78 {
79 return addr[nr / BITS_PER_LONG] & __bit__(nr);
80 }
81
82 static unsigned long
__find_next_bit(const unsigned long * addr,unsigned long nbits,unsigned long start,unsigned long invert)83 __find_next_bit(const unsigned long *addr,
84 unsigned long nbits, unsigned long start,
85 unsigned long invert)
86 {
87 unsigned long tmp;
88
89 if (!nbits || start >= nbits)
90 return nbits;
91
92 tmp = addr[start / BITS_PER_LONG] ^ invert;
93
94 /* Handle 1st word. */
95 tmp &= BITMAP_FIRST_WORD_MASK(start);
96 start = round_down(start, BITS_PER_LONG);
97
98 while (!tmp) {
99 start += BITS_PER_LONG;
100 if (start >= nbits)
101 return nbits;
102
103 tmp = addr[start / BITS_PER_LONG] ^ invert;
104 }
105
106 return min(start + __builtin_ffsl(tmp) - 1, nbits);
107 }
108
find_next_bit(const unsigned long * addr,unsigned long size,unsigned long offset)109 static unsigned long find_next_bit(const unsigned long *addr,
110 unsigned long size,
111 unsigned long offset)
112 {
113 return __find_next_bit(addr, size, offset, 0UL);
114 }
115
slow_next_prime_number(unsigned long x)116 static unsigned long slow_next_prime_number(unsigned long x)
117 {
118 for (;;) {
119 unsigned long y = sqrt(++x) + 1;
120 while (y > 1) {
121 if ((x % y) == 0)
122 break;
123 y--;
124 }
125 if (y == 1)
126 return x;
127 }
128 }
129
mark_multiples(unsigned long x,unsigned long * primes,unsigned long start,unsigned long end)130 static unsigned long mark_multiples(unsigned long x,
131 unsigned long *primes,
132 unsigned long start,
133 unsigned long end)
134 {
135 unsigned long m;
136
137 m = 2*x;
138 if (m < start)
139 m = (start / x + 1) * x;
140
141 while (m < end) {
142 clear_bit(m, primes);
143 m += x;
144 }
145
146 return x;
147 }
148
igt_next_prime_number(unsigned long x)149 unsigned long igt_next_prime_number(unsigned long x)
150 {
151 static unsigned long *primes;
152 static unsigned long last, last_sz;
153
154 if (x == 0)
155 return 1; /* a white lie for for_each_prime_number() */
156 if (x == 1)
157 return 2;
158
159 if (x >= last) {
160 unsigned long sz, y;
161 unsigned long *nprimes;
162
163 sz = x*x;
164 if (sz < x)
165 return slow_next_prime_number(x);
166
167 sz = round_up(sz, BITS_PER_LONG);
168 nprimes = realloc(primes, sz / sizeof(long));
169 if (!nprimes)
170 return slow_next_prime_number(x);
171
172 /* Where memory permits, track the primes using the
173 * Sieve of Eratosthenes.
174 */
175 memset(nprimes + last_sz / BITS_PER_LONG,
176 0xff, (sz - last_sz) / sizeof(long));
177 for (y = 2UL; y < sz; y = find_next_bit(nprimes, sz, y + 1))
178 last = mark_multiples(y, nprimes, last_sz, sz);
179
180 primes = nprimes;
181 last_sz = sz;
182 }
183
184 return find_next_bit(primes, last, x + 1);
185 }
186