1 /* Big numbers for Emacs.
2
3 Copyright 2018-2021 Free Software Foundation, Inc.
4
5 This file is part of GNU Emacs.
6
7 GNU Emacs is free software: you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation, either version 3 of the License, or (at
10 your option) any later version.
11
12 GNU Emacs is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with GNU Emacs. If not, see <https://www.gnu.org/licenses/>. */
19
20 #include <config.h>
21
22 #include "bignum.h"
23
24 #include "lisp.h"
25
26 #include <math.h>
27 #include <stdlib.h>
28
29 /* mpz global temporaries. Making them global saves the trouble of
30 properly using mpz_init and mpz_clear on temporaries even when
31 storage is exhausted. Admittedly this is not ideal. An mpz value
32 in a temporary is made permanent by mpz_swapping it with a bignum's
33 value. Although typically at most two temporaries are needed,
34 rounddiv_q and rounding_driver both need four and time_arith needs
35 five. */
36
37 mpz_t mpz[5];
38
39 static void *
xrealloc_for_gmp(void * ptr,size_t ignore,size_t size)40 xrealloc_for_gmp (void *ptr, size_t ignore, size_t size)
41 {
42 return xrealloc (ptr, size);
43 }
44
45 static void
xfree_for_gmp(void * ptr,size_t ignore)46 xfree_for_gmp (void *ptr, size_t ignore)
47 {
48 xfree (ptr);
49 }
50
51 void
init_bignum(void)52 init_bignum (void)
53 {
54 eassert (mp_bits_per_limb == GMP_NUMB_BITS);
55 integer_width = 1 << 16;
56
57 /* FIXME: The Info node `(gmp) Custom Allocation' states: "No error
58 return is allowed from any of these functions, if they return
59 then they must have performed the specified operation. [...]
60 There's currently no defined way for the allocation functions to
61 recover from an error such as out of memory, they must terminate
62 program execution. A 'longjmp' or throwing a C++ exception will
63 have undefined results." But xmalloc and xrealloc do call
64 'longjmp'. */
65 mp_set_memory_functions (xmalloc, xrealloc_for_gmp, xfree_for_gmp);
66
67 for (int i = 0; i < ARRAYELTS (mpz); i++)
68 mpz_init (mpz[i]);
69 }
70
71 /* Return the value of the Lisp bignum N, as a double. */
72 double
bignum_to_double(Lisp_Object n)73 bignum_to_double (Lisp_Object n)
74 {
75 return mpz_get_d_rounded (*xbignum_val (n));
76 }
77
78 /* Return D, converted to a Lisp integer. Discard any fraction.
79 Signal an error if D cannot be converted. */
80 Lisp_Object
double_to_integer(double d)81 double_to_integer (double d)
82 {
83 if (!isfinite (d))
84 overflow_error ();
85 mpz_set_d (mpz[0], d);
86 return make_integer_mpz ();
87 }
88
89 /* Return a Lisp integer equal to mpz[0], which has BITS bits and which
90 must not be in fixnum range. Set mpz[0] to a junk value. */
91 static Lisp_Object
make_bignum_bits(size_t bits)92 make_bignum_bits (size_t bits)
93 {
94 /* The documentation says integer-width should be nonnegative, so
95 comparing it to BITS works even though BITS is unsigned. Treat
96 integer-width as if it were at least twice the machine integer width,
97 so that timefns.c can safely use bignums for double-precision
98 timestamps. */
99 if (integer_width < bits && 2 * max (INTMAX_WIDTH, UINTMAX_WIDTH) < bits)
100 overflow_error ();
101
102 struct Lisp_Bignum *b = ALLOCATE_PLAIN_PSEUDOVECTOR (struct Lisp_Bignum,
103 PVEC_BIGNUM);
104 mpz_init (b->value);
105 mpz_swap (b->value, mpz[0]);
106 return make_lisp_ptr (b, Lisp_Vectorlike);
107 }
108
109 /* Return a Lisp integer equal to mpz[0], which must not be in fixnum range.
110 Set mpz[0] to a junk value. */
111 static Lisp_Object
make_bignum(void)112 make_bignum (void)
113 {
114 return make_bignum_bits (mpz_sizeinbase (mpz[0], 2));
115 }
116
117 /* Return a Lisp integer equal to N, which must not be in fixnum range. */
118 Lisp_Object
make_bigint(intmax_t n)119 make_bigint (intmax_t n)
120 {
121 eassert (FIXNUM_OVERFLOW_P (n));
122 mpz_set_intmax (mpz[0], n);
123 return make_bignum ();
124 }
125 Lisp_Object
make_biguint(uintmax_t n)126 make_biguint (uintmax_t n)
127 {
128 eassert (FIXNUM_OVERFLOW_P (n));
129 mpz_set_uintmax (mpz[0], n);
130 return make_bignum ();
131 }
132
133 /* Return a Lisp integer equal to -N, which must not be in fixnum range. */
134 Lisp_Object
make_neg_biguint(uintmax_t n)135 make_neg_biguint (uintmax_t n)
136 {
137 eassert (-MOST_NEGATIVE_FIXNUM < n);
138 mpz_set_uintmax (mpz[0], n);
139 mpz_neg (mpz[0], mpz[0]);
140 return make_bignum ();
141 }
142
143 /* Return a Lisp integer with value taken from mpz[0].
144 Set mpz[0] to a junk value. */
145 Lisp_Object
make_integer_mpz(void)146 make_integer_mpz (void)
147 {
148 size_t bits = mpz_sizeinbase (mpz[0], 2);
149
150 if (bits <= FIXNUM_BITS)
151 {
152 EMACS_INT v = 0;
153 int i = 0, shift = 0;
154
155 do
156 {
157 EMACS_INT limb = mpz_getlimbn (mpz[0], i++);
158 v += limb << shift;
159 shift += GMP_NUMB_BITS;
160 }
161 while (shift < bits);
162
163 if (mpz_sgn (mpz[0]) < 0)
164 v = -v;
165
166 if (!FIXNUM_OVERFLOW_P (v))
167 return make_fixnum (v);
168 }
169
170 return make_bignum_bits (bits);
171 }
172
173 /* Set RESULT to V. This code is for when intmax_t is wider than long. */
174 void
mpz_set_intmax_slow(mpz_t result,intmax_t v)175 mpz_set_intmax_slow (mpz_t result, intmax_t v)
176 {
177 int maxlimbs = (INTMAX_WIDTH + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS;
178 mp_limb_t *limb = mpz_limbs_write (result, maxlimbs);
179 int n = 0;
180 uintmax_t u = v;
181 bool negative = v < 0;
182 if (negative)
183 {
184 uintmax_t two = 2;
185 u = -u & ((two << (UINTMAX_WIDTH - 1)) - 1);
186 }
187
188 do
189 {
190 limb[n++] = u;
191 u = GMP_NUMB_BITS < UINTMAX_WIDTH ? u >> GMP_NUMB_BITS : 0;
192 }
193 while (u != 0);
194
195 mpz_limbs_finish (result, negative ? -n : n);
196 }
197 void
mpz_set_uintmax_slow(mpz_t result,uintmax_t v)198 mpz_set_uintmax_slow (mpz_t result, uintmax_t v)
199 {
200 int maxlimbs = (UINTMAX_WIDTH + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS;
201 mp_limb_t *limb = mpz_limbs_write (result, maxlimbs);
202 int n = 0;
203
204 do
205 {
206 limb[n++] = v;
207 v = GMP_NUMB_BITS < INTMAX_WIDTH ? v >> GMP_NUMB_BITS : 0;
208 }
209 while (v != 0);
210
211 mpz_limbs_finish (result, n);
212 }
213
214 /* If Z fits into *PI, store its value there and return true.
215 Return false otherwise. */
216 bool
mpz_to_intmax(mpz_t const z,intmax_t * pi)217 mpz_to_intmax (mpz_t const z, intmax_t *pi)
218 {
219 ptrdiff_t bits = mpz_sizeinbase (z, 2);
220 bool negative = mpz_sgn (z) < 0;
221
222 if (bits < INTMAX_WIDTH)
223 {
224 intmax_t v = 0;
225 int i = 0, shift = 0;
226
227 do
228 {
229 intmax_t limb = mpz_getlimbn (z, i++);
230 v += limb << shift;
231 shift += GMP_NUMB_BITS;
232 }
233 while (shift < bits);
234
235 *pi = negative ? -v : v;
236 return true;
237 }
238 if (bits == INTMAX_WIDTH && INTMAX_MIN < -INTMAX_MAX && negative
239 && mpz_scan1 (z, 0) == INTMAX_WIDTH - 1)
240 {
241 *pi = INTMAX_MIN;
242 return true;
243 }
244 return false;
245 }
246 bool
mpz_to_uintmax(mpz_t const z,uintmax_t * pi)247 mpz_to_uintmax (mpz_t const z, uintmax_t *pi)
248 {
249 if (mpz_sgn (z) < 0)
250 return false;
251 ptrdiff_t bits = mpz_sizeinbase (z, 2);
252 if (UINTMAX_WIDTH < bits)
253 return false;
254
255 uintmax_t v = 0;
256 int i = 0, shift = 0;
257
258 do
259 {
260 uintmax_t limb = mpz_getlimbn (z, i++);
261 v += limb << shift;
262 shift += GMP_NUMB_BITS;
263 }
264 while (shift < bits);
265
266 *pi = v;
267 return true;
268 }
269
270 /* Return the value of the bignum X if it fits, 0 otherwise.
271 A bignum cannot be zero, so 0 indicates failure reliably. */
272 intmax_t
bignum_to_intmax(Lisp_Object x)273 bignum_to_intmax (Lisp_Object x)
274 {
275 intmax_t i;
276 return mpz_to_intmax (*xbignum_val (x), &i) ? i : 0;
277 }
278 uintmax_t
bignum_to_uintmax(Lisp_Object x)279 bignum_to_uintmax (Lisp_Object x)
280 {
281 uintmax_t i;
282 return mpz_to_uintmax (*xbignum_val (x), &i) ? i : 0;
283 }
284
285
286 /* Multiply and exponentiate mpz_t values without aborting due to size
287 limits. */
288
289 /* GMP tests for this value and aborts (!) if it is exceeded.
290 This is as of GMP 6.1.2 (2016); perhaps future versions will differ. */
291 enum { GMP_NLIMBS_MAX = min (INT_MAX, ULONG_MAX / GMP_NUMB_BITS) };
292
293 /* An upper bound on limb counts, needed to prevent libgmp and/or
294 Emacs from aborting or otherwise misbehaving. This bound applies
295 to estimates of mpz_t sizes before the mpz_t objects are created,
296 as opposed to integer-width which operates on mpz_t values after
297 creation and before conversion to Lisp bignums. */
298 enum
299 {
300 NLIMBS_LIMIT = min (min (/* libgmp needs to store limb counts. */
301 GMP_NLIMBS_MAX,
302
303 /* Size calculations need to work. */
304 min (PTRDIFF_MAX, SIZE_MAX) / sizeof (mp_limb_t)),
305
306 /* Emacs puts bit counts into fixnums. */
307 MOST_POSITIVE_FIXNUM / GMP_NUMB_BITS)
308 };
309
310 /* Like mpz_size, but tell the compiler the result is a nonnegative int. */
311
312 static int
emacs_mpz_size(mpz_t const op)313 emacs_mpz_size (mpz_t const op)
314 {
315 mp_size_t size = mpz_size (op);
316 eassume (0 <= size && size <= INT_MAX);
317 return size;
318 }
319
320 /* Wrappers to work around GMP limitations. As of GMP 6.1.2 (2016),
321 the library code aborts when a number is too large. These wrappers
322 avoid the problem for functions that can return numbers much larger
323 than their arguments. For slowly-growing numbers, the integer
324 width checks in bignum.c should suffice. */
325
326 void
emacs_mpz_mul(mpz_t rop,mpz_t const op1,mpz_t const op2)327 emacs_mpz_mul (mpz_t rop, mpz_t const op1, mpz_t const op2)
328 {
329 if (NLIMBS_LIMIT - emacs_mpz_size (op1) < emacs_mpz_size (op2))
330 overflow_error ();
331 mpz_mul (rop, op1, op2);
332 }
333
334 void
emacs_mpz_mul_2exp(mpz_t rop,mpz_t const op1,EMACS_INT op2)335 emacs_mpz_mul_2exp (mpz_t rop, mpz_t const op1, EMACS_INT op2)
336 {
337 /* Fudge factor derived from GMP 6.1.2, to avoid an abort in
338 mpz_mul_2exp (look for the '+ 1' in its source code). */
339 enum { mul_2exp_extra_limbs = 1 };
340 enum { lim = min (NLIMBS_LIMIT, GMP_NLIMBS_MAX - mul_2exp_extra_limbs) };
341
342 EMACS_INT op2limbs = op2 / GMP_NUMB_BITS;
343 if (lim - emacs_mpz_size (op1) < op2limbs)
344 overflow_error ();
345 mpz_mul_2exp (rop, op1, op2);
346 }
347
348 void
emacs_mpz_pow_ui(mpz_t rop,mpz_t const base,unsigned long exp)349 emacs_mpz_pow_ui (mpz_t rop, mpz_t const base, unsigned long exp)
350 {
351 /* This fudge factor is derived from GMP 6.1.2, to avoid an abort in
352 mpz_n_pow_ui (look for the '5' in its source code). */
353 enum { pow_ui_extra_limbs = 5 };
354 enum { lim = min (NLIMBS_LIMIT, GMP_NLIMBS_MAX - pow_ui_extra_limbs) };
355
356 int nbase = emacs_mpz_size (base), n;
357 if (INT_MULTIPLY_WRAPV (nbase, exp, &n) || lim < n)
358 overflow_error ();
359 mpz_pow_ui (rop, base, exp);
360 }
361
362
363 /* Yield an upper bound on the buffer size needed to contain a C
364 string representing the NUM in base BASE. This includes any
365 preceding '-' and the terminating null. */
366 static ptrdiff_t
mpz_bufsize(mpz_t const num,int base)367 mpz_bufsize (mpz_t const num, int base)
368 {
369 return mpz_sizeinbase (num, base) + 2;
370 }
371 ptrdiff_t
bignum_bufsize(Lisp_Object num,int base)372 bignum_bufsize (Lisp_Object num, int base)
373 {
374 return mpz_bufsize (*xbignum_val (num), base);
375 }
376
377 /* Convert NUM to a nearest double, as opposed to mpz_get_d which
378 truncates toward zero. */
379 double
mpz_get_d_rounded(mpz_t const num)380 mpz_get_d_rounded (mpz_t const num)
381 {
382 ptrdiff_t size = mpz_bufsize (num, 10);
383
384 /* Use mpz_get_d as a shortcut for a bignum so small that rounding
385 errors cannot occur, which is possible if EMACS_INT (not counting
386 sign) has fewer bits than a double significand. */
387 if (! ((FLT_RADIX == 2 && DBL_MANT_DIG <= FIXNUM_BITS - 1)
388 || (FLT_RADIX == 16 && DBL_MANT_DIG * 4 <= FIXNUM_BITS - 1))
389 && size <= DBL_DIG + 2)
390 return mpz_get_d (num);
391
392 USE_SAFE_ALLOCA;
393 char *buf = SAFE_ALLOCA (size);
394 mpz_get_str (buf, 10, num);
395 double result = strtod (buf, NULL);
396 SAFE_FREE ();
397 return result;
398 }
399
400 /* Store into BUF (of size SIZE) the value of NUM as a base-BASE string.
401 If BASE is negative, use upper-case digits in base -BASE.
402 Return the string's length.
403 SIZE must equal bignum_bufsize (NUM, abs (BASE)). */
404 ptrdiff_t
bignum_to_c_string(char * buf,ptrdiff_t size,Lisp_Object num,int base)405 bignum_to_c_string (char *buf, ptrdiff_t size, Lisp_Object num, int base)
406 {
407 eassert (bignum_bufsize (num, abs (base)) == size);
408 mpz_get_str (buf, base, *xbignum_val (num));
409 ptrdiff_t n = size - 2;
410 return !buf[n - 1] ? n - 1 : n + !!buf[n];
411 }
412
413 /* Convert NUM to a base-BASE Lisp string.
414 If BASE is negative, use upper-case digits in base -BASE. */
415
416 Lisp_Object
bignum_to_string(Lisp_Object num,int base)417 bignum_to_string (Lisp_Object num, int base)
418 {
419 ptrdiff_t size = bignum_bufsize (num, abs (base));
420 USE_SAFE_ALLOCA;
421 char *str = SAFE_ALLOCA (size);
422 ptrdiff_t len = bignum_to_c_string (str, size, num, base);
423 Lisp_Object result = make_unibyte_string (str, len);
424 SAFE_FREE ();
425 return result;
426 }
427
428 /* Create a bignum by scanning NUM, with digits in BASE.
429 NUM must consist of an optional '-', a nonempty sequence
430 of base-BASE digits, and a terminating null byte, and
431 the represented number must not be in fixnum range. */
432
433 Lisp_Object
make_bignum_str(char const * num,int base)434 make_bignum_str (char const *num, int base)
435 {
436 struct Lisp_Bignum *b = ALLOCATE_PLAIN_PSEUDOVECTOR (struct Lisp_Bignum,
437 PVEC_BIGNUM);
438 mpz_init (b->value);
439 int check = mpz_set_str (b->value, num, base);
440 eassert (check == 0);
441 return make_lisp_ptr (b, Lisp_Vectorlike);
442 }
443
444 /* Check that X is a Lisp integer in the range LO..HI.
445 Return X's value as an intmax_t. */
446
447 intmax_t
check_integer_range(Lisp_Object x,intmax_t lo,intmax_t hi)448 check_integer_range (Lisp_Object x, intmax_t lo, intmax_t hi)
449 {
450 CHECK_INTEGER (x);
451 intmax_t i;
452 if (! (integer_to_intmax (x, &i) && lo <= i && i <= hi))
453 args_out_of_range_3 (x, make_int (lo), make_int (hi));
454 return i;
455 }
456
457 /* Check that X is a Lisp integer in the range 0..HI.
458 Return X's value as an uintmax_t. */
459
460 uintmax_t
check_uinteger_max(Lisp_Object x,uintmax_t hi)461 check_uinteger_max (Lisp_Object x, uintmax_t hi)
462 {
463 CHECK_INTEGER (x);
464 uintmax_t i;
465 if (! (integer_to_uintmax (x, &i) && i <= hi))
466 args_out_of_range_3 (x, make_fixnum (0), make_uint (hi));
467 return i;
468 }
469
470 /* Check that X is a Lisp integer no greater than INT_MAX,
471 and return its value or zero, whichever is greater. */
472
473 int
check_int_nonnegative(Lisp_Object x)474 check_int_nonnegative (Lisp_Object x)
475 {
476 CHECK_INTEGER (x);
477 return NILP (Fnatnump (x)) ? 0 : check_integer_range (x, 0, INT_MAX);
478 }
479