bignum.c 8.74 KB
Newer Older
Paul Eggert's avatar
Paul Eggert committed
1 2
/* Big numbers for Emacs.

Paul Eggert's avatar
Paul Eggert committed
3
Copyright 2018-2019 Free Software Foundation, Inc.
Paul Eggert's avatar
Paul Eggert committed
4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25

This file is part of GNU Emacs.

GNU Emacs is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or (at
your option) any later version.

GNU Emacs is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with GNU Emacs.  If not, see <https://www.gnu.org/licenses/>.  */

#include <config.h>

#include "bignum.h"

#include "lisp.h"

Paul Eggert's avatar
Paul Eggert committed
26
#include <math.h>
27 28
#include <stdlib.h>

29 30 31 32 33
/* mpz global temporaries.  Making them global saves the trouble of
   properly using mpz_init and mpz_clear on temporaries even when
   storage is exhausted.  Admittedly this is not ideal.  An mpz value
   in a temporary is made permanent by mpz_swapping it with a bignum's
   value.  Although typically at most two temporaries are needed,
34
   time_arith, rounddiv_q and rounding_driver each need four.  */
35 36 37

mpz_t mpz[4];

Paul Eggert's avatar
Paul Eggert committed
38 39 40 41 42 43 44 45 46 47 48 49
static void *
xrealloc_for_gmp (void *ptr, size_t ignore, size_t size)
{
  return xrealloc (ptr, size);
}

static void
xfree_for_gmp (void *ptr, size_t ignore)
{
  xfree (ptr);
}

50
void
Paul Eggert's avatar
Paul Eggert committed
51
init_bignum (void)
52
{
Paul Eggert's avatar
Paul Eggert committed
53 54 55 56
  eassert (mp_bits_per_limb == GMP_NUMB_BITS);
  integer_width = 1 << 16;
  mp_set_memory_functions (xmalloc, xrealloc_for_gmp, xfree_for_gmp);

57 58 59 60
  for (int i = 0; i < ARRAYELTS (mpz); i++)
    mpz_init (mpz[i]);
}

Paul Eggert's avatar
Paul Eggert committed
61 62 63 64
/* Return the value of the Lisp bignum N, as a double.  */
double
bignum_to_double (Lisp_Object n)
{
65
  return mpz_get_d_rounded (XBIGNUM (n)->value);
Paul Eggert's avatar
Paul Eggert committed
66 67
}

Paul Eggert's avatar
Paul Eggert committed
68 69
/* Return D, converted to a Lisp integer.  Discard any fraction.
   Signal an error if D cannot be converted.  */
Paul Eggert's avatar
Paul Eggert committed
70
Lisp_Object
71
double_to_integer (double d)
Paul Eggert's avatar
Paul Eggert committed
72
{
Paul Eggert's avatar
Paul Eggert committed
73 74
  if (!isfinite (d))
    overflow_error ();
75 76
  mpz_set_d (mpz[0], d);
  return make_integer_mpz ();
Paul Eggert's avatar
Paul Eggert committed
77 78
}

79 80
/* Return a Lisp integer equal to mpz[0], which has BITS bits and which
   must not be in fixnum range.  Set mpz[0] to a junk value.  */
Paul Eggert's avatar
Paul Eggert committed
81
static Lisp_Object
82
make_bignum_bits (size_t bits)
Paul Eggert's avatar
Paul Eggert committed
83 84 85 86
{
  /* The documentation says integer-width should be nonnegative, so
     a single comparison suffices even though 'bits' is unsigned.  */
  if (integer_width < bits)
87
    overflow_error ();
Paul Eggert's avatar
Paul Eggert committed
88

89 90
  struct Lisp_Bignum *b = ALLOCATE_PLAIN_PSEUDOVECTOR (struct Lisp_Bignum,
						       PVEC_BIGNUM);
91 92
  mpz_init (b->value);
  mpz_swap (b->value, mpz[0]);
Paul Eggert's avatar
Paul Eggert committed
93 94 95
  return make_lisp_ptr (b, Lisp_Vectorlike);
}

96 97
/* Return a Lisp integer equal to mpz[0], which must not be in fixnum range.
   Set mpz[0] to a junk value.  */
Paul Eggert's avatar
Paul Eggert committed
98
static Lisp_Object
99
make_bignum (void)
Paul Eggert's avatar
Paul Eggert committed
100
{
101
  return make_bignum_bits (mpz_sizeinbase (mpz[0], 2));
Paul Eggert's avatar
Paul Eggert committed
102 103 104 105 106 107 108
}

/* Return a Lisp integer equal to N, which must not be in fixnum range.  */
Lisp_Object
make_bigint (intmax_t n)
{
  eassert (FIXNUM_OVERFLOW_P (n));
109 110
  mpz_set_intmax (mpz[0], n);
  return make_bignum ();
Paul Eggert's avatar
Paul Eggert committed
111
}
112 113 114 115
Lisp_Object
make_biguint (uintmax_t n)
{
  eassert (FIXNUM_OVERFLOW_P (n));
116 117
  mpz_set_uintmax (mpz[0], n);
  return make_bignum ();
118
}
Paul Eggert's avatar
Paul Eggert committed
119

120 121 122 123 124 125 126 127 128 129
/* Return a Lisp integer equal to -N, which must not be in fixnum range.  */
Lisp_Object
make_neg_biguint (uintmax_t n)
{
  eassert (-MOST_NEGATIVE_FIXNUM < n);
  mpz_set_uintmax (mpz[0], n);
  mpz_neg (mpz[0], mpz[0]);
  return make_bignum ();
}

130 131
/* Return a Lisp integer with value taken from mpz[0].
   Set mpz[0] to a junk value.  */
Paul Eggert's avatar
Paul Eggert committed
132
Lisp_Object
133
make_integer_mpz (void)
Paul Eggert's avatar
Paul Eggert committed
134
{
135
  size_t bits = mpz_sizeinbase (mpz[0], 2);
Paul Eggert's avatar
Paul Eggert committed
136 137 138 139 140 141 142 143

  if (bits <= FIXNUM_BITS)
    {
      EMACS_INT v = 0;
      int i = 0, shift = 0;

      do
	{
144
	  EMACS_INT limb = mpz_getlimbn (mpz[0], i++);
Paul Eggert's avatar
Paul Eggert committed
145 146 147 148 149
	  v += limb << shift;
	  shift += GMP_NUMB_BITS;
	}
      while (shift < bits);

150
      if (mpz_sgn (mpz[0]) < 0)
Paul Eggert's avatar
Paul Eggert committed
151 152 153 154 155 156
	v = -v;

      if (!FIXNUM_OVERFLOW_P (v))
	return make_fixnum (v);
    }

157
  return make_bignum_bits (bits);
Paul Eggert's avatar
Paul Eggert committed
158 159
}

160
/* Set RESULT to V.  This code is for when intmax_t is wider than long.  */
Paul Eggert's avatar
Paul Eggert committed
161 162 163
void
mpz_set_intmax_slow (mpz_t result, intmax_t v)
{
164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183
  int maxlimbs = (INTMAX_WIDTH + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS;
  mp_limb_t *limb = mpz_limbs_write (result, maxlimbs);
  int n = 0;
  uintmax_t u = v;
  bool negative = v < 0;
  if (negative)
    {
      uintmax_t two = 2;
      u = -u & ((two << (UINTMAX_WIDTH - 1)) - 1);
    }

  do
    {
      limb[n++] = u;
      u = GMP_NUMB_BITS < UINTMAX_WIDTH ? u >> GMP_NUMB_BITS : 0;
    }
  while (u != 0);

  mpz_limbs_finish (result, negative ? -n : n);
}
184
void
185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
mpz_set_uintmax_slow (mpz_t result, uintmax_t v)
{
  int maxlimbs = (UINTMAX_WIDTH + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS;
  mp_limb_t *limb = mpz_limbs_write (result, maxlimbs);
  int n = 0;

  do
    {
      limb[n++] = v;
      v = GMP_NUMB_BITS < INTMAX_WIDTH ? v >> GMP_NUMB_BITS : 0;
    }
  while (v != 0);

  mpz_limbs_finish (result, n);
}

201 202 203 204
/* If Z fits into *PI, store its value there and return true.
   Return false otherwise.  */
bool
mpz_to_intmax (mpz_t const z, intmax_t *pi)
205
{
206 207
  ptrdiff_t bits = mpz_sizeinbase (z, 2);
  bool negative = mpz_sgn (z) < 0;
208 209 210 211 212 213 214 215

  if (bits < INTMAX_WIDTH)
    {
      intmax_t v = 0;
      int i = 0, shift = 0;

      do
	{
216
	  intmax_t limb = mpz_getlimbn (z, i++);
217 218 219 220 221
	  v += limb << shift;
	  shift += GMP_NUMB_BITS;
	}
      while (shift < bits);

222 223 224 225 226 227 228 229
      *pi = negative ? -v : v;
      return true;
    }
  if (bits == INTMAX_WIDTH && INTMAX_MIN < -INTMAX_MAX && negative
      && mpz_scan1 (z, 0) == INTMAX_WIDTH - 1)
    {
      *pi = INTMAX_MIN;
      return true;
230
    }
231
  return false;
232
}
233 234
bool
mpz_to_uintmax (mpz_t const z, uintmax_t *pi)
235
{
236 237 238 239 240 241
  if (mpz_sgn (z) < 0)
    return false;
  ptrdiff_t bits = mpz_sizeinbase (z, 2);
  if (UINTMAX_WIDTH < bits)
    return false;

242
  uintmax_t v = 0;
243 244 245
  int i = 0, shift = 0;

  do
246
    {
247 248 249
      uintmax_t limb = mpz_getlimbn (z, i++);
      v += limb << shift;
      shift += GMP_NUMB_BITS;
250
    }
251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269
  while (shift < bits);

  *pi = v;
  return true;
}

/* Return the value of the bignum X if it fits, 0 otherwise.
   A bignum cannot be zero, so 0 indicates failure reliably.  */
intmax_t
bignum_to_intmax (Lisp_Object x)
{
  intmax_t i;
  return mpz_to_intmax (XBIGNUM (x)->value, &i) ? i : 0;
}
uintmax_t
bignum_to_uintmax (Lisp_Object x)
{
  uintmax_t i;
  return mpz_to_uintmax (XBIGNUM (x)->value, &i) ? i : 0;
Paul Eggert's avatar
Paul Eggert committed
270 271
}

272
/* Yield an upper bound on the buffer size needed to contain a C
273
   string representing the NUM in base BASE.  This includes any
274
   preceding '-' and the terminating NUL.  */
275 276 277 278 279
static ptrdiff_t
mpz_bufsize (mpz_t const num, int base)
{
  return mpz_sizeinbase (num, base) + 2;
}
280 281 282
ptrdiff_t
bignum_bufsize (Lisp_Object num, int base)
{
283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306
  return mpz_bufsize (XBIGNUM (num)->value, base);
}

/* Convert NUM to a nearest double, as opposed to mpz_get_d which
   truncates toward zero.  */
double
mpz_get_d_rounded (mpz_t const num)
{
  ptrdiff_t size = mpz_bufsize (num, 10);

  /* Use mpz_get_d as a shortcut for a bignum so small that rounding
     errors cannot occur, which is possible if EMACS_INT (not counting
     sign) has fewer bits than a double significand.  */
  if (! ((FLT_RADIX == 2 && DBL_MANT_DIG <= FIXNUM_BITS - 1)
	 || (FLT_RADIX == 16 && DBL_MANT_DIG * 4 <= FIXNUM_BITS - 1))
      && size <= DBL_DIG + 2)
    return mpz_get_d (num);

  USE_SAFE_ALLOCA;
  char *buf = SAFE_ALLOCA (size);
  mpz_get_str (buf, 10, num);
  double result = strtod (buf, NULL);
  SAFE_FREE ();
  return result;
307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323
}

/* Store into BUF (of size SIZE) the value of NUM as a base-BASE string.
   If BASE is negative, use upper-case digits in base -BASE.
   Return the string's length.
   SIZE must equal bignum_bufsize (NUM, abs (BASE)).  */
ptrdiff_t
bignum_to_c_string (char *buf, ptrdiff_t size, Lisp_Object num, int base)
{
  eassert (bignum_bufsize (num, abs (base)) == size);
  mpz_get_str (buf, base, XBIGNUM (num)->value);
  ptrdiff_t n = size - 2;
  return !buf[n - 1] ? n - 1 : n + !!buf[n];
}

/* Convert NUM to a base-BASE Lisp string.
   If BASE is negative, use upper-case digits in base -BASE.  */
Paul Eggert's avatar
Paul Eggert committed
324 325 326 327

Lisp_Object
bignum_to_string (Lisp_Object num, int base)
{
328
  ptrdiff_t size = bignum_bufsize (num, abs (base));
Paul Eggert's avatar
Paul Eggert committed
329
  USE_SAFE_ALLOCA;
330 331 332
  char *str = SAFE_ALLOCA (size);
  ptrdiff_t len = bignum_to_c_string (str, size, num, base);
  Lisp_Object result = make_unibyte_string (str, len);
Paul Eggert's avatar
Paul Eggert committed
333 334 335 336 337 338
  SAFE_FREE ();
  return result;
}

/* Create a bignum by scanning NUM, with digits in BASE.
   NUM must consist of an optional '-', a nonempty sequence
339
   of base-BASE digits, and a terminating NUL byte, and
Paul Eggert's avatar
Paul Eggert committed
340 341 342 343 344
   the represented number must not be in fixnum range.  */

Lisp_Object
make_bignum_str (char const *num, int base)
{
345 346
  struct Lisp_Bignum *b = ALLOCATE_PLAIN_PSEUDOVECTOR (struct Lisp_Bignum,
						       PVEC_BIGNUM);
Paul Eggert's avatar
Paul Eggert committed
347 348 349 350 351
  mpz_init (b->value);
  int check = mpz_set_str (b->value, num, base);
  eassert (check == 0);
  return make_lisp_ptr (b, Lisp_Vectorlike);
}