1 /* gen-trialdivtab.c
2
3 Contributed to the GNU project by Torbjorn Granlund.
4
5 Copyright 2009, 2012, 2013 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 /*
34 Generate tables for fast, division-free trial division for GMP.
35
36 There is one main table, ptab. It contains primes, multiplied together, and
37 several types of pre-computed inverses. It refers to tables of the type
38 dtab, via the last two indices. That table contains the individual primes in
39 the range, except that the primes are not actually included in the table (see
40 the P macro; it sneakingly excludes the primes themselves). Instead, the
41 dtab tables contains tuples for each prime (modular-inverse, limit) used for
42 divisibility checks.
43
44 This interface is not intended for division of very many primes, since then
45 other algorithms apply.
46 */
47
48 #include <stdlib.h>
49 #include <stdio.h>
50 #include "bootstrap.c"
51
52 int sumspills (mpz_t, mpz_t *, int);
53 void mpn_mod_1s_4p_cps (mpz_t [7], mpz_t);
54
55 int limb_bits;
56
57 mpz_t B;
58
59 int
main(int argc,char * argv[])60 main (int argc, char *argv[])
61 {
62 unsigned long t, p;
63 mpz_t ppp, acc, inv, gmp_numb_max, tmp, Bhalf;
64 mpz_t pre[7];
65 int i;
66 int start_p, end_p, interval_start, interval_end, omitted_p;
67 const char *endtok;
68 int stop;
69 int np, start_idx;
70
71 if (argc < 2)
72 {
73 fprintf (stderr, "usage: %s bits endprime\n", argv[0]);
74 exit (1);
75 }
76
77 limb_bits = atoi (argv[1]);
78
79 end_p = 1290; /* default end prime */
80 if (argc == 3)
81 end_p = atoi (argv[2]);
82
83 printf ("#if GMP_LIMB_BITS != %d\n", limb_bits);
84 printf ("#error This table is for GMP_LIMB_BITS = %d\n", limb_bits);
85 printf ("#endif\n\n");
86
87 printf ("#if GMP_NAIL_BITS != 0\n");
88 printf ("#error This table does not support nails\n");
89 printf ("#endif\n\n");
90
91 for (i = 0; i < 7; i++)
92 mpz_init (pre[i]);
93
94 mpz_init_set_ui (gmp_numb_max, 1);
95 mpz_mul_2exp (gmp_numb_max, gmp_numb_max, limb_bits);
96 mpz_sub_ui (gmp_numb_max, gmp_numb_max, 1);
97
98 mpz_init (tmp);
99 mpz_init (inv);
100
101 mpz_init_set_ui (B, 1); mpz_mul_2exp (B, B, limb_bits);
102 mpz_init_set_ui (Bhalf, 1); mpz_mul_2exp (Bhalf, Bhalf, limb_bits - 1);
103
104 start_p = 3;
105
106 mpz_init_set_ui (ppp, 1);
107 mpz_init (acc);
108 interval_start = start_p;
109 omitted_p = 3;
110 interval_end = 0;
111
112 /* printf ("static struct gmp_primes_dtab gmp_primes_dtab[] = {\n"); */
113
114 printf ("#ifdef WANT_dtab\n");
115
116 for (t = start_p; t <= end_p; t += 2)
117 {
118 if (! isprime (t))
119 continue;
120
121 mpz_mul_ui (acc, ppp, t);
122 stop = mpz_cmp (acc, Bhalf) >= 0;
123 if (!stop)
124 {
125 mpn_mod_1s_4p_cps (pre, acc);
126 stop = sumspills (acc, pre + 2, 5);
127 }
128
129 if (stop)
130 {
131 for (p = interval_start; p <= interval_end; p += 2)
132 {
133 if (! isprime (p))
134 continue;
135
136 printf (" P(%d,", (int) p);
137 mpz_invert_ui_2exp (inv, p, limb_bits);
138 printf ("CNST_LIMB(0x"); mpz_out_str (stdout, 16, inv); printf ("),");
139
140 mpz_tdiv_q_ui (tmp, gmp_numb_max, p);
141 printf ("CNST_LIMB(0x"); mpz_out_str (stdout, 16, tmp);
142 printf (")),\n");
143 }
144 mpz_set_ui (ppp, t);
145 interval_start = t;
146 omitted_p = t;
147 }
148 else
149 {
150 mpz_set (ppp, acc);
151 }
152 interval_end = t;
153 }
154 printf ("#define SMALLEST_OMITTED_PRIME %d\n", (int) omitted_p);
155 printf ("#endif\n");
156
157 printf ("#ifdef WANT_ptab\n");
158
159 /* printf ("static struct gmp_primes_ptab gmp_primes_ptab[] = {\n"); */
160
161 endtok = "";
162
163 mpz_set_ui (ppp, 1);
164 interval_start = start_p;
165 interval_end = 0;
166 np = 0;
167 start_idx = 0;
168 for (t = start_p; t <= end_p; t += 2)
169 {
170 if (! isprime (t))
171 continue;
172
173 mpz_mul_ui (acc, ppp, t);
174
175 stop = mpz_cmp (acc, Bhalf) >= 0;
176 if (!stop)
177 {
178 mpn_mod_1s_4p_cps (pre, acc);
179 stop = sumspills (acc, pre + 2, 5);
180 }
181
182 if (stop)
183 {
184 mpn_mod_1s_4p_cps (pre, ppp);
185 printf ("%s", endtok);
186 printf (" {CNST_LIMB(0x"); mpz_out_str (stdout, 16, ppp);
187 printf ("),{CNST_LIMB(0x"); mpz_out_str (stdout, 16, pre[0]);
188 printf ("),%d", (int) PTR(pre[1])[0]);
189 for (i = 0; i < 5; i++)
190 {
191 printf (",");
192 printf ("CNST_LIMB(0x"); mpz_out_str (stdout, 16, pre[2 + i]);
193 printf (")");
194 }
195 printf ("},");
196 printf ("%d,", start_idx);
197 printf ("%d}", np - start_idx);
198
199 endtok = ",\n";
200 mpz_set_ui (ppp, t);
201 interval_start = t;
202 start_idx = np;
203 }
204 else
205 {
206 mpz_set (ppp, acc);
207 }
208 interval_end = t;
209 np++;
210 }
211
212 printf ("\n");
213 printf ("#endif\n");
214
215 return 0;
216 }
217
218 unsigned long
mpz_log2(mpz_t x)219 mpz_log2 (mpz_t x)
220 {
221 return mpz_sgn (x) ? mpz_sizeinbase (x, 2) : 0;
222 }
223
224 void
mpn_mod_1s_4p_cps(mpz_t cps[7],mpz_t bparm)225 mpn_mod_1s_4p_cps (mpz_t cps[7], mpz_t bparm)
226 {
227 mpz_t b, bi;
228 mpz_t B1modb, B2modb, B3modb, B4modb, B5modb;
229 mpz_t t;
230 int cnt;
231
232 mpz_init_set (b, bparm);
233
234 cnt = limb_bits - mpz_log2 (b);
235
236 mpz_init (bi);
237 mpz_init (t);
238 mpz_init (B1modb);
239 mpz_init (B2modb);
240 mpz_init (B3modb);
241 mpz_init (B4modb);
242 mpz_init (B5modb);
243
244 mpz_set_ui (t, 1);
245 mpz_mul_2exp (t, t, limb_bits - cnt);
246 mpz_sub (t, t, b);
247 mpz_mul_2exp (t, t, limb_bits);
248 mpz_tdiv_q (bi, t, b); /* bi = B^2/b, except msb */
249
250 mpz_set_ui (t, 1);
251 mpz_mul_2exp (t, t, limb_bits); /* t = B */
252 mpz_tdiv_r (B1modb, t, b);
253
254 mpz_mul_2exp (t, B1modb, limb_bits);
255 mpz_tdiv_r (B2modb, t, b);
256
257 mpz_mul_2exp (t, B2modb, limb_bits);
258 mpz_tdiv_r (B3modb, t, b);
259
260 mpz_mul_2exp (t, B3modb, limb_bits);
261 mpz_tdiv_r (B4modb, t, b);
262
263 mpz_mul_2exp (t, B4modb, limb_bits);
264 mpz_tdiv_r (B5modb, t, b);
265
266 mpz_set (cps[0], bi);
267 mpz_set_ui (cps[1], cnt);
268 mpz_tdiv_q_2exp (cps[2], B1modb, 0);
269 mpz_tdiv_q_2exp (cps[3], B2modb, 0);
270 mpz_tdiv_q_2exp (cps[4], B3modb, 0);
271 mpz_tdiv_q_2exp (cps[5], B4modb, 0);
272 mpz_tdiv_q_2exp (cps[6], B5modb, 0);
273
274 mpz_clear (b);
275 mpz_clear (bi);
276 mpz_clear (t);
277 mpz_clear (B1modb);
278 mpz_clear (B2modb);
279 mpz_clear (B3modb);
280 mpz_clear (B4modb);
281 mpz_clear (B5modb);
282 }
283
284 int
sumspills(mpz_t ppp,mpz_t * a,int n)285 sumspills (mpz_t ppp, mpz_t *a, int n)
286 {
287 mpz_t s;
288 int i, ret;
289
290 mpz_init_set (s, a[0]);
291
292 for (i = 1; i < n; i++)
293 {
294 mpz_add (s, s, a[i]);
295 }
296 ret = mpz_cmp (s, B) >= 0;
297 mpz_clear (s);
298
299 return ret;
300 }
301