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