1 /* mpfr-mini-gmp.c -- Interface functions for mini-gmp. 2 3 Copyright 2014-2020 Free Software Foundation, Inc. 4 Contributed by the AriC and Caramba projects, INRIA. 5 6 This file is part of the GNU MPFR Library. 7 8 The GNU MPFR Library is free software; you can redistribute it and/or modify 9 it under the terms of the GNU Lesser General Public License as published by 10 the Free Software Foundation; either version 3 of the License, or (at your 11 option) any later version. 12 13 The GNU MPFR Library is distributed in the hope that it will be useful, but 14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16 License for more details. 17 18 You should have received a copy of the GNU Lesser General Public License 19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23 /* The following include will do 2 things: include the config.h 24 if there is one (as it may define MPFR_USE_MINI_GMP), and avoid 25 an empty translation unit (see ISO C99, 6.9). */ 26 #include "mpfr-impl.h" 27 28 #ifdef MPFR_USE_MINI_GMP 29 30 /************************ random generation functions ************************/ 31 32 #ifdef WANT_gmp_randinit_default 33 void 34 gmp_randinit_default (gmp_randstate_t state) 35 { 36 } 37 #endif 38 39 #ifdef WANT_gmp_randseed_ui 40 void 41 gmp_randseed_ui (gmp_randstate_t state, unsigned long int seed) 42 { 43 /* With a portable version of the conversion from unsigned long to long 44 (at least GCC and Clang optimize this expression to identity). */ 45 srand48 (seed > LONG_MAX ? -1 - (long) ~seed : (long) seed); 46 } 47 #endif 48 49 #ifdef WANT_gmp_randclear 50 void 51 gmp_randclear (gmp_randstate_t state) 52 { 53 } 54 #endif 55 56 #ifdef WANT_gmp_randinit_set 57 void 58 gmp_randinit_set (gmp_randstate_t s1, gmp_randstate_t s2) 59 { 60 } 61 #endif 62 63 static mp_limb_t 64 random_limb (void) 65 { 66 /* lrand48() returns a random number in [0, 2^31-1], 67 but the low 15 bits do not depend on the random seed, 68 thus it is safer to use the upper bits */ 69 #if GMP_NUMB_BITS < 32 70 /* use the upper GMP_NUMB_BITS bits from lrand48 () */ 71 return (mp_limb_t) (lrand48 () >> (31 - GMP_NUMB_BITS)); 72 #elif GMP_NUMB_BITS == 32 73 /* use the upper 16 bits from two lrand48 calls */ 74 return (lrand48 () >> 15) + ((lrand48 () >> 15) << 16); 75 #elif GMP_NUMB_BITS == 64 76 /* use the upper 16 bits from four lrand48 calls */ 77 return (lrand48 () >> 15) + ((((mp_limb_t) lrand48 ()) >> 15) << 16) 78 + ((((mp_limb_t) lrand48 ()) >> 15) << 32) 79 + ((((mp_limb_t) lrand48 ()) >> 15) << 48); 80 #else 81 #error "GMP_NUMB_BITS should be 8, 16, 32 or >= 64" 82 #endif 83 } 84 85 #ifdef WANT_mpz_urandomb 86 void 87 mpz_urandomb (mpz_t rop, gmp_randstate_t state, mp_bitcnt_t nbits) 88 { 89 unsigned long n, i; 90 91 mpz_realloc2 (rop, nbits); 92 n = (nbits - 1) / GMP_NUMB_BITS + 1; /* number of limbs */ 93 for (i = n; i-- > 0;) 94 rop->_mp_d[i] = random_limb (); 95 i = n * GMP_NUMB_BITS - nbits; 96 /* mask the upper i bits */ 97 if (i) 98 rop->_mp_d[n-1] = MPFR_LIMB_LSHIFT(rop->_mp_d[n-1], i) >> i; 99 while (n > 0 && (rop->_mp_d[n-1] == 0)) 100 n--; 101 rop->_mp_size = n; 102 } 103 #endif 104 105 #ifdef WANT_gmp_urandomm_ui 106 /* generates a random unsigned long */ 107 static unsigned long 108 random_ulong (void) 109 { 110 #ifdef MPFR_LONG_WITHIN_LIMB 111 /* we assume a limb and an unsigned long have both a number of different 112 values that is a power of two, thus when we cast a random limb into 113 an unsigned long, we still get an uniform distribution */ 114 return random_limb (); 115 #else 116 /* with the same assumption as above, we need to generate as many random 117 limbs needed to "fill" an unsigned long */ 118 unsigned long u, v; 119 120 v = MPFR_LIMB_MAX; 121 u = random_limb (); 122 while (v < ULONG_MAX) 123 { 124 v = (v << GMP_NUMB_BITS) + MPFR_LIMB_MAX; 125 u = (u << GMP_NUMB_BITS) + random_limb (); 126 } 127 return u; 128 #endif 129 } 130 131 unsigned long 132 gmp_urandomm_ui (gmp_randstate_t state, unsigned long n) 133 { 134 unsigned long p, q, r; 135 136 MPFR_ASSERTN (n > 0); 137 p = random_ulong (); /* p is in [0, ULONG_MAX], thus p is uniform among 138 ULONG_MAX+1 values */ 139 q = n * (ULONG_MAX / n); 140 r = ULONG_MAX % n; 141 if (r != n - 1) /* ULONG_MAX+1 is not multiple of n, will happen whenever 142 n is not a power of two */ 143 while (p >= q) 144 p = random_ulong (); 145 return p % n; 146 } 147 #endif 148 149 #ifdef WANT_gmp_urandomb_ui 150 unsigned long 151 gmp_urandomb_ui (gmp_randstate_t state, unsigned long n) 152 { 153 #ifdef MPFR_LONG_WITHIN_LIMB 154 return random_limb () % (1UL << n); 155 #else 156 unsigned long res = 0; 157 int m = n; /* remaining bits to generate */ 158 while (m >= GMP_NUMB_BITS) 159 { 160 /* we can generate a full limb */ 161 res = (res << GMP_NUMB_BITS) | (unsigned long) random_limb (); 162 m -= GMP_NUMB_BITS; 163 } 164 /* now m < GMP_NUMB_BITS */ 165 if (m) /* generate m extra bits */ 166 res = (res << m) | (unsigned long) (random_limb () % (1UL << m)); 167 return res; 168 #endif 169 } 170 #endif 171 172 /************************* division functions ********************************/ 173 174 #ifdef WANT_mpn_divrem_1 175 mp_limb_t 176 mpn_divrem_1 (mp_limb_t *qp, mp_size_t qxn, mp_limb_t *np, mp_size_t nn, 177 mp_limb_t d0) 178 { 179 mpz_t q, r, n, d; 180 mp_limb_t ret, dd[1]; 181 182 d->_mp_d = dd; 183 d->_mp_d[0] = d0; 184 d->_mp_size = 1; 185 mpz_init (q); 186 mpz_init (r); 187 if (qxn == 0) 188 { 189 n->_mp_d = np; 190 n->_mp_size = nn; 191 } 192 else 193 { 194 mpz_init2 (n, (nn + qxn) * GMP_NUMB_BITS); 195 mpn_copyi (n->_mp_d + qxn, np, nn); 196 mpn_zero (n->_mp_d, qxn); 197 n->_mp_size = nn + qxn; 198 } 199 mpz_tdiv_qr (q, r, n, d); 200 if (q->_mp_size > 0) 201 mpn_copyi (qp, q->_mp_d, q->_mp_size); 202 if (q->_mp_size < nn + qxn) 203 mpn_zero (qp + q->_mp_size, nn + qxn - q->_mp_size); 204 ret = (r->_mp_size == 1) ? r->_mp_d[0] : 0; 205 mpz_clear (q); 206 mpz_clear (r); 207 if (qxn != 0) 208 mpz_clear (n); 209 return ret; 210 } 211 #endif 212 213 #ifdef WANT_mpn_divrem 214 mp_limb_t 215 mpn_divrem (mp_limb_t *qp, mp_size_t qn, mp_limb_t *np, 216 mp_size_t nn, const mp_limb_t *dp, mp_size_t dn) 217 { 218 mpz_t q, r, n, d; 219 mp_limb_t ret; 220 221 MPFR_ASSERTN(qn == 0); 222 qn = nn - dn; 223 n->_mp_d = np; 224 n->_mp_size = nn; 225 d->_mp_d = (mp_limb_t*) dp; 226 d->_mp_size = dn; 227 mpz_init (q); 228 mpz_init (r); 229 mpz_tdiv_qr (q, r, n, d); 230 MPFR_ASSERTN(q->_mp_size == qn || q->_mp_size == qn + 1); 231 mpn_copyi (qp, q->_mp_d, qn); 232 ret = (q->_mp_size == qn) ? 0 : q->_mp_d[qn]; 233 if (r->_mp_size > 0) 234 mpn_copyi (np, r->_mp_d, r->_mp_size); 235 if (r->_mp_size < dn) 236 mpn_zero (np + r->_mp_size, dn - r->_mp_size); 237 mpz_clear (q); 238 mpz_clear (r); 239 return ret; 240 } 241 #endif 242 243 #ifdef WANT_mpn_tdiv_qr 244 void 245 mpn_tdiv_qr (mp_limb_t *qp, mp_limb_t *rp, mp_size_t qxn, 246 const mp_limb_t *np, mp_size_t nn, 247 const mp_limb_t *dp, mp_size_t dn) 248 { 249 mpz_t q, r, n, d; 250 251 MPFR_ASSERTN(qxn == 0); 252 n->_mp_d = (mp_limb_t*) np; 253 n->_mp_size = nn; 254 d->_mp_d = (mp_limb_t*) dp; 255 d->_mp_size = dn; 256 mpz_init (q); 257 mpz_init (r); 258 mpz_tdiv_qr (q, r, n, d); 259 MPFR_ASSERTN(q->_mp_size > 0); 260 mpn_copyi (qp, q->_mp_d, q->_mp_size); 261 if (q->_mp_size < nn - dn + 1) 262 mpn_zero (qp + q->_mp_size, nn - dn + 1 - q->_mp_size); 263 if (r->_mp_size > 0) 264 mpn_copyi (rp, r->_mp_d, r->_mp_size); 265 if (r->_mp_size < dn) 266 mpn_zero (rp + r->_mp_size, dn - r->_mp_size); 267 mpz_clear (q); 268 mpz_clear (r); 269 } 270 #endif 271 272 #if 0 /* this function is useful for debugging, thus please keep it here */ 273 void 274 mpz_dump (mpz_t z) 275 { 276 mp_size_t n = z->_mp_size; 277 278 MPFR_STAT_STATIC_ASSERT ((GMP_NUMB_BITS % 4) == 0); 279 280 if (n == 0) 281 printf ("0"); 282 else 283 { 284 int first = 1; 285 if (n < 0) 286 { 287 printf ("-"); 288 n = -n; 289 } 290 while (n > 0) 291 { 292 if (first) 293 { 294 printf ("%lx", (unsigned long) z->_mp_d[n-1]); 295 first = 0; 296 } 297 else 298 { 299 char s[17]; 300 int len; 301 sprintf (s, "%lx", (unsigned long) z->_mp_d[n-1]); 302 len = strlen (s); 303 /* one character should be printed for 4 bits */ 304 while (len++ < GMP_NUMB_BITS / 4) 305 printf ("0"); 306 printf ("%lx", (unsigned long) z->_mp_d[n-1]); 307 } 308 n--; 309 } 310 } 311 printf ("\n"); 312 } 313 #endif 314 315 #endif /* MPFR_USE_MINI_GMP */ 316