root/src/bignum.c

DEFINITIONS

1. xrealloc_for_gmp
2. xfree_for_gmp
3. init_bignum
4. bignum_to_double
5. double_to_integer
6. make_bignum_bits
7. make_bignum
8. make_bigint
9. make_biguint
10. make_neg_biguint
11. make_integer_mpz
12. mpz_set_intmax_slow
13. mpz_set_uintmax_slow
14. mpz_to_intmax
15. mpz_to_uintmax
16. bignum_to_intmax
17. bignum_to_uintmax
18. emacs_mpz_size
19. emacs_mpz_mul
20. emacs_mpz_mul_2exp
21. emacs_mpz_pow_ui
22. mpz_bufsize
23. bignum_bufsize
24. mpz_get_d_rounded
25. bignum_to_c_string
26. bignum_to_string
27. make_bignum_str
28. check_integer_range
29. check_uinteger_max
30. check_int_nonnegative
31. get_random_limb
32. get_random_limb_lim
33. get_random_bignum

     1 /* Big numbers for Emacs.
     2 
     3 Copyright 2018-2023 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 *
    40 xrealloc_for_gmp (void *ptr, size_t ignore, size_t size)
    41 {
    42   return xrealloc (ptr, size);
    43 }
    44 
    45 static void
    46 xfree_for_gmp (void *ptr, size_t ignore)
    47 {
    48   xfree (ptr);
    49 }
    50 
    51 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
    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
    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
    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
   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
   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
   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
   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
   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
   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
   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
   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
   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
   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
   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
   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
   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
   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
   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 (ckd_mul (&n, nbase, exp) || 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
   367 mpz_bufsize (mpz_t const num, int base)
   368 {
   369   return mpz_sizeinbase (num, base) + 2;
   370 }
   371 ptrdiff_t
   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
   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
   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
   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
   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
   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
   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
   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 
   480 /* Return a random mp_limb_t.  */
   481 
   482 static mp_limb_t
   483 get_random_limb (void)
   484 {
   485   if (GMP_NUMB_BITS <= ULONG_WIDTH)
   486     return get_random_ulong ();
   487 
   488   /* Work around GCC -Wshift-count-overflow false alarm.  */
   489   int shift = GMP_NUMB_BITS <= ULONG_WIDTH ? 0 : ULONG_WIDTH;
   490 
   491   /* This is in case someone builds GMP with unusual definitions for
   492      MINI_GMP_LIMB_TYPE or _LONG_LONG_LIMB.  */
   493   mp_limb_t r = 0;
   494   for (int i = 0; i < GMP_NUMB_BITS; i += ULONG_WIDTH)
   495     r = (r << shift) | get_random_ulong ();
   496   return r;
   497 }
   498 
   499 /* Return a random mp_limb_t I in the range 0 <= I < LIM.
   500    If LIM is zero, simply return a random mp_limb_t.  */
   501 
   502 static mp_limb_t
   503 get_random_limb_lim (mp_limb_t lim)
   504 {
   505   /* Return the remainder of a random mp_limb_t R divided by LIM,
   506      except reject the rare case where R is so close to the maximum
   507      mp_limb_t that the remainder isn't random.  */
   508   mp_limb_t difflim = - lim, diff, remainder;
   509   do
   510     {
   511       mp_limb_t r = get_random_limb ();
   512       if (lim == 0)
   513         return r;
   514       remainder = r % lim;
   515       diff = r - remainder;
   516     }
   517   while (difflim < diff);
   518 
   519   return remainder;
   520 }
   521 
   522 /* Return a random Lisp integer I in the range 0 <= I < LIMIT,
   523    where LIMIT is a positive bignum.  */
   524 
   525 Lisp_Object
   526 get_random_bignum (struct Lisp_Bignum const *limit)
   527 {
   528   mpz_t const *lim = bignum_val (limit);
   529   mp_size_t nlimbs = mpz_size (*lim);
   530   eassume (0 < nlimbs);
   531   mp_limb_t *r_limb = mpz_limbs_write (mpz[0], nlimbs);
   532   mp_limb_t const *lim_limb = mpz_limbs_read (*lim);
   533   mp_limb_t limhi = lim_limb[nlimbs - 1];
   534   eassert (limhi);
   535   bool edgy;
   536 
   537   do
   538     {
   539       /* Generate the result one limb at a time, most significant first.
   540          Choose the most significant limb RHI randomly from 0..LIMHI,
   541          where LIMHI is the LIM's first limb, except choose from
   542          0..(LIMHI-1) if there is just one limb.  RHI == LIMHI is an
   543          unlucky edge case as later limbs might cause the result to be
   544          exceed or equal LIM; if this happens, it causes another
   545          iteration in the outer loop.  */
   546 
   547       mp_limb_t rhi = get_random_limb_lim (limhi + (1 < nlimbs));
   548       edgy = rhi == limhi;
   549       r_limb[nlimbs - 1] = rhi;
   550 
   551       for (mp_size_t i = nlimbs - 1; 0 < i--; )
   552         {
   553           /* get_random_limb_lim (edgy ? limb_lim[i] + 1 : 0)
   554              would be wrong here, as the full mp_limb_t range is
   555              needed in later limbs for the edge case to have the
   556              proper weighting.  */
   557           mp_limb_t ri = get_random_limb ();
   558           if (edgy)
   559             {
   560               if (lim_limb[i] < ri)
   561                 break;
   562               edgy = lim_limb[i] == ri;
   563             }
   564           r_limb[i] = ri;
   565         }
   566     }
   567   while (edgy);
   568 
   569   mpz_limbs_finish (mpz[0], nlimbs);
   570   return make_integer_mpz ();
   571 }