Home | History | Annotate | Line # | Download | only in dist 
      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
     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 (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
    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