1 /* gen-trialdivtab.c 2 3 Contributed to the GNU project by Torbjorn Granlund. 4 5 Copyright 2009 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 "dumbmp.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 mpz_t y; 207 unsigned long cnt; 208 209 mpz_init (y); 210 mpz_set (y, x); 211 cnt = 0; 212 while (mpz_sgn (y) != 0) 213 { 214 mpz_tdiv_q_2exp (y, y, 1); 215 cnt++; 216 } 217 mpz_clear (y); 218 219 return cnt; 220 } 221 222 void 223 mpn_mod_1s_4p_cps (mpz_t cps[7], mpz_t bparm) 224 { 225 mpz_t b, bi; 226 mpz_t B1modb, B2modb, B3modb, B4modb, B5modb; 227 mpz_t t; 228 int cnt; 229 230 mpz_init_set (b, bparm); 231 232 cnt = limb_bits - mpz_log2 (b); 233 234 mpz_init (bi); 235 mpz_init (t); 236 mpz_init (B1modb); 237 mpz_init (B2modb); 238 mpz_init (B3modb); 239 mpz_init (B4modb); 240 mpz_init (B5modb); 241 242 mpz_set_ui (t, 1); 243 mpz_mul_2exp (t, t, limb_bits - cnt); 244 mpz_sub (t, t, b); 245 mpz_mul_2exp (t, t, limb_bits); 246 mpz_tdiv_q (bi, t, b); /* bi = B^2/b, except msb */ 247 248 mpz_set_ui (t, 1); 249 mpz_mul_2exp (t, t, limb_bits); /* t = B */ 250 mpz_tdiv_r (B1modb, t, b); 251 252 mpz_mul_2exp (t, B1modb, limb_bits); 253 mpz_tdiv_r (B2modb, t, b); 254 255 mpz_mul_2exp (t, B2modb, limb_bits); 256 mpz_tdiv_r (B3modb, t, b); 257 258 mpz_mul_2exp (t, B3modb, limb_bits); 259 mpz_tdiv_r (B4modb, t, b); 260 261 mpz_mul_2exp (t, B4modb, limb_bits); 262 mpz_tdiv_r (B5modb, t, b); 263 264 mpz_set (cps[0], bi); 265 mpz_set_ui (cps[1], cnt); 266 mpz_tdiv_q_2exp (cps[2], B1modb, 0); 267 mpz_tdiv_q_2exp (cps[3], B2modb, 0); 268 mpz_tdiv_q_2exp (cps[4], B3modb, 0); 269 mpz_tdiv_q_2exp (cps[5], B4modb, 0); 270 mpz_tdiv_q_2exp (cps[6], B5modb, 0); 271 272 mpz_clear (b); 273 mpz_clear (bi); 274 mpz_clear (t); 275 mpz_clear (B1modb); 276 mpz_clear (B2modb); 277 mpz_clear (B3modb); 278 mpz_clear (B4modb); 279 mpz_clear (B5modb); 280 } 281 282 int 283 sumspills (mpz_t ppp, mpz_t *a, int n) 284 { 285 mpz_t s; 286 int i, ret; 287 288 mpz_init_set (s, a[0]); 289 290 for (i = 1; i < n; i++) 291 { 292 mpz_add (s, s, a[i]); 293 } 294 ret = mpz_cmp (s, B) >= 0; 295 mpz_clear (s); 296 297 return ret; 298 } 299