1 /* Generate mpz_fac_ui data.
3 Copyright 2002 Free Software Foundation, Inc.
5 This file is part of the GNU MP Library.
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.
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.
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/. */
26 /* sets x=y*(y+2)*(y+4)*....*(y+2*(z-1)) */
28 odd_products (mpz_t x, mpz_t y, int z)
43 /* returns 0 on success */
45 gen_consts (int numb, int nail, int limb)
48 unsigned long a, b, first = 1;
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);
54 printf ("#if GMP_LIMB_BITS != %d\n", limb);
55 printf ("Error , error this data is for %d GMP_LIMB_BITS only\n", limb);
59 ("/* This table is 0!,1!,2!,3!,...,n! where n! has <= GMP_NUMB_BITS bits */\n");
61 ("#define ONE_LIMB_FACTORIAL_TABLE CNST_LIMB(0x1),CNST_LIMB(0x1),CNST_LIMB(0x2),");
62 mpz_init_set_ui (x, 2);
65 mpz_mul_ui (x, x, b); /* so b!=a */
66 if (mpz_sizeinbase (x, 2) > numb)
76 printf ("CNST_LIMB(0x");
77 mpz_out_str (stdout, 16, x);
83 mpz_mul_2exp (x, x, limb + 1); /* x=2^(limb+1) */
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);
96 ("\n/* FACMULn is largest odd x such that x*(x+2)*...*(x+2(n-1))<=2^GMP_NUMB_BITS-1 */\n\n");
99 for (a = 2; a <= 4; a++)
102 mpz_mul_2exp (x, x, numb);
104 /* so x is approx sol */
106 mpz_sub_ui (x, x, 1);
108 mpz_mul_2exp (y, y, numb);
109 mpz_sub_ui (y, y, 1);
110 /* decrement x until we are <= real sol */
113 mpz_sub_ui (x, x, 2);
114 odd_products (t, x, a);
115 if (mpz_cmp (t, y) <= 0)
119 /* increment x until > real sol */
122 mpz_add_ui (x, x, 2);
123 odd_products (t, x, a);
124 if (mpz_cmp (t, y) > 0)
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);
139 main (int argc, char *argv[])
141 int nail_bits, limb_bits, numb_bits;
145 fprintf (stderr, "Usage: gen-fac_ui limbbits nailbits\n");
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)
153 fprintf (stderr, "Invalid limb/nail bits %d,%d\n", limb_bits,
157 gen_consts (numb_bits, nail_bits, limb_bits);