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