1 /* Generate primesieve data.
2
3 Contributed to the GNU project by Marco Bodrato.
4
5 Copyright 2021, 2022 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 #include <stdio.h>
34 #include "bootstrap.c"
35
36 static int
37 bit_to_n (int bit) { return (bit*3+4)|1; }
38
39 int
40 generate (int limb_bits, int limit)
41 {
42 mpz_t limb;
43 int i, lb, pc, c, totpc, maxprime;
44
45 mpz_init (limb);
46
47 printf ("/* This file generated by gen-sieve.c - DO NOT EDIT. */\n");
48 printf ("\n");
49 printf ("#if GMP_LIMB_BITS != %d\n", limb_bits);
50 printf ("Error, error, this data is for %d bits\n", limb_bits);
51 printf ("#endif\n");
52 printf ("\n");
53 printf ("#define PRIMESIEVE_INIT_TABLE ");
54
55 maxprime = 3;
56 lb = pc = c = totpc = 0;
57 for (i = 0; i < limit; i++)
58 {
59 if (! isprime (bit_to_n (i)))
60 mpz_setbit (limb, lb);
61 else
62 maxprime = bit_to_n (i), ++pc;
63 ++lb;
64 if (lb == limb_bits)
65 {
66 ++c;
67 printf ("\\\n\tCNST_LIMB (0x");
68 mpz_out_str (stdout, -16, limb);
69 printf ("),\t/* %d - %d (%d primes) */\t", bit_to_n (i + 1 - limb_bits),
70 bit_to_n (i + 1) - 1, pc);
71 totpc += pc;
72 lb = pc = 0;
73 mpz_set_ui (limb, 0);
74 }
75 }
76
77 if ((mpz_sgn (limb) | lb | pc) != 0)
78 {
79 printf ("\ngen-sieve: Internal error, during generate (%d, %d).\n", limb_bits, limit);
80 abort();
81 }
82
83 mpz_clear (limb);
84
85 printf ("\n");
86 printf ("#define PRIMESIEVE_NUMBEROF_TABLE %d\n", c);
87
88 printf ("/* #define PRIMESIEVE_PRIMES_IN_TABLE %d */\n", totpc);
89 printf ("#define PRIMESIEVE_HIGHEST_PRIME %d\n", maxprime);
90 printf ("/* #define PRIMESIEVE_FIRST_UNCHECKED %d */\n", bit_to_n (limit));
91
92 return c;
93 }
94
95 void
96 setmask (mpz_t mask, int a, int b)
97 {
98 mpz_set_ui (mask, 0);
99 for (unsigned i = 0; i < 2 * a * b; ++i)
100 if ((bit_to_n (i) % a == 0) || (bit_to_n (i) % b == 0))
101 mpz_setbit (mask, i);
102 }
103
104 void
105 gen_sieve_masks (int limb_bits) {
106 mpz_t mask, limb;
107
108 mpz_init (mask);
109 mpz_init (limb);
110
111 printf ("\n");
112 if (limb_bits > 60 && limb_bits < 91)
113 {
114 setmask (mask, 5, 11);
115
116 mpz_tdiv_r_2exp (limb, mask, limb_bits);
117 printf ("#define SIEVE_MASK1 CNST_LIMB(0x");
118 mpz_out_str (stdout, -16, limb);
119 printf (")\n");
120 mpz_tdiv_q_2exp (limb, mask, limb_bits);
121 printf ("#define SIEVE_MASKT CNST_LIMB(0x");
122 mpz_out_str (stdout, -16, limb);
123 printf (")\n");
124
125 setmask (mask, 7, 13);
126
127 mpz_tdiv_r_2exp (limb, mask, limb_bits);
128 printf ("#define SIEVE_2MSK1 CNST_LIMB(0x");
129 mpz_out_str (stdout, -16, limb);
130 printf (")\n");
131 mpz_tdiv_q_2exp (mask, mask, limb_bits);
132 mpz_tdiv_r_2exp (limb, mask, limb_bits);
133 printf ("#define SIEVE_2MSK2 CNST_LIMB(0x");
134 mpz_out_str (stdout, -16, limb);
135 printf (")\n");
136 mpz_tdiv_q_2exp (limb, mask, limb_bits);
137 printf ("#define SIEVE_2MSKT CNST_LIMB(0x");
138 mpz_out_str (stdout, -16, limb);
139 printf (")\n");
140 }
141 else if (limb_bits > 23 && limb_bits < 36)
142 {
143 setmask (mask, 5, 7);
144
145 mpz_tdiv_r_2exp (limb, mask, limb_bits);
146 printf ("#define SIEVE_MASK1 CNST_LIMB(0x");
147 mpz_out_str (stdout, -16, limb);
148 printf (")\n");
149 mpz_tdiv_q_2exp (mask, mask, limb_bits);
150 mpz_tdiv_r_2exp (limb, mask, limb_bits);
151 printf ("#define SIEVE_MASK2 CNST_LIMB(0x");
152 mpz_out_str (stdout, -16, limb);
153 printf (")\n");
154 mpz_tdiv_q_2exp (limb, mask, limb_bits);
155 printf ("#define SIEVE_MASKT CNST_LIMB(0x");
156 mpz_out_str (stdout, -16, limb);
157 printf (")\n");
158 }
159 printf ("\n");
160
161 mpz_clear (limb);
162 mpz_clear (mask);
163 }
164
165 /* 5*2 = 10
166 7*2 = 14
167 5*7*2 = 70 (2*35, 3*24, 4*18, 5*14...)
168 5*11*2 = 110 (2*55, 3*37, 4*28, 5*22...)
169 5*13*2 = 130 (2*65, 3*44, 4*33, 5*26...)
170 7*11*2 = 154 (2*77, 3*52, 4*39, 5*31...)
171 7*13*2 = 182 (2*91, 3*61, 4*46, 5*37...)
172 */
173
174 int
175 main (int argc, char *argv[])
176 {
177 int limb_bits, limit;
178
179 if (argc != 2)
180 {
181 fprintf (stderr, "Usage: gen-sieve <limbbits>\n");
182 exit (1);
183 }
184
185 limb_bits = atoi (argv[1]);
186
187 limit = 64 * 28; /* bits in the presieved sieve */
188 if (limit % limb_bits != 0)
189 limit += limb_bits - limit % limb_bits;
190 generate (limb_bits, limit);
191 gen_sieve_masks (limb_bits);
192
193 return 0;
194 }