1/* Generate mpz_fac_ui data.
2
3Copyright 2002 Free Software Foundation, Inc.
4
5This file is part of the GNU MP Library.
6
7The GNU MP Library is free software; you can redistribute it and/or modify
8it under the terms of the GNU Lesser General Public License as published by
9the Free Software Foundation; either version 3 of the License, or (at your
10option) any later version.
11
12The GNU MP Library is distributed in the hope that it will be useful, but
13WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
15License for more details.
16
17You should have received a copy of the GNU Lesser General Public License
18along 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))	*/
27void
28odd_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		*/
44int
45gen_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
138int
139main (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