src/bignum.c

   1 /* Big numbers for Emacs.
   2
   3 Copyright 2018 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    rounding_driver and rounddiv_q need four altogther.  */
  35
  36 mpz_t mpz[4];
  37
  38 static void *
  39 xrealloc_for_gmp (void *ptr, size_t ignore, size_t size)
  40 {
  41   return xrealloc (ptr, size);
  42 }
  43
  44 static void
  45 xfree_for_gmp (void *ptr, size_t ignore)
  46 {
  47   xfree (ptr);
  48 }
  49
  50 void
  51 init_bignum (void)
  52 {
  53   eassert (mp_bits_per_limb == GMP_NUMB_BITS);
  54   integer_width = 1 << 16;
  55   mp_set_memory_functions (xmalloc, xrealloc_for_gmp, xfree_for_gmp);
  56
  57   for (int i = 0; i < ARRAYELTS (mpz); i++)
  58     mpz_init (mpz[i]);
  59 }
  60
  61 /* Return the value of the Lisp bignum N, as a double.  */
  62 double
  63 bignum_to_double (Lisp_Object n)
  64 {
  65   return mpz_get_d_rounded (XBIGNUM (n)->value);
  66 }
  67
  68 /* Return D, converted to a Lisp integer.  Discard any fraction.
  69    Signal an error if D cannot be converted.  */
  70 Lisp_Object
  71 double_to_integer (double d)
  72 {
  73   if (!isfinite (d))
  74     overflow_error ();
  75   mpz_set_d (mpz[0], d);
  76   return make_integer_mpz ();
  77 }
  78
  79 /* Return a Lisp integer equal to mpz[0], which has BITS bits and which
  80    must not be in fixnum range.  Set mpz[0] to a junk value.  */
  81 static Lisp_Object
  82 make_bignum_bits (size_t bits)
  83 {
  84   /* The documentation says integer-width should be nonnegative, so
  85      a single comparison suffices even though 'bits' is unsigned.  */
  86   if (integer_width < bits)
  87     overflow_error ();
  88
  89   struct Lisp_Bignum *b = ALLOCATE_PSEUDOVECTOR (struct Lisp_Bignum, value,
  90                                                  PVEC_BIGNUM);
  91   mpz_init (b->value);
  92   mpz_swap (b->value, mpz[0]);
  93   return make_lisp_ptr (b, Lisp_Vectorlike);
  94 }
  95
  96 /* Return a Lisp integer equal to mpz[0], which must not be in fixnum range.
  97    Set mpz[0] to a junk value.  */
  98 static Lisp_Object
  99 make_bignum (void)
 100 {
 101   return make_bignum_bits (mpz_sizeinbase (mpz[0], 2));
 102 }
 103
 104 static void mpz_set_uintmax_slow (mpz_t, uintmax_t);
 105
 106 /* Set RESULT to V.  */
 107 static void
 108 mpz_set_uintmax (mpz_t result, uintmax_t v)
 109 {
 110   if (v <= ULONG_MAX)
 111     mpz_set_ui (result, v);
 112   else
 113     mpz_set_uintmax_slow (result, v);
 114 }
 115
 116 /* Return a Lisp integer equal to N, which must not be in fixnum range.  */
 117 Lisp_Object
 118 make_bigint (intmax_t n)
 119 {
 120   eassert (FIXNUM_OVERFLOW_P (n));
 121   mpz_set_intmax (mpz[0], n);
 122   return make_bignum ();
 123 }
 124 Lisp_Object
 125 make_biguint (uintmax_t n)
 126 {
 127   eassert (FIXNUM_OVERFLOW_P (n));
 128   mpz_set_uintmax (mpz[0], n);
 129   return make_bignum ();
 130 }
 131
 132 /* Return a Lisp integer with value taken from mpz[0].
 133    Set mpz[0] to a junk value.  */
 134 Lisp_Object
 135 make_integer_mpz (void)
 136 {
 137   size_t bits = mpz_sizeinbase (mpz[0], 2);
 138
 139   if (bits <= FIXNUM_BITS)
 140     {
 141       EMACS_INT v = 0;
 142       int i = 0, shift = 0;
 143
 144       do
 145         {
 146           EMACS_INT limb = mpz_getlimbn (mpz[0], i++);
 147           v += limb << shift;
 148           shift += GMP_NUMB_BITS;
 149         }
 150       while (shift < bits);
 151
 152       if (mpz_sgn (mpz[0]) < 0)
 153         v = -v;
 154
 155       if (!FIXNUM_OVERFLOW_P (v))
 156         return make_fixnum (v);
 157     }
 158
 159   return make_bignum_bits (bits);
 160 }
 161
 162 /* Set RESULT to V.  This code is for when intmax_t is wider than long.  */
 163 void
 164 mpz_set_intmax_slow (mpz_t result, intmax_t v)
 165 {
 166   int maxlimbs = (INTMAX_WIDTH + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS;
 167   mp_limb_t *limb = mpz_limbs_write (result, maxlimbs);
 168   int n = 0;
 169   uintmax_t u = v;
 170   bool negative = v < 0;
 171   if (negative)
 172     {
 173       uintmax_t two = 2;
 174       u = -u & ((two << (UINTMAX_WIDTH - 1)) - 1);
 175     }
 176
 177   do
 178     {
 179       limb[n++] = u;
 180       u = GMP_NUMB_BITS < UINTMAX_WIDTH ? u >> GMP_NUMB_BITS : 0;
 181     }
 182   while (u != 0);
 183
 184   mpz_limbs_finish (result, negative ? -n : n);
 185 }
 186 static void
 187 mpz_set_uintmax_slow (mpz_t result, uintmax_t v)
 188 {
 189   int maxlimbs = (UINTMAX_WIDTH + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS;
 190   mp_limb_t *limb = mpz_limbs_write (result, maxlimbs);
 191   int n = 0;
 192
 193   do
 194     {
 195       limb[n++] = v;
 196       v = GMP_NUMB_BITS < INTMAX_WIDTH ? v >> GMP_NUMB_BITS : 0;
 197     }
 198   while (v != 0);
 199
 200   mpz_limbs_finish (result, n);
 201 }
 202
 203 /* Return the value of the bignum X if it fits, 0 otherwise.
 204    A bignum cannot be zero, so 0 indicates failure reliably.  */
 205 intmax_t
 206 bignum_to_intmax (Lisp_Object x)
 207 {
 208   ptrdiff_t bits = mpz_sizeinbase (XBIGNUM (x)->value, 2);
 209   bool negative = mpz_sgn (XBIGNUM (x)->value) < 0;
 210
 211   if (bits < INTMAX_WIDTH)
 212     {
 213       intmax_t v = 0;
 214       int i = 0, shift = 0;
 215
 216       do
 217         {
 218           intmax_t limb = mpz_getlimbn (XBIGNUM (x)->value, i++);
 219           v += limb << shift;
 220           shift += GMP_NUMB_BITS;
 221         }
 222       while (shift < bits);
 223
 224       return negative ? -v : v;
 225     }
 226   return ((bits == INTMAX_WIDTH && INTMAX_MIN < -INTMAX_MAX && negative
 227            && mpz_scan1 (XBIGNUM (x)->value, 0) == INTMAX_WIDTH - 1)
 228           ? INTMAX_MIN : 0);
 229 }
 230 uintmax_t
 231 bignum_to_uintmax (Lisp_Object x)
 232 {
 233   uintmax_t v = 0;
 234   if (0 <= mpz_sgn (XBIGNUM (x)->value))
 235     {
 236       ptrdiff_t bits = mpz_sizeinbase (XBIGNUM (x)->value, 2);
 237       if (bits <= UINTMAX_WIDTH)
 238         {
 239           int i = 0, shift = 0;
 240
 241           do
 242             {
 243               uintmax_t limb = mpz_getlimbn (XBIGNUM (x)->value, i++);
 244               v += limb << shift;
 245               shift += GMP_NUMB_BITS;
 246             }
 247           while (shift < bits);
 248         }
 249     }
 250   return v;
 251 }
 252
 253 /* Yield an upper bound on the buffer size needed to contain a C
 254    string representing the NUM in base BASE.  This includes any
 255    preceding '-' and the terminating null.  */
 256 static ptrdiff_t
 257 mpz_bufsize (mpz_t const num, int base)
 258 {
 259   return mpz_sizeinbase (num, base) + 2;
 260 }
 261 ptrdiff_t
 262 bignum_bufsize (Lisp_Object num, int base)
 263 {
 264   return mpz_bufsize (XBIGNUM (num)->value, base);
 265 }
 266
 267 /* Convert NUM to a nearest double, as opposed to mpz_get_d which
 268    truncates toward zero.  */
 269 double
 270 mpz_get_d_rounded (mpz_t const num)
 271 {
 272   ptrdiff_t size = mpz_bufsize (num, 10);
 273
 274   /* Use mpz_get_d as a shortcut for a bignum so small that rounding
 275      errors cannot occur, which is possible if EMACS_INT (not counting
 276      sign) has fewer bits than a double significand.  */
 277   if (! ((FLT_RADIX == 2 && DBL_MANT_DIG <= FIXNUM_BITS - 1)
 278          || (FLT_RADIX == 16 && DBL_MANT_DIG * 4 <= FIXNUM_BITS - 1))
 279       && size <= DBL_DIG + 2)
 280     return mpz_get_d (num);
 281
 282   USE_SAFE_ALLOCA;
 283   char *buf = SAFE_ALLOCA (size);
 284   mpz_get_str (buf, 10, num);
 285   double result = strtod (buf, NULL);
 286   SAFE_FREE ();
 287   return result;
 288 }
 289
 290 /* Store into BUF (of size SIZE) the value of NUM as a base-BASE string.
 291    If BASE is negative, use upper-case digits in base -BASE.
 292    Return the string's length.
 293    SIZE must equal bignum_bufsize (NUM, abs (BASE)).  */
 294 ptrdiff_t
 295 bignum_to_c_string (char *buf, ptrdiff_t size, Lisp_Object num, int base)
 296 {
 297   eassert (bignum_bufsize (num, abs (base)) == size);
 298   mpz_get_str (buf, base, XBIGNUM (num)->value);
 299   ptrdiff_t n = size - 2;
 300   return !buf[n - 1] ? n - 1 : n + !!buf[n];
 301 }
 302
 303 /* Convert NUM to a base-BASE Lisp string.
 304    If BASE is negative, use upper-case digits in base -BASE.  */
 305
 306 Lisp_Object
 307 bignum_to_string (Lisp_Object num, int base)
 308 {
 309   ptrdiff_t size = bignum_bufsize (num, abs (base));
 310   USE_SAFE_ALLOCA;
 311   char *str = SAFE_ALLOCA (size);
 312   ptrdiff_t len = bignum_to_c_string (str, size, num, base);
 313   Lisp_Object result = make_unibyte_string (str, len);
 314   SAFE_FREE ();
 315   return result;
 316 }
 317
 318 /* Create a bignum by scanning NUM, with digits in BASE.
 319    NUM must consist of an optional '-', a nonempty sequence
 320    of base-BASE digits, and a terminating null byte, and
 321    the represented number must not be in fixnum range.  */
 322
 323 Lisp_Object
 324 make_bignum_str (char const *num, int base)
 325 {
 326   struct Lisp_Bignum *b = ALLOCATE_PSEUDOVECTOR (struct Lisp_Bignum, value,
 327                                                  PVEC_BIGNUM);
 328   mpz_init (b->value);
 329   int check = mpz_set_str (b->value, num, base);
 330   eassert (check == 0);
 331   return make_lisp_ptr (b, Lisp_Vectorlike);
 332 }