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