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