1 /* Generate data for combinatorics: fac_ui, bin_uiui, ...
2
3 Copyright 2002, 2011-2016 Free Software Foundation, Inc.
4
5 This file is part of the GNU MP Library.
6
7 The GNU MP Library is free software; you can redistribute it and/or modify
8 it under the terms of either:
9
10 * the GNU Lesser General Public License as published by the Free
11 Software Foundation; either version 3 of the License, or (at your
12 option) any later version.
13
14 or
15
16 * the GNU General Public License as published by the Free Software
17 Foundation; either version 2 of the License, or (at your option) any
18 later version.
19
20 or both in parallel, as here.
21
22 The GNU MP Library is distributed in the hope that it will be useful, but
23 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
24 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
25 for more details.
26
27 You should have received copies of the GNU General Public License and the
28 GNU Lesser General Public License along with the GNU MP Library. If not,
29 see https://www.gnu.org/licenses/. */
30
31 #include <stdio.h>
32 #include <stdlib.h>
33
34 #include "bootstrap.c"
35
36 int
mpz_remove_twos(mpz_t x)37 mpz_remove_twos (mpz_t x)
38 {
39 mp_bitcnt_t r = mpz_scan1(x, 0);
40 mpz_tdiv_q_2exp (x, x, r);
41 return r;
42 }
43
44 /* returns 0 on success */
45 int
gen_consts(unsigned numb,unsigned limb)46 gen_consts (unsigned numb, unsigned limb)
47 {
48 mpz_t x, mask, y, last;
49 unsigned long a, b;
50 unsigned long ofl, ofe;
51
52 printf ("/* This file is automatically generated by gen-fac.c */\n\n");
53 printf ("#if GMP_NUMB_BITS != %u\n", numb);
54 printf ("Error , error this data is for %u GMP_NUMB_BITS only\n", numb);
55 printf ("#endif\n");
56 #if 0
57 printf ("#if GMP_LIMB_BITS != %u\n", limb);
58 printf ("Error , error this data is for %u GMP_LIMB_BITS only\n", limb);
59 printf ("#endif\n");
60 #endif
61
62 printf
63 ("/* This table is 0!,1!,2!,3!,...,n! where n! has <= GMP_NUMB_BITS bits */\n");
64 printf
65 ("#define ONE_LIMB_FACTORIAL_TABLE CNST_LIMB(0x1),CNST_LIMB(0x1");
66 mpz_init_set_ui (x, 1);
67 mpz_init (last);
68 for (b = 2;; b++)
69 {
70 mpz_mul_ui (x, x, b); /* so b!=a */
71 if (mpz_sizeinbase (x, 2) > numb)
72 break;
73 printf ("),CNST_LIMB(0x");
74 mpz_out_str (stdout, 16, x);
75 }
76 printf (")\n");
77
78 printf
79 ("\n/* This table is 0!,1!,2!/2,3!/2,...,n!/2^sn where n!/2^sn is an */\n");
80 printf
81 ("/* odd integer for each n, and n!/2^sn has <= GMP_NUMB_BITS bits */\n");
82 printf
83 ("#define ONE_LIMB_ODD_FACTORIAL_TABLE CNST_LIMB(0x1),CNST_LIMB(0x1),CNST_LIMB(0x1");
84 mpz_set_ui (x, 1);
85 for (b = 3;; b++)
86 {
87 for (a = b; (a & 1) == 0; a >>= 1);
88 mpz_swap (last, x);
89 mpz_mul_ui (x, last, a);
90 if (mpz_sizeinbase (x, 2) > numb)
91 break;
92 printf ("),CNST_LIMB(0x");
93 mpz_out_str (stdout, 16, x);
94 }
95 printf (")\n");
96 printf
97 ("#define ODD_FACTORIAL_TABLE_MAX CNST_LIMB(0x");
98 mpz_out_str (stdout, 16, last);
99 printf (")\n");
100
101 ofl = b - 1;
102 printf
103 ("#define ODD_FACTORIAL_TABLE_LIMIT (%lu)\n", ofl);
104 mpz_init (mask);
105 mpz_setbit (mask, numb);
106 mpz_sub_ui (mask, mask, 1);
107 printf
108 ("\n/* Previous table, continued, values modulo 2^GMP_NUMB_BITS */\n");
109 printf
110 ("#define ONE_LIMB_ODD_FACTORIAL_EXTTABLE CNST_LIMB(0x");
111 mpz_and (x, x, mask);
112 mpz_out_str (stdout, 16, x);
113 mpz_init (y);
114 mpz_bin_uiui (y, b, b/2);
115 b++;
116 for (;; b++)
117 {
118 for (a = b; (a & 1) == 0; a >>= 1);
119 if (a == b) {
120 mpz_divexact_ui (y, y, a/2+1);
121 mpz_mul_ui (y, y, a);
122 } else
123 mpz_mul_2exp (y, y, 1);
124 if (mpz_sizeinbase (y, 2) > numb)
125 break;
126 mpz_mul_ui (x, x, a);
127 mpz_and (x, x, mask);
128 printf ("),CNST_LIMB(0x");
129 mpz_out_str (stdout, 16, x);
130 }
131 printf (")\n");
132 ofe = b - 1;
133 printf
134 ("#define ODD_FACTORIAL_EXTTABLE_LIMIT (%lu)\n", ofe);
135
136 printf
137 ("\n/* This table is 1!!,3!!,...,(2n+1)!! where (2n+1)!! has <= GMP_NUMB_BITS bits */\n");
138 printf
139 ("#define ONE_LIMB_ODD_DOUBLEFACTORIAL_TABLE CNST_LIMB(0x1");
140 mpz_set_ui (x, 1);
141 for (b = 3;; b+=2)
142 {
143 mpz_swap (last, x);
144 mpz_mul_ui (x, last, b);
145 if (mpz_sizeinbase (x, 2) > numb)
146 break;
147 printf ("),CNST_LIMB(0x");
148 mpz_out_str (stdout, 16, x);
149 }
150 printf (")\n");
151 printf
152 ("#define ODD_DOUBLEFACTORIAL_TABLE_MAX CNST_LIMB(0x");
153 mpz_out_str (stdout, 16, last);
154 printf (")\n");
155
156 printf
157 ("#define ODD_DOUBLEFACTORIAL_TABLE_LIMIT (%lu)\n", b - 2);
158
159 printf
160 ("\n/* This table x_1, x_2,... contains values s.t. x_n^n has <= GMP_NUMB_BITS bits */\n");
161 printf
162 ("#define NTH_ROOT_NUMB_MASK_TABLE (GMP_NUMB_MASK");
163 for (b = 2;b <= 8; b++)
164 {
165 mpz_root (x, mask, b);
166 printf ("),CNST_LIMB(0x");
167 mpz_out_str (stdout, 16, x);
168 }
169 printf (")\n");
170
171 mpz_add_ui (mask, mask, 1);
172 printf
173 ("\n/* This table contains inverses of odd factorials, modulo 2^GMP_NUMB_BITS */\n");
174 printf
175 ("\n/* It begins with (2!/2)^-1=1 */\n");
176 printf
177 ("#define ONE_LIMB_ODD_FACTORIAL_INVERSES_TABLE CNST_LIMB(0x1");
178 mpz_set_ui (x, 1);
179 for (b = 3;b <= ofe - 2; b++)
180 {
181 for (a = b; (a & 1) == 0; a >>= 1);
182 mpz_mul_ui (x, x, a);
183 mpz_invert (y, x, mask);
184 printf ("),CNST_LIMB(0x");
185 mpz_out_str (stdout, 16, y);
186 }
187 printf (")\n");
188
189 ofe = (ofe / 16 + 1) * 16;
190
191 printf
192 ("\n/* This table contains 2n-popc(2n) for small n */\n");
193 printf
194 ("\n/* It begins with 2-1=1 (n=1) */\n");
195 printf
196 ("#define TABLE_2N_MINUS_POPC_2N 1");
197 for (b = 4; b <= ofe; b += 2)
198 {
199 mpz_set_ui (x, b);
200 printf (",%lu",b - mpz_popcount (x));
201 }
202 printf ("\n");
203 printf
204 ("#define TABLE_LIMIT_2N_MINUS_POPC_2N %lu\n", ofe + 1);
205
206
207 ofl = (ofl + 1) / 2;
208 printf
209 ("#define ODD_CENTRAL_BINOMIAL_OFFSET (%lu)\n", ofl);
210 printf
211 ("\n/* This table contains binomial(2k,k)/2^t */\n");
212 printf
213 ("\n/* It begins with ODD_CENTRAL_BINOMIAL_TABLE_MIN */\n");
214 printf
215 ("#define ONE_LIMB_ODD_CENTRAL_BINOMIAL_TABLE ");
216 for (b = ofl;; b++)
217 {
218 mpz_bin_uiui (x, 2 * b, b);
219 mpz_remove_twos (x);
220 if (mpz_sizeinbase (x, 2) > numb)
221 break;
222 if (b != ofl)
223 printf ("),");
224 printf("CNST_LIMB(0x");
225 mpz_out_str (stdout, 16, x);
226 }
227 printf (")\n");
228
229 ofe = b - 1;
230 printf
231 ("#define ODD_CENTRAL_BINOMIAL_TABLE_LIMIT (%lu)\n", ofe);
232
233 printf
234 ("\n/* This table contains the inverses of elements in the previous table. */\n");
235 printf
236 ("#define ONE_LIMB_ODD_CENTRAL_BINOMIAL_INVERSE_TABLE CNST_LIMB(0x");
237 for (b = ofl; b <= ofe; b++)
238 {
239 mpz_bin_uiui (x, 2 * b, b);
240 mpz_remove_twos (x);
241 mpz_invert (x, x, mask);
242 mpz_out_str (stdout, 16, x);
243 if (b != ofe)
244 printf ("),CNST_LIMB(0x");
245 }
246 printf (")\n");
247
248 printf
249 ("\n/* This table contains the values t in the formula binomial(2k,k)/2^t */\n");
250 printf
251 ("#define CENTRAL_BINOMIAL_2FAC_TABLE ");
252 for (b = ofl; b <= ofe; b++)
253 {
254 mpz_bin_uiui (x, 2 * b, b);
255 printf ("%d", mpz_remove_twos (x));
256 if (b != ofe)
257 printf (",");
258 }
259 printf ("\n");
260
261 return 0;
262 }
263
264 int
main(int argc,char * argv[])265 main (int argc, char *argv[])
266 {
267 int nail_bits, limb_bits, numb_bits;
268
269 if (argc != 3)
270 {
271 fprintf (stderr, "Usage: gen-fac limbbits nailbits\n");
272 exit (1);
273 }
274 limb_bits = atoi (argv[1]);
275 nail_bits = atoi (argv[2]);
276 numb_bits = limb_bits - nail_bits;
277 if (limb_bits < 2 || nail_bits < 0 || numb_bits < 1)
278 {
279 fprintf (stderr, "Invalid limb/nail bits %d,%d\n", limb_bits,
280 nail_bits);
281 exit (1);
282 }
283 gen_consts (numb_bits, limb_bits);
284 return 0;
285 }
286