1 /* gen-trialdivtab.c 2 3 Contributed to the GNU project by Torbjorn Granlund. 4 5 Copyright 2009, 2012, 2013, 2016, 2018 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 60 main (int argc, char *argv[]) 61 { 62 int 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 (B); 95 mpz_setbit (B, limb_bits); 96 mpz_init_set (gmp_numb_max, B); 97 mpz_sub_ui (gmp_numb_max, gmp_numb_max, 1); 98 99 mpz_init (tmp); 100 mpz_init (inv); 101 102 mpz_init (Bhalf); 103 mpz_setbit (Bhalf, limb_bits - 1); 104 105 start_p = 3; 106 107 mpz_init_set_ui (ppp, 1); 108 mpz_init (acc); 109 interval_start = start_p; 110 omitted_p = 3; 111 interval_end = 0; 112 113 /* printf ("static struct gmp_primes_dtab gmp_primes_dtab[] = {\n"); */ 114 115 printf ("#ifdef WANT_dtab\n"); 116 117 for (t = start_p; t <= end_p; t += 2) 118 { 119 if (! isprime (t)) 120 continue; 121 122 mpz_mul_ui (acc, ppp, t); 123 stop = mpz_cmp (acc, Bhalf) >= 0; 124 if (!stop) 125 { 126 mpn_mod_1s_4p_cps (pre, acc); 127 stop = sumspills (acc, pre + 2, 5); 128 } 129 130 if (stop) 131 { 132 for (p = interval_start; p <= interval_end; p += 2) 133 { 134 if (! isprime (p)) 135 continue; 136 137 printf (" P(%d,", (int) p); 138 mpz_invert_ui_2exp (inv, p, limb_bits); 139 printf ("CNST_LIMB(0x"); mpz_out_str (stdout, 16, inv); printf ("),"); 140 141 mpz_tdiv_q_ui (tmp, gmp_numb_max, p); 142 printf ("CNST_LIMB(0x"); mpz_out_str (stdout, 16, tmp); 143 printf (")),\n"); 144 } 145 mpz_set_ui (ppp, t); 146 interval_start = t; 147 omitted_p = t; 148 } 149 else 150 { 151 mpz_set (ppp, acc); 152 } 153 interval_end = t; 154 } 155 printf ("#define SMALLEST_OMITTED_PRIME %d\n", (int) omitted_p); 156 printf ("#endif\n"); 157 158 printf ("#ifdef WANT_ptab\n"); 159 160 /* printf ("static struct gmp_primes_ptab gmp_primes_ptab[] = {\n"); */ 161 162 endtok = ""; 163 164 mpz_set_ui (ppp, 1); 165 interval_start = start_p; 166 interval_end = 0; 167 np = 0; 168 start_idx = 0; 169 for (t = start_p; t <= end_p; t += 2) 170 { 171 if (! isprime (t)) 172 continue; 173 174 mpz_mul_ui (acc, ppp, t); 175 176 stop = mpz_cmp (acc, Bhalf) >= 0; 177 if (!stop) 178 { 179 mpn_mod_1s_4p_cps (pre, acc); 180 stop = sumspills (acc, pre + 2, 5); 181 } 182 183 if (stop) 184 { 185 mpn_mod_1s_4p_cps (pre, ppp); 186 printf ("%s", endtok); 187 printf (" {CNST_LIMB(0x"); mpz_out_str (stdout, 16, ppp); 188 printf ("),{CNST_LIMB(0x"); mpz_out_str (stdout, 16, pre[0]); 189 printf ("),%d", (int) PTR(pre[1])[0]); 190 for (i = 0; i < 5; i++) 191 { 192 printf (","); 193 printf ("CNST_LIMB(0x"); mpz_out_str (stdout, 16, pre[2 + i]); 194 printf (")"); 195 } 196 printf ("},"); 197 printf ("%d,", start_idx); 198 printf ("%d}", np - start_idx); 199 200 endtok = ",\n"; 201 mpz_set_ui (ppp, t); 202 interval_start = t; 203 start_idx = np; 204 } 205 else 206 { 207 mpz_set (ppp, acc); 208 } 209 interval_end = t; 210 np++; 211 } 212 213 printf ("\n"); 214 printf ("#endif\n"); 215 216 return 0; 217 } 218 219 unsigned long 220 mpz_log2 (mpz_t x) 221 { 222 return mpz_sgn (x) ? mpz_sizeinbase (x, 2) : 0; 223 } 224 225 void 226 mpn_mod_1s_4p_cps (mpz_t cps[7], mpz_t bparm) 227 { 228 mpz_t b, bi; 229 mpz_t B1modb, B2modb, B3modb, B4modb, B5modb; 230 mpz_t t; 231 int cnt; 232 233 mpz_init_set (b, bparm); 234 235 cnt = limb_bits - mpz_log2 (b); 236 237 mpz_init (bi); 238 mpz_init (t); 239 mpz_init (B1modb); 240 mpz_init (B2modb); 241 mpz_init (B3modb); 242 mpz_init (B4modb); 243 mpz_init (B5modb); 244 245 mpz_set_ui (t, 1); 246 mpz_mul_2exp (t, t, limb_bits - cnt); 247 mpz_sub (t, t, b); 248 mpz_mul_2exp (t, t, limb_bits); 249 mpz_tdiv_q (bi, t, b); /* bi = B^2/b, except msb */ 250 251 mpz_set_ui (t, 1); 252 mpz_mul_2exp (t, t, limb_bits); /* t = B */ 253 mpz_tdiv_r (B1modb, t, b); 254 255 mpz_mul_2exp (t, B1modb, limb_bits); 256 mpz_tdiv_r (B2modb, t, b); 257 258 mpz_mul_2exp (t, B2modb, limb_bits); 259 mpz_tdiv_r (B3modb, t, b); 260 261 mpz_mul_2exp (t, B3modb, limb_bits); 262 mpz_tdiv_r (B4modb, t, b); 263 264 mpz_mul_2exp (t, B4modb, limb_bits); 265 mpz_tdiv_r (B5modb, t, b); 266 267 mpz_set (cps[0], bi); 268 mpz_set_ui (cps[1], cnt); 269 mpz_tdiv_q_2exp (cps[2], B1modb, 0); 270 mpz_tdiv_q_2exp (cps[3], B2modb, 0); 271 mpz_tdiv_q_2exp (cps[4], B3modb, 0); 272 mpz_tdiv_q_2exp (cps[5], B4modb, 0); 273 mpz_tdiv_q_2exp (cps[6], B5modb, 0); 274 275 mpz_clear (b); 276 mpz_clear (bi); 277 mpz_clear (t); 278 mpz_clear (B1modb); 279 mpz_clear (B2modb); 280 mpz_clear (B3modb); 281 mpz_clear (B4modb); 282 mpz_clear (B5modb); 283 } 284 285 int 286 sumspills (mpz_t ppp, mpz_t *a, int n) 287 { 288 mpz_t s; 289 int i, ret; 290 291 mpz_init_set (s, a[0]); 292 293 for (i = 1; i < n; i++) 294 { 295 mpz_add (s, s, a[i]); 296 } 297 ret = mpz_cmp (s, B) >= 0; 298 mpz_clear (s); 299 300 return ret; 301 } 302