source/libs/gmp/gmp-src/gen-fac.c

   1 /* Generate data for combinatorics: fac_ui, bin_uiui, ...
   2
   3 Copyright 2002, 2011-2015 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
  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
  46 gen_consts (int numb, int nail, int 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 != %d\n", numb);
  54   printf ("Error , error this data is for %d GMP_NUMB_BITS only\n", numb);
  55   printf ("#endif\n");
  56 #if 0
  57   printf ("#if GMP_LIMB_BITS != %d\n", limb);
  58   printf ("Error , error this data is for %d 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_init2 (mask, numb + 1);
 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
 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, nail_bits, limb_bits);
 284   return 0;
 285 }