File Coverage

util.h
Criterion Covered Total %
statement 17 17 100.0
branch 8 12 66.6
condition n/a
subroutine n/a
pod n/a
total 25 29 86.2


line stmt bran cond sub pod time code
1             #ifndef MPU_UTIL_H
2             #define MPU_UTIL_H
3              
4             #include "ptypes.h"
5              
6             extern int _XS_get_verbose(void);
7             extern void _XS_set_verbose(int v);
8             extern int _XS_get_callgmp(void);
9             extern void _XS_set_callgmp(int v);
10             /* Disable all manual seeding */
11             extern int _XS_get_secure(void);
12             extern void _XS_set_secure(void);
13              
14             extern int is_prime(UV x);
15             extern UV next_prime(UV x);
16             extern UV prev_prime(UV x);
17              
18             extern void print_primes(UV low, UV high, int fd);
19              
20             extern int powerof(UV n);
21             extern int is_power(UV n, UV a);
22             extern UV rootof(UV n, UV k);
23             extern int primepower(UV n, UV* prime);
24             extern UV valuation(UV n, UV k);
25             extern UV logint(UV n, UV b);
26             extern UV mpu_popcount_string(const char* ptr, uint32_t len);
27              
28             extern signed char* _moebius_range(UV low, UV high);
29             extern UV* _totient_range(UV low, UV high);
30             extern IV mertens(UV n);
31             extern long double chebyshev_function(UV n, int which); /* 0 = theta, 1 = psi */
32              
33             extern long double Ei(long double x);
34             extern long double Li(long double x);
35             extern long double ld_riemann_zeta(long double x);
36             extern long double RiemannR(long double x);
37             extern long double lambertw(long double k);
38             extern UV inverse_li(UV x);
39             extern UV inverse_R(UV x);
40              
41             extern int kronecker_uu(UV a, UV b);
42             extern int kronecker_su(IV a, UV b);
43             extern int kronecker_ss(IV a, IV b);
44              
45             extern UV factorial(UV n);
46             extern UV binomial(UV n, UV k);
47             extern IV gcdext(IV a, IV b, IV* u, IV* v, IV* s, IV* t); /* Ext Euclidean */
48             extern UV modinverse(UV a, UV p); /* Returns 1/a mod p */
49             extern UV divmod(UV a, UV b, UV n); /* Returns a/b mod n */
50             extern int sqrtmod(UV* s, UV a, UV p); /* sqrt(a) mod p */
51             extern int sqrtmod_composite(UV* s, UV a,UV n);/* sqrt(a) mod n */
52             extern UV chinese(UV* a, UV* n, UV num, int *status); /* Chinese Remainder */
53              
54             extern UV totient(UV n);
55             extern int moebius(UV n);
56             extern UV exp_mangoldt(UV n);
57             extern UV carmichael_lambda(UV n);
58             extern UV jordan_totient(UV k, UV n);
59             extern UV znprimroot(UV n);
60             extern UV znorder(UV a, UV n);
61             extern int is_primitive_root(UV a, UV n, int nprime);
62             extern UV factorialmod(UV n, UV m);
63             #define is_square_free(n) (moebius(n) != 0)
64             extern int is_fundamental(UV n, int neg);
65             extern int is_totient(UV n);
66             extern int is_semiprime(UV n);
67             extern int is_carmichael(UV n);
68             extern UV is_quasi_carmichael(UV n); /* Returns number of bases */
69             extern UV pillai_v(UV n); /* v: v! % n == n-1 && n % v != 1 */
70              
71             extern UV stirling3(UV n, UV m);
72             extern IV stirling2(UV n, UV m);
73             extern IV stirling1(UV n, UV m);
74              
75             extern IV hclassno(UV n);
76             extern IV ramanujan_tau(UV n);
77              
78             extern char* pidigits(int digits);
79              
80             extern int strnum_minmax(int min, char* a, STRLEN alen, char* b, STRLEN blen);
81              
82             extern int from_digit_string(UV* n, const char* s, int base);
83             extern int from_digit_to_UV(UV* rn, UV* r, int len, int base);
84             extern int from_digit_to_str(char** rstr, UV* r, int len, int base);
85             extern int to_digit_array(int* bits, UV n, int base, int length);
86             extern int to_digit_string(char *s, UV n, int base, int length);
87             extern int to_string_128(char s[40], IV hi, UV lo);
88              
89             extern int is_catalan_pseudoprime(UV n);
90              
91             extern UV polygonal_root(UV n, UV k, int* overflow);
92              
93             extern int num_to_perm(UV rank, int n, int *vec);
94             extern int perm_to_num(int n, int *vec, UV *rank);
95             extern void randperm(void* ctx, UV n, UV k, UV *S);
96              
97             extern UV gcdz(UV x, UV y);
98              
99             #if defined(FUNC_isqrt) || defined(FUNC_is_perfect_square)
100 20           static UV isqrt(UV n) {
101             UV root;
102             #if BITS_PER_WORD == 32
103             if (n >= UVCONST(4294836225)) return UVCONST(65535);
104             #else
105 20 50         if (n >= UVCONST(18446744065119617025)) return UVCONST(4294967295);
106             #endif
107 20           root = (UV) sqrt((double)n);
108 20 50         while (root*root > n) root--;
109 20 50         while ((root+1)*(root+1) <= n) root++;
110 20           return root;
111             }
112             #endif
113              
114             #if defined(FUNC_icbrt) || defined(FUNC_is_perfect_cube)
115 20           static UV icbrt(UV n) {
116 20           UV b, root = 0;
117             #if BITS_PER_WORD == 32
118             int s = 30;
119             if (n >= UVCONST(4291015625)) return UVCONST(1625);
120             #else
121 20           int s = 63;
122 20 50         if (n >= UVCONST(18446724184312856125)) return UVCONST(2642245);
123             #endif
124 460 100         for ( ; s >= 0; s -= 3) {
125 440           root += root;
126 440           b = 3*root*(root+1)+1;
127 440 100         if ((n >> s) >= b) {
128 124           n -= b << s;
129 124           root++;
130             }
131             }
132 20           return root;
133             }
134             #endif
135              
136             #if defined(FUNC_ipow)
137             static UV ipow(UV n, UV k) {
138             UV p = 1;
139             while (k) {
140             if (k & 1) p *= n;
141             k >>= 1;
142             if (k) n *= n;
143             }
144             return p;
145             }
146             #endif
147              
148             #if defined(FUNC_gcd_ui) || defined(FUNC_lcm_ui)
149             /* If we have a very fast ctz, then use the fast FLINT version of gcd */
150             #if defined(__GNUC__) && (__GNUC__ >= 4 || (__GNUC__ == 3 && __GNUC_MINOR__ >= 4))
151             #define gcd_ui(x,y) gcdz(x,y)
152             #else
153             static UV gcd_ui(UV x, UV y) {
154             UV t;
155             if (y < x) { t = x; x = y; y = t; }
156             while (y > 0) {
157             t = y; y = x % y; x = t; /* y1 <- x0 % y0 ; x1 <- y0 */
158             }
159             return x;
160             }
161             #endif
162             #endif
163              
164             #ifdef FUNC_lcm_ui
165             static UV lcm_ui(UV x, UV y) {
166             /* Can overflow if lcm(x,y) > 2^64 (e.g. two primes each > 2^32) */
167             return x * (y / gcd_ui(x,y));
168             }
169             #endif
170              
171             #ifdef FUNC_is_perfect_square
172             /* See: http://mersenneforum.org/showpost.php?p=110896 */
173             static int is_perfect_square(UV n)
174             {
175             /* Step 1, reduce to 18% of inputs */
176             uint32_t m = n & 127;
177             if ((m*0x8bc40d7d) & (m*0xa1e2f5d1) & 0x14020a) return 0;
178             /* Step 2, reduce to 7% of inputs (mod 99 reduces to 4% but slower) */
179             m = n %240; if ((m*0xfa445556) & (m*0x8021feb1) & 0x614aaa0f) return 0;
180             /* m = n % 99; if ((m*0x5411171d) & (m*0xe41dd1c7) & 0x80028a80) return 0; */
181             /* Step 3, do the square root instead of any more rejections */
182             m = isqrt(n);
183             return (UV)m*(UV)m == n;
184             }
185             #endif
186              
187             #ifdef FUNC_is_perfect_cube
188             static int is_perfect_cube(UV n)
189             {
190             UV m;
191             if ((n & 3) == 2) return 0;
192             /* m = n & 511; if ((m*5016427) & (m*95638165) & 438) return 0; */
193             m = n % 117; if ((m*833230740) & (m*120676722) & 813764715) return 0;
194             m = n % 133; if ((m*76846229) & (m*305817297) & 306336544) return 0;
195             m = icbrt(n);
196             return m*m*m == n;
197             }
198             #endif
199              
200             #ifdef FUNC_is_perfect_fifth
201             static int is_perfect_fifth(UV n)
202             {
203             UV m;
204             if ((n & 3) == 2) return 0;
205             m = n % 88; if ((m*85413603) & (m*76260301) & 26476550) return 0;
206             m = n % 31; if ((m*80682551) & (m*73523539) & 45414528) return 0;
207             m = n % 41; if ((m*92806493) & (m*130690042) & 35668129) return 0;
208             /* m = n % 25; if ((m*109794298) & (m*105535723) & 16097553) return 0; */
209             m = rootof(n, 5);
210             return m*m*m*m*m == n;
211             }
212             #endif
213              
214             #ifdef FUNC_is_perfect_seventh
215             static int is_perfect_seventh(UV n)
216             {
217             UV m;
218             /* if ((n & 3) == 2) return 0; */
219             m = n & 511; if ((m*97259473) & (m*51311663) & 894) return 0;
220             m = n % 49; if ((m*109645301) & (m*76482737) & 593520192) return 0;
221             m = n % 71; if ((m*71818386) & (m*38821587) & 35299393) return 0;
222             /* m = n % 43; if ((m*101368253) & (m*814158665) & 142131408) return 0; */
223             /* m = n % 29; if ((m*81935611) & (m*84736134) & 37831965) return 0; */
224             /* m = n % 116; if ((m*348163737) & (m*1539055705) & 2735997248) return 0; */
225             m = rootof(n, 7);
226             return m*m*m*m*m*m*m == n;
227             }
228             #endif
229              
230             #if defined(FUNC_clz) || defined(FUNC_ctz) || defined(FUNC_log2floor)
231             /* log2floor(n) gives the location of the first set bit (starting from left)
232             * ctz(n) gives the number of times n is divisible by 2
233             * clz(n) gives the number of zeros on the left */
234             #if defined(__GNUC__) && (__GNUC__ >= 4 || (__GNUC__ == 3 && __GNUC_MINOR__ >= 4))
235             #if BITS_PER_WORD == 64
236             #define ctz(n) ((n) ? __builtin_ctzll(n) : 64)
237             #define clz(n) ((n) ? __builtin_clzll(n) : 64)
238             #define log2floor(n) ((n) ? 63-__builtin_clzll(n) : 0)
239             #else
240             #define ctz(n) ((n) ? __builtin_ctzl(n) : 32)
241             #define clz(n) ((n) ? __builtin_clzl(n) : 32)
242             #define log2floor(n) ((n) ? 31-__builtin_clzl(n) : 0)
243             #endif
244              
245             /* For MSC, we need to use _BitScanForward and _BitScanReverse. The way to
246             * get to them has changed, so we're going to only use them on new systems.
247             * The performance of these functions are not super critical.
248             * What is: popcnt, mulmod, and muladd.
249             */
250             #elif defined (_MSC_VER) && _MSC_VER >= 1400 && !defined(__clang__) && !defined(_WIN32_WCE)
251             #include
252             #ifdef FUNC_ctz
253             static int ctz(UV n) {
254             UV tz = 0;
255             #if BITS_PER_WORD == 64
256             if (_BitScanForward64(&tz, n)) return tz; else return 64;
257             #else
258             if (_BitScanForward(&tz, n)) return tz; else return 32;
259             #endif
260             }
261             #endif
262             #if defined(FUNC_clz) || defined(FUNC_log2floor)
263             static int log2floor(UV n) {
264             UV lz = 0;
265             #if BITS_PER_WORD == 64
266             if (_BitScanReverse64(&lz, n)) return lz; else return 0;
267             #else
268             if (_BitScanReverse(&lz, n)) return lz; else return 0;
269             #endif
270             }
271             #endif
272             #elif BITS_PER_WORD == 64
273             static const unsigned char _debruijn64[64] = {
274             63, 0,58, 1,59,47,53, 2, 60,39,48,27,54,33,42, 3, 61,51,37,40,49,18,28,20,
275             55,30,34,11,43,14,22, 4, 62,57,46,52,38,26,32,41, 50,36,17,19,29,10,13,21,
276             56,45,25,31,35,16, 9,12, 44,24,15, 8,23, 7, 6, 5 };
277             #ifdef FUNC_ctz
278             static unsigned int ctz(UV n) {
279             return n ? _debruijn64[((n & -n)*UVCONST(0x07EDD5E59A4E28C2)) >> 58] : 64;
280             }
281             #endif
282             #if defined(FUNC_clz) || defined(FUNC_log2floor)
283             static unsigned int log2floor(UV n) {
284             if (n == 0) return 0;
285             n |= n >> 1; n |= n >> 2; n |= n >> 4;
286             n |= n >> 8; n |= n >> 16; n |= n >> 32;
287             return _debruijn64[((n-(n>>1))*UVCONST(0x07EDD5E59A4E28C2)) >> 58];
288             }
289             #endif
290             #else
291             #ifdef FUNC_ctz
292             static const unsigned char _trail_debruijn32[32] = {
293             0, 1,28, 2,29,14,24, 3,30,22,20,15,25,17, 4, 8,
294             31,27,13,23,21,19,16, 7,26,12,18, 6,11, 5,10, 9 };
295             static unsigned int ctz(UV n) {
296             return n ? _trail_debruijn32[((n & -n) * UVCONST(0x077CB531)) >> 27] : 32;
297             }
298             #endif
299             #if defined(FUNC_clz) || defined(FUNC_log2floor)
300             static const unsigned char _lead_debruijn32[32] = {
301             0, 9, 1, 10, 13, 21, 2, 29, 11, 14, 16, 18, 22, 25, 3, 30,
302             8, 12, 20, 28, 15, 17, 24, 7, 19, 27, 23, 6, 26, 5, 4, 31 };
303             static unsigned int log2floor(UV n) {
304             if (n == 0) return 0;
305             n |= n >> 1; n |= n >> 2; n |= n >> 4; n |= n >> 8; n |= n >> 16;
306             return _lead_debruijn32[(n * UVCONST(0x07C4ACDD)) >> 27];
307             }
308             #endif
309             #endif
310             #if defined(FUNC_clz) && !defined(clz)
311             #define clz(n) ( (n) ? BITS_PER_WORD-1-log2floor(n) : BITS_PER_WORD )
312             #endif
313             #endif /* End of log2floor, clz, and ctz */
314              
315             #ifdef FUNC_popcnt
316             /* GCC 3.4 - 4.1 has broken 64-bit popcount.
317             * GCC 4.2+ can generate awful code when it doesn't have asm (GCC bug 36041).
318             * When the asm is present (e.g. compile with -march=native on a platform that
319             * has them, like Nahelem+), then it is almost as fast as manually written asm. */
320             #if BITS_PER_WORD == 64
321             #if defined(__POPCNT__) && defined(__GNUC__) && (__GNUC__> 4 || (__GNUC__== 4 && __GNUC_MINOR__> 1))
322             #define popcnt(b) __builtin_popcountll(b)
323             #else
324             static UV popcnt(UV b) {
325             b -= (b >> 1) & 0x5555555555555555;
326             b = (b & 0x3333333333333333) + ((b >> 2) & 0x3333333333333333);
327             b = (b + (b >> 4)) & 0x0f0f0f0f0f0f0f0f;
328             return (b * 0x0101010101010101) >> 56;
329             }
330             #endif
331             #else
332             static UV popcnt(UV b) {
333             b -= (b >> 1) & 0x55555555;
334             b = (b & 0x33333333) + ((b >> 2) & 0x33333333);
335             b = (b + (b >> 4)) & 0x0f0f0f0f;
336             return (b * 0x01010101) >> 24;
337             }
338             #endif
339             #endif
340              
341             #endif