1 /* Generate data for combinatorics: fac_ui, bin_uiui, ...
3 Copyright 2002, 2011-2013 Free Software Foundation, Inc.
5 This file is part of the GNU MP Library.
7 The GNU MP Library is free software; you can redistribute it and/or modify
8 it under the terms of either:
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.
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
20 or both in parallel, as here.
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
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/. */
34 #include "bootstrap.c"
37 mpz_remove_twos (mpz_t x
)
39 mp_bitcnt_t r
= mpz_scan1(x
, 0);
40 mpz_tdiv_q_2exp (x
, x
, r
);
44 /* returns 0 on success */
46 gen_consts (int numb
, int nail
, int limb
)
48 mpz_t x
, mask
, y
, last
;
50 unsigned long ofl
, ofe
;
52 printf ("/* This file is automatically generated by gen-fac.c */\n\n");
53 printf ("#if GMP_NUMB_BITS != %d\n", numb
);
54 printf ("Error , error this data is for %d GMP_NUMB_BITS only\n", numb
);
57 printf ("#if GMP_LIMB_BITS != %d\n", limb
);
58 printf ("Error , error this data is for %d GMP_LIMB_BITS only\n", limb
);
63 ("/* This table is 0!,1!,2!,3!,...,n! where n! has <= GMP_NUMB_BITS bits */\n");
65 ("#define ONE_LIMB_FACTORIAL_TABLE CNST_LIMB(0x1),CNST_LIMB(0x1");
66 mpz_init_set_ui (x
, 1);
70 mpz_mul_ui (x
, x
, b
); /* so b!=a */
71 if (mpz_sizeinbase (x
, 2) > numb
)
73 printf ("),CNST_LIMB(0x");
74 mpz_out_str (stdout
, 16, x
);
79 ("\n/* This table is 0!,1!,2!/2,3!/2,...,n!/2^sn where n!/2^sn is an */\n");
81 ("/* odd integer for each n, and n!/2^sn has <= GMP_NUMB_BITS bits */\n");
83 ("#define ONE_LIMB_ODD_FACTORIAL_TABLE CNST_LIMB(0x1),CNST_LIMB(0x1),CNST_LIMB(0x1");
87 for (a
= b
; (a
& 1) == 0; a
>>= 1);
89 mpz_mul_ui (x
, last
, a
);
90 if (mpz_sizeinbase (x
, 2) > numb
)
92 printf ("),CNST_LIMB(0x");
93 mpz_out_str (stdout
, 16, x
);
97 ("#define ODD_FACTORIAL_TABLE_MAX CNST_LIMB(0x");
98 mpz_out_str (stdout
, 16, last
);
103 ("#define ODD_FACTORIAL_TABLE_LIMIT (%lu)\n", ofl
);
104 mpz_init2 (mask
, numb
);
105 mpz_setbit (mask
, numb
);
106 mpz_sub_ui (mask
, mask
, 1);
108 ("\n/* Previous table, continued, values modulo 2^GMP_NUMB_BITS */\n");
110 ("#define ONE_LIMB_ODD_FACTORIAL_EXTTABLE CNST_LIMB(0x");
111 mpz_and (x
, x
, mask
);
112 mpz_out_str (stdout
, 16, x
);
114 mpz_bin_uiui (y
, b
, b
/2);
118 for (a
= b
; (a
& 1) == 0; a
>>= 1);
120 mpz_divexact_ui (y
, y
, a
/2+1);
121 mpz_mul_ui (y
, y
, a
);
123 mpz_mul_2exp (y
, y
, 1);
124 if (mpz_sizeinbase (y
, 2) > numb
)
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
);
134 ("#define ODD_FACTORIAL_EXTTABLE_LIMIT (%lu)\n", ofe
);
137 ("\n/* This table is 1!!,3!!,...,(2n+1)!! where (2n+1)!! has <= GMP_NUMB_BITS bits */\n");
139 ("#define ONE_LIMB_ODD_DOUBLEFACTORIAL_TABLE CNST_LIMB(0x1");
144 mpz_mul_ui (x
, last
, b
);
145 if (mpz_sizeinbase (x
, 2) > numb
)
147 printf ("),CNST_LIMB(0x");
148 mpz_out_str (stdout
, 16, x
);
152 ("#define ODD_DOUBLEFACTORIAL_TABLE_MAX CNST_LIMB(0x");
153 mpz_out_str (stdout
, 16, last
);
157 ("#define ODD_DOUBLEFACTORIAL_TABLE_LIMIT (%lu)\n", b
- 2);
160 ("\n/* This table x_1, x_2,... contains values s.t. x_n^n has <= GMP_NUMB_BITS bits */\n");
162 ("#define NTH_ROOT_NUMB_MASK_TABLE (GMP_NUMB_MASK");
163 for (b
= 2;b
<= 8; b
++)
165 mpz_root (x
, mask
, b
);
166 printf ("),CNST_LIMB(0x");
167 mpz_out_str (stdout
, 16, x
);
171 mpz_add_ui (mask
, mask
, 1);
173 ("\n/* This table contains inverses of odd factorials, modulo 2^GMP_NUMB_BITS */\n");
175 ("\n/* It begins with (2!/2)^-1=1 */\n");
177 ("#define ONE_LIMB_ODD_FACTORIAL_INVERSES_TABLE CNST_LIMB(0x1");
179 for (b
= 3;b
<= ofe
- 2; b
++)
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
);
189 ofe
= (ofe
/ 16 + 1) * 16;
192 ("\n/* This table contains 2n-popc(2n) for small n */\n");
194 ("\n/* It begins with 2-1=1 (n=1) */\n");
196 ("#define TABLE_2N_MINUS_POPC_2N 1");
197 for (b
= 4; b
<= ofe
; b
+= 2)
200 printf (",%lu",b
- mpz_popcount (x
));
204 ("#define TABLE_LIMIT_2N_MINUS_POPC_2N %lu\n", ofe
+ 1);
209 ("#define ODD_CENTRAL_BINOMIAL_OFFSET (%lu)\n", ofl
);
211 ("\n/* This table contains binomial(2k,k)/2^t */\n");
213 ("\n/* It begins with ODD_CENTRAL_BINOMIAL_TABLE_MIN */\n");
215 ("#define ONE_LIMB_ODD_CENTRAL_BINOMIAL_TABLE ");
218 mpz_bin_uiui (x
, 2 * b
, b
);
220 if (mpz_sizeinbase (x
, 2) > numb
)
224 printf("CNST_LIMB(0x");
225 mpz_out_str (stdout
, 16, x
);
231 ("#define ODD_CENTRAL_BINOMIAL_TABLE_LIMIT (%lu)\n", ofe
);
234 ("\n/* This table contains the inverses of elements in the previous table. */\n");
236 ("#define ONE_LIMB_ODD_CENTRAL_BINOMIAL_INVERSE_TABLE CNST_LIMB(0x");
237 for (b
= ofl
; b
<= ofe
; b
++)
239 mpz_bin_uiui (x
, 2 * b
, b
);
241 mpz_invert (x
, x
, mask
);
242 mpz_out_str (stdout
, 16, x
);
244 printf ("),CNST_LIMB(0x");
249 ("\n/* This table contains the values t in the formula binomial(2k,k)/2^t */\n");
251 ("#define CENTRAL_BINOMIAL_2FAC_TABLE ");
252 for (b
= ofl
; b
<= ofe
; b
++)
254 mpz_bin_uiui (x
, 2 * b
, b
);
255 printf ("%d", mpz_remove_twos (x
));
265 main (int argc
, char *argv
[])
267 int nail_bits
, limb_bits
, numb_bits
;
271 fprintf (stderr
, "Usage: gen-fac_ui limbbits nailbits\n");
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)
279 fprintf (stderr
, "Invalid limb/nail bits %d,%d\n", limb_bits
,
283 gen_consts (numb_bits
, nail_bits
, limb_bits
);