gen-fac.c

   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 }