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 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 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 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