1/* gen-trialdivtab.c
2
3   Contributed to the GNU project by Torbjorn Granlund.
4
5Copyright 2009 Free Software Foundation, Inc.
6
7This file is part of the GNU MP Library.
8
9The GNU MP Library is free software; you can redistribute it and/or modify
10it under the terms of the GNU Lesser General Public License as published by
11the Free Software Foundation; either version 3 of the License, or (at your
12option) any later version.
13
14The GNU MP Library is distributed in the hope that it will be useful, but
15WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
17License for more details.
18
19You should have received a copy of the GNU Lesser General Public License
20along 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
41int sumspills (mpz_t, mpz_t *, int);
42void mpn_mod_1s_4p_cps (mpz_t [7], mpz_t);
43
44int limb_bits;
45
46mpz_t B;
47
48int
49main (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
203unsigned long
204mpz_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
222void
223mpn_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
282int
283sumspills (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