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