/src/gmp-6.2.1/mpz/primorial_ui.c
Line | Count | Source |
1 | | /* mpz_primorial_ui(RESULT, N) -- Set RESULT to N# the product of primes <= N. |
2 | | |
3 | | Contributed to the GNU project by Marco Bodrato. |
4 | | |
5 | | Copyright 2012, 2015, 2016 Free Software Foundation, Inc. |
6 | | |
7 | | This file is part of the GNU MP Library. |
8 | | |
9 | | The GNU MP Library is free software; you can redistribute it and/or modify |
10 | | it under the terms of either: |
11 | | |
12 | | * the GNU Lesser General Public License as published by the Free |
13 | | Software Foundation; either version 3 of the License, or (at your |
14 | | option) any later version. |
15 | | |
16 | | or |
17 | | |
18 | | * the GNU General Public License as published by the Free Software |
19 | | Foundation; either version 2 of the License, or (at your option) any |
20 | | later version. |
21 | | |
22 | | or both in parallel, as here. |
23 | | |
24 | | The GNU MP Library is distributed in the hope that it will be useful, but |
25 | | WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
26 | | or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
27 | | for more details. |
28 | | |
29 | | You should have received copies of the GNU General Public License and the |
30 | | GNU Lesser General Public License along with the GNU MP Library. If not, |
31 | | see https://www.gnu.org/licenses/. */ |
32 | | |
33 | | #include "gmp-impl.h" |
34 | | |
35 | | /* TODO: Remove duplicated constants / macros / static functions... |
36 | | */ |
37 | | |
38 | | /*************************************************************/ |
39 | | /* Section macros: common macros, for swing/fac/bin (&sieve) */ |
40 | | /*************************************************************/ |
41 | | |
42 | | #define FACTOR_LIST_STORE(P, PR, MAX_PR, VEC, I) \ |
43 | 116 | do { \ |
44 | 116 | if ((PR) > (MAX_PR)) { \ |
45 | 5 | (VEC)[(I)++] = (PR); \ |
46 | 5 | (PR) = (P); \ |
47 | 5 | } else \ |
48 | 116 | (PR) *= (P); \ |
49 | 116 | } while (0) |
50 | | |
51 | | #define LOOP_ON_SIEVE_CONTINUE(prime,end,sieve) \ |
52 | 7 | __max_i = (end); \ |
53 | 7 | \ |
54 | 155 | do { \ |
55 | 155 | ++__i; \ |
56 | 155 | if (((sieve)[__index] & __mask) == 0) \ |
57 | 155 | { \ |
58 | 116 | mp_limb_t prime; \ |
59 | 116 | prime = id_to_n(__i) |
60 | | |
61 | | #define LOOP_ON_SIEVE_BEGIN(prime,start,end,off,sieve) \ |
62 | 7 | do { \ |
63 | 7 | mp_limb_t __mask, __index, __max_i, __i; \ |
64 | 7 | \ |
65 | 7 | __i = (start)-(off); \ |
66 | 7 | __index = __i / GMP_LIMB_BITS; \ |
67 | 7 | __mask = CNST_LIMB(1) << (__i % GMP_LIMB_BITS); \ |
68 | 7 | __i += (off); \ |
69 | 7 | \ |
70 | 7 | LOOP_ON_SIEVE_CONTINUE(prime,end,sieve) |
71 | | |
72 | | #define LOOP_ON_SIEVE_STOP \ |
73 | 116 | } \ |
74 | 155 | __mask = __mask << 1 | __mask >> (GMP_LIMB_BITS-1); \ |
75 | 155 | __index += __mask & 1; \ |
76 | 155 | } while (__i <= __max_i) |
77 | | |
78 | | #define LOOP_ON_SIEVE_END \ |
79 | 116 | LOOP_ON_SIEVE_STOP; \ |
80 | 7 | } while (0) |
81 | | |
82 | | /*********************************************************/ |
83 | | /* Section sieve: sieving functions and tools for primes */ |
84 | | /*********************************************************/ |
85 | | |
86 | | #if 0 |
87 | | static mp_limb_t |
88 | | bit_to_n (mp_limb_t bit) { return (bit*3+4)|1; } |
89 | | #endif |
90 | | |
91 | | /* id_to_n (x) = bit_to_n (x-1) = (id*3+1)|1*/ |
92 | | static mp_limb_t |
93 | 116 | id_to_n (mp_limb_t id) { return id*3+1+(id&1); } |
94 | | |
95 | | /* n_to_bit (n) = ((n-1)&(-CNST_LIMB(2)))/3U-1 */ |
96 | | static mp_limb_t |
97 | 21 | n_to_bit (mp_limb_t n) { return ((n-5)|1)/3U; } |
98 | | |
99 | | #if WANT_ASSERT |
100 | | static mp_size_t |
101 | 7 | primesieve_size (mp_limb_t n) { return n_to_bit(n) / GMP_LIMB_BITS + 1; } |
102 | | #endif |
103 | | |
104 | | /*********************************************************/ |
105 | | /* Section primorial: implementation */ |
106 | | /*********************************************************/ |
107 | | |
108 | | void |
109 | | mpz_primorial_ui (mpz_ptr x, unsigned long n) |
110 | 8 | { |
111 | 8 | ASSERT (n <= GMP_NUMB_MAX); |
112 | | |
113 | 8 | if (n < 5) |
114 | 1 | { |
115 | 1 | MPZ_NEWALLOC (x, 1)[0] = (066211 >> (n*3)) & 7; |
116 | 1 | SIZ (x) = 1; |
117 | 1 | } |
118 | 7 | else |
119 | 7 | { |
120 | 7 | mp_limb_t *sieve, *factors; |
121 | 7 | mp_size_t size, j; |
122 | 7 | mp_limb_t prod; |
123 | 7 | TMP_DECL; |
124 | | |
125 | 7 | size = n / GMP_NUMB_BITS; |
126 | 7 | size = size + (size >> 1) + 1; |
127 | 7 | ASSERT (size >= primesieve_size (n)); |
128 | 7 | sieve = MPZ_NEWALLOC (x, size); |
129 | 7 | size = (gmp_primesieve (sieve, n) + 1) / log_n_max (n) + 1; |
130 | | |
131 | 7 | TMP_MARK; |
132 | 7 | factors = TMP_ALLOC_LIMBS (size); |
133 | | |
134 | 7 | j = 0; |
135 | | |
136 | 7 | prod = 6; |
137 | | |
138 | | /* Store primes from 5 to n */ |
139 | 7 | { |
140 | 7 | mp_limb_t max_prod; |
141 | | |
142 | 7 | max_prod = GMP_NUMB_MAX / n; |
143 | | |
144 | 278 | LOOP_ON_SIEVE_BEGIN (prime, n_to_bit(5), n_to_bit (n), 0, sieve); |
145 | 278 | FACTOR_LIST_STORE (prime, prod, max_prod, factors, j); |
146 | 278 | LOOP_ON_SIEVE_END; |
147 | 7 | } |
148 | | |
149 | 7 | if (j != 0) |
150 | 5 | { |
151 | 5 | factors[j++] = prod; |
152 | 5 | mpz_prodlimbs (x, factors, j); |
153 | 5 | } |
154 | 2 | else |
155 | 2 | { |
156 | 2 | PTR (x)[0] = prod; |
157 | 2 | SIZ (x) = 1; |
158 | 2 | } |
159 | | |
160 | 7 | TMP_FREE; |
161 | 7 | } |
162 | 8 | } |