contrib/gmp/gen-fac_ui.c

   1 /* Generate mpz_fac_ui data.
   2
   3 Copyright 2002 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 "dumbmp.c"
  24
  25
  26 /* sets x=y*(y+2)*(y+4)*....*(y+2*(z-1))        */
  27 void
  28 odd_products (mpz_t x, mpz_t y, int z)
  29 {
  30   mpz_t t;
  31
  32   mpz_init_set (t, y);
  33   mpz_set_ui (x, 1);
  34   for (; z != 0; z--)
  35     {
  36       mpz_mul (x, x, t);
  37       mpz_add_ui (t, t, 2);
  38     }
  39   mpz_clear (t);
  40   return;
  41 }
  42
  43 /* returns 0 on success         */
  44 int
  45 gen_consts (int numb, int nail, int limb)
  46 {
  47   mpz_t x, y, z, t;
  48   unsigned long a, b, first = 1;
  49
  50   printf ("/* This file is automatically generated by gen-fac_ui.c */\n\n");
  51   printf ("#if GMP_NUMB_BITS != %d\n", numb);
  52   printf ("Error , error this data is for %d GMP_NUMB_BITS only\n", numb);
  53   printf ("#endif\n");
  54   printf ("#if GMP_LIMB_BITS != %d\n", limb);
  55   printf ("Error , error this data is for %d GMP_LIMB_BITS only\n", limb);
  56   printf ("#endif\n");
  57
  58   printf
  59     ("/* This table is 0!,1!,2!,3!,...,n! where n! has <= GMP_NUMB_BITS bits */\n");
  60   printf
  61     ("#define ONE_LIMB_FACTORIAL_TABLE CNST_LIMB(0x1),CNST_LIMB(0x1),CNST_LIMB(0x2),");
  62   mpz_init_set_ui (x, 2);
  63   for (b = 3;; b++)
  64     {
  65       mpz_mul_ui (x, x, b);     /* so b!=a       */
  66       if (mpz_sizeinbase (x, 2) > numb)
  67         break;
  68       if (first)
  69         {
  70           first = 0;
  71         }
  72       else
  73         {
  74           printf ("),");
  75         }
  76       printf ("CNST_LIMB(0x");
  77       mpz_out_str (stdout, 16, x);
  78     }
  79   printf (")\n");
  80
  81
  82   mpz_set_ui (x, 1);
  83   mpz_mul_2exp (x, x, limb + 1);        /* x=2^(limb+1)        */
  84   mpz_init (y);
  85   mpz_set_ui (y, 10000);
  86   mpz_mul (x, x, y);            /* x=2^(limb+1)*10^4     */
  87   mpz_set_ui (y, 27182);        /* exp(1)*10^4      */
  88   mpz_tdiv_q (x, x, y);         /* x=2^(limb+1)/exp(1)        */
  89   printf ("\n/* is 2^(GMP_LIMB_BITS+1)/exp(1) */\n");
  90   printf ("#define FAC2OVERE CNST_LIMB(0x");
  91   mpz_out_str (stdout, 16, x);
  92   printf (")\n");
  93
  94
  95   printf
  96     ("\n/* FACMULn is largest odd x such that x*(x+2)*...*(x+2(n-1))<=2^GMP_NUMB_BITS-1 */\n\n");
  97   mpz_init (z);
  98   mpz_init (t);
  99   for (a = 2; a <= 4; a++)
 100     {
 101       mpz_set_ui (x, 1);
 102       mpz_mul_2exp (x, x, numb);
 103       mpz_root (x, x, a);
 104       /* so x is approx sol       */
 105       if (mpz_even_p (x))
 106         mpz_sub_ui (x, x, 1);
 107       mpz_set_ui (y, 1);
 108       mpz_mul_2exp (y, y, numb);
 109       mpz_sub_ui (y, y, 1);
 110       /* decrement x until we are <= real sol     */
 111       do
 112         {
 113           mpz_sub_ui (x, x, 2);
 114           odd_products (t, x, a);
 115           if (mpz_cmp (t, y) <= 0)
 116             break;
 117         }
 118       while (1);
 119       /* increment x until > real sol     */
 120       do
 121         {
 122           mpz_add_ui (x, x, 2);
 123           odd_products (t, x, a);
 124           if (mpz_cmp (t, y) > 0)
 125             break;
 126         }
 127       while (1);
 128       /* dec once to get real sol */
 129       mpz_sub_ui (x, x, 2);
 130       printf ("#define FACMUL%lu CNST_LIMB(0x", a);
 131       mpz_out_str (stdout, 16, x);
 132       printf (")\n");
 133     }
 134
 135   return 0;
 136 }
 137
 138 int
 139 main (int argc, char *argv[])
 140 {
 141   int nail_bits, limb_bits, numb_bits;
 142
 143   if (argc != 3)
 144     {
 145       fprintf (stderr, "Usage: gen-fac_ui limbbits nailbits\n");
 146       exit (1);
 147     }
 148   limb_bits = atoi (argv[1]);
 149   nail_bits = atoi (argv[2]);
 150   numb_bits = limb_bits - nail_bits;
 151   if (limb_bits < 0 || nail_bits < 0 || numb_bits < 0)
 152     {
 153       fprintf (stderr, "Invalid limb/nail bits %d,%d\n", limb_bits,
 154                nail_bits);
 155       exit (1);
 156     }
 157   gen_consts (numb_bits, nail_bits, limb_bits);
 158   return 0;
 159 }