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
odd_products(mpz_t x,mpz_t y,int z)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
gen_consts(int numb,int nail,int limb)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
main(int argc,char * argv[])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 }
160