mirror/guile
mirror/guile
1
Fork 0
mirror of https://git.savannah.gnu.org/git/guile.git synced 2025-04-29 19:30:36 +02:00
Code Activity
guile/libguile/integers.c
Michael Gran 1ddc4eb964 Fixes arg type for scm_integer_from_mpz 
The definition and the declaration for scm_integer_from_mpz
do not match

* libguile/integers.c (scm_integer_from_mpz): takes const mpz_t arg
  No callers need to be changed.
2022-10-15 10:35:54 -07:00

3424 lines
74 KiB
C
Raw Permalink Blame History

This file contains invisible Unicode characters

This file contains invisible Unicode characters that are indistinguishable to humans but may be processed differently by a computer. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

/* Copyright 1995-2016,2018-2022
Free Software Foundation, Inc.
This file is part of Guile.
Guile is free software: you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published
by the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
Guile 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 Lesser General Public
License for more details.
You should have received a copy of the GNU Lesser General Public
License along with Guile. If not, see
<https://www.gnu.org/licenses/>. */
#ifdef HAVE_CONFIG_H
# include <config.h>
#endif
#include <math.h>
#include <stdlib.h>
#include <stdint.h>
#include <stdio.h>
#include <string.h>
#include <verify.h>
#include "boolean.h"
#include "numbers.h"
#include "strings.h"
#include "integers.h"
/* Some functions that use GMP's mpn functions assume that a
non-negative fixnum will always fit in a 'mp_limb_t'. */
verify (SCM_MOST_POSITIVE_FIXNUM <= (mp_limb_t) -1);
#define NLIMBS_MAX (SSIZE_MAX / sizeof(mp_limb_t))
#ifndef NDEBUG
#define ASSERT(x) \
do { \
if (!(x)) \
{ \
fprintf (stderr, "%s:%d: assertion failed\n", __FILE__, __LINE__); \
abort(); \
} \
} while (0)
#else
#define ASSERT(x) do { } while (0)
#endif
struct scm_bignum
{
scm_t_bits tag;
/* FIXME: In Guile 3.2, replace this union with just a "size" member.
Digits are always allocated inline. */
union {
mpz_t mpz;
struct {
int zero;
int size;
mp_limb_t *limbs;
} z;
} u;
mp_limb_t limbs[];
};
static int
bignum_size (struct scm_bignum *z)
{
return z->u.z.size;
}
static int
bignum_is_negative (struct scm_bignum *z)
{
return bignum_size (z) < 0;
}
static int
bignum_is_positive (struct scm_bignum *z)
{
return bignum_size (z) > 0;
}
static size_t
bignum_limb_count (struct scm_bignum *z)
{
return bignum_is_negative (z) ? -bignum_size (z) : bignum_size (z);
}
static mp_limb_t*
bignum_limbs (struct scm_bignum *z)
{
// FIXME: In the future we can just return z->limbs.
return z->u.z.limbs;
}
static inline unsigned long
long_magnitude (long l)
{
unsigned long mag = l;
return l < 0 ? ~mag + 1 : mag;
}
static inline long
negative_long (unsigned long mag)
{
ASSERT (mag <= (unsigned long) LONG_MIN);
return ~mag + 1;
}
static inline int64_t
negative_int64 (uint64_t mag)
{
ASSERT (mag <= (uint64_t) INT64_MIN);
return ~mag + 1;
}
static inline uint64_t
int64_magnitude (int64_t i)
{
uint64_t mag = i;
if (i < 0)
mag = ~mag + 1;
return mag;
}
static inline scm_t_bits
inum_magnitude (scm_t_inum i)
{
scm_t_bits mag = i;
if (i < 0)
mag = ~mag + 1;
return mag;
}
static struct scm_bignum *
allocate_bignum (size_t nlimbs)
{
ASSERT (nlimbs <= (size_t)INT_MAX);
ASSERT (nlimbs <= NLIMBS_MAX);
size_t size = sizeof (struct scm_bignum) + nlimbs * sizeof(mp_limb_t);
struct scm_bignum *z = scm_gc_malloc_pointerless (size, "bignum");
z->tag = scm_tc16_big;
z->u.z.zero = 0;
z->u.z.size = nlimbs;
z->u.z.limbs = z->limbs;
// _mp_alloc == 0 means GMP will never try to free this memory.
ASSERT (z->u.mpz[0]._mp_alloc == 0);
// Our "size" field should alias the mpz's _mp_size field.
ASSERT (z->u.mpz[0]._mp_size == nlimbs);
// Limbs are always allocated inline.
ASSERT (z->u.mpz[0]._mp_d == z->limbs);
// z->limbs left uninitialized.
return z;
}
static struct scm_bignum *
bignum_trim1 (struct scm_bignum *z)
{
ASSERT (z->u.z.size > 0);
z->u.z.size -= (z->limbs[z->u.z.size - 1] == 0);
return z;
}
static struct scm_bignum *
bignum_trimn (struct scm_bignum *z)
{
ASSERT (z->u.z.size > 0);
while (z->u.z.size > 0 && z->limbs[z->u.z.size - 1] == 0)
z->u.z.size--;
return z;
}
static struct scm_bignum *
negate_bignum (struct scm_bignum *z)
{
z->u.z.size = -z->u.z.size;
return z;
}
static struct scm_bignum *
bignum_negate_if (int negate, struct scm_bignum *z)
{
return negate ? negate_bignum (z) : z;
}
static struct scm_bignum *
make_bignum_0 (void)
{
return allocate_bignum (0);
}
static struct scm_bignum *
make_bignum_1 (int is_negative, mp_limb_t limb)
{
struct scm_bignum *z = allocate_bignum (1);
z->limbs[0] = limb;
return is_negative ? negate_bignum(z) : z;
}
static struct scm_bignum *
make_bignum_2 (int is_negative, mp_limb_t lo, mp_limb_t hi)
{
struct scm_bignum *z = allocate_bignum (2);
z->limbs[0] = lo;
z->limbs[1] = hi;
return is_negative ? negate_bignum(z) : z;
}
static struct scm_bignum *
make_bignum_from_uint64 (uint64_t val)
{
#if SCM_SIZEOF_LONG == 4
if (val > UINT32_MAX)
return make_bignum_2 (0, val, val >> 32);
#endif
return val == 0 ? make_bignum_0 () : make_bignum_1 (0, val);
}
static struct scm_bignum *
make_bignum_from_int64 (int64_t val)
{
return val < 0
? negate_bignum (make_bignum_from_uint64 (int64_magnitude (val)))
: make_bignum_from_uint64 (val);
}
static struct scm_bignum *
ulong_to_bignum (unsigned long u)
{
return u == 0 ? make_bignum_0 () : make_bignum_1 (0, u);
};
static struct scm_bignum *
long_to_bignum (long i)
{
if (i > 0)
return ulong_to_bignum (i);
return i == 0 ? make_bignum_0 () : make_bignum_1 (1, long_magnitude (i));
};
static inline SCM
scm_from_bignum (struct scm_bignum *x)
{
return SCM_PACK (x);
}
static SCM
long_to_scm (long i)
{
if (SCM_FIXABLE (i))
return SCM_I_MAKINUM (i);
return scm_from_bignum (long_to_bignum (i));
}
static SCM
ulong_to_scm (unsigned long i)
{
if (SCM_POSFIXABLE (i))
return SCM_I_MAKINUM (i);
return scm_from_bignum (ulong_to_bignum (i));
}
static struct scm_bignum *
clone_bignum (struct scm_bignum *z)
{
struct scm_bignum *ret = allocate_bignum (bignum_limb_count (z));
mpn_copyi (bignum_limbs (ret), bignum_limbs (z), bignum_limb_count (z));
return bignum_is_negative (z) ? negate_bignum (ret) : ret;
}
static void
alias_bignum_to_mpz (struct scm_bignum *z, mpz_ptr mpz)
{
// No need to clear this mpz.
mpz->_mp_alloc = 0;
mpz->_mp_size = bignum_size (z);
// Gotta be careful to keep z alive.
mpz->_mp_d = bignum_limbs (z);
}
static struct scm_bignum *
make_bignum_from_mpz (mpz_srcptr mpz)
{
size_t nlimbs = mpz_size (mpz);
struct scm_bignum *ret = allocate_bignum (nlimbs);
mpn_copyi (bignum_limbs (ret), mpz_limbs_read (mpz), nlimbs);
return mpz_sgn (mpz) < 0 ? negate_bignum (ret) : ret;
}
static SCM
normalize_bignum (struct scm_bignum *z)
{
switch (bignum_size (z))
{
case -1:
if (bignum_limbs (z)[0] <= inum_magnitude (SCM_MOST_NEGATIVE_FIXNUM))
return SCM_I_MAKINUM (negative_long (bignum_limbs (z)[0]));
break;
case 0:
return SCM_INUM0;
case 1:
if (bignum_limbs (z)[0] <= SCM_MOST_POSITIVE_FIXNUM)
return SCM_I_MAKINUM (bignum_limbs (z)[0]);
break;
default:
break;
}
return scm_from_bignum (z);
}
static SCM
take_mpz (mpz_ptr mpz)
{
SCM ret;
if (mpz_fits_slong_p (mpz))
ret = long_to_scm (mpz_get_si (mpz));
else
ret = scm_from_bignum (make_bignum_from_mpz (mpz));
mpz_clear (mpz);
return ret;
}
static int
long_sign (long l)
{
if (l < 0) return -1;
if (l == 0) return 0;
return 1;
}
static int
negative_uint64_to_int64 (uint64_t magnitude, int64_t *val)
{
if (magnitude > int64_magnitude (INT64_MIN))
return 0;
*val = negative_int64 (magnitude);
return 1;
}
static int
positive_uint64_to_int64 (uint64_t magnitude, int64_t *val)
{
if (magnitude > INT64_MAX)
return 0;
*val = magnitude;
return 1;
}
static int
bignum_to_int64 (struct scm_bignum *z, int64_t *val)
{
switch (bignum_size (z))
{
#if SCM_SIZEOF_LONG == 4
case -2:
{
uint64_t mag = bignum_limbs (z)[0];
mag |= ((uint64_t) bignum_limbs (z)[1]) << 32;
return negative_uint64_to_int64 (mag, val);
}
#endif
case -1:
return negative_uint64_to_int64 (bignum_limbs (z)[0], val);
case 0:
*val = 0;
return 1;
case 1:
return positive_uint64_to_int64 (bignum_limbs (z)[0], val);
#if SCM_SIZEOF_LONG == 4
case 2:
{
uint64_t mag = bignum_limbs (z)[0];
mag |= ((uint64_t) bignum_limbs (z)[1]) << 32;
return positive_uint64_to_int64 (mag, val);
}
#endif
default:
return 0;
}
}
static int
bignum_to_uint64 (struct scm_bignum *z, uint64_t *val)
{
switch (bignum_size (z))
{
case 0:
*val = 0;
return 1;
case 1:
*val = bignum_limbs (z)[0];
return 1;
#if SCM_SIZEOF_LONG == 4
case 2:
{
uint64_t mag = bignum_limbs (z)[0];
mag |= ((uint64_t) bignum_limbs (z)[1]) << 32;
*val = mag;
return 1;
}
#endif
default:
return 0;
}
}
#if SCM_SIZEOF_LONG == 4
static int
negative_uint32_to_int32 (uint32_t magnitude, int32_t *val)
{
if (magnitude > long_magnitude (INT32_MIN))
return 0;
*val = negative_long (magnitude);
return 1;
}
static int
positive_uint32_to_int32 (uint32_t magnitude, int32_t *val)
{
if (magnitude > INT32_MAX)
return 0;
*val = magnitude;
return 1;
}
static int
bignum_to_int32 (struct scm_bignum *z, int32_t *val)
{
switch (bignum_size (z))
{
case -1:
return negative_uint32_to_int32 (bignum_limbs (z)[0], val);
case 0:
*val = 0;
return 1;
case 1:
return positive_uint32_to_int32 (bignum_limbs (z)[0], val);
default:
return 0;
}
}
static int
bignum_to_uint32 (struct scm_bignum *z, uint32_t *val)
{
switch (bignum_size (z))
{
case 0:
*val = 0;
return 1;
case 1:
*val = bignum_limbs (z)[0];
return 1;
default:
return 0;
}
}
#endif
static int
bignum_cmp_long (struct scm_bignum *z, long l)
{
switch (bignum_size (z))
{
case -1:
if (l >= 0)
return -1;
return long_sign (long_magnitude (l) - bignum_limbs (z)[0]);
case 0:
return long_sign (l);
case 1:
if (l <= 0)
return 1;
return long_sign (bignum_limbs (z)[0] - (unsigned long) l);
default:
return long_sign (bignum_size (z));
}
}
SCM
scm_integer_from_mpz (const mpz_t mpz)
{
return normalize_bignum (make_bignum_from_mpz (mpz));
}
int
scm_is_integer_odd_i (scm_t_inum i)
{
return i & 1;
}
int
scm_is_integer_odd_z (struct scm_bignum *z)
{
return bignum_limbs (z)[0] & 1;
}
SCM
scm_integer_abs_i (scm_t_inum i)
{
if (i >= 0)
return SCM_I_MAKINUM (i);
return ulong_to_scm (long_magnitude (i));
}
SCM
scm_integer_abs_z (struct scm_bignum *z)
{
if (!bignum_is_negative (z))
return scm_from_bignum (z);
return scm_integer_negate_z (z);
}
SCM
scm_integer_floor_quotient_ii (scm_t_inum x, scm_t_inum y)
{
if (y > 0)
{
if (x < 0)
x = x - y + 1;
}
else if (y == 0)
scm_num_overflow ("floor-quotient");
else if (x > 0)
x = x - y - 1;
scm_t_inum q = x / y;
return long_to_scm (q);
}
SCM
scm_integer_floor_quotient_iz (scm_t_inum x, struct scm_bignum *y)
{
if (x == 0 || ((x < 0) == bignum_is_negative (y)))
return SCM_INUM0;
return SCM_I_MAKINUM (-1);
}
SCM
scm_integer_floor_quotient_zi (struct scm_bignum *x, scm_t_inum y)
{
if (y == 0)
scm_num_overflow ("floor-quotient");
else if (y == 1)
return scm_from_bignum (x);
mpz_t zx, q;
alias_bignum_to_mpz (x, zx);
mpz_init (q);
if (y > 0)
mpz_fdiv_q_ui (q, zx, y);
else
{
mpz_cdiv_q_ui (q, zx, -y);
mpz_neg (q, q);
}
scm_remember_upto_here_1 (x);
return take_mpz (q);
}
SCM
scm_integer_floor_quotient_zz (struct scm_bignum *x, struct scm_bignum *y)
{
mpz_t zx, zy, q;
alias_bignum_to_mpz (x, zx);
alias_bignum_to_mpz (y, zy);
mpz_init (q);
mpz_fdiv_q (q, zx, zy);
scm_remember_upto_here_2 (x, y);
return take_mpz (q);
}
SCM
scm_integer_floor_remainder_ii (scm_t_inum x, scm_t_inum y)
{
if (y == 0)
scm_num_overflow ("floor-remainder");
scm_t_inum r = x % y;
int needs_adjustment = (y > 0) ? (r < 0) : (r > 0);
if (needs_adjustment)
r += y;
return SCM_I_MAKINUM (r);
}
SCM
scm_integer_floor_remainder_iz (scm_t_inum x, struct scm_bignum *y)
{
if (bignum_is_positive (y))
{
if (x < 0)
{
mpz_t r, zy;
mpz_init (r);
alias_bignum_to_mpz (y, zy);
mpz_sub_ui (r, zy, -x);
scm_remember_upto_here_1 (y);
return take_mpz (r);
}
else
return SCM_I_MAKINUM (x);
}
else if (x <= 0)
return SCM_I_MAKINUM (x);
else
{
mpz_t r, zy;
mpz_init (r);
alias_bignum_to_mpz (y, zy);
mpz_add_ui (r, zy, x);
scm_remember_upto_here_1 (y);
return take_mpz (r);
}
}
SCM
scm_integer_floor_remainder_zi (struct scm_bignum *x, scm_t_inum y)
{
if (y == 0)
scm_num_overflow ("floor-remainder");
else
{
scm_t_inum r;
mpz_t zx;
alias_bignum_to_mpz (x, zx);
if (y > 0)
r = mpz_fdiv_ui (zx, y);
else
r = -mpz_cdiv_ui (zx, -y);
scm_remember_upto_here_1 (x);
return SCM_I_MAKINUM (r);
}
}
SCM
scm_integer_floor_remainder_zz (struct scm_bignum *x, struct scm_bignum *y)
{
mpz_t zx, zy, r;
alias_bignum_to_mpz (x, zx);
alias_bignum_to_mpz (y, zy);
mpz_init (r);
mpz_fdiv_r (r, zx, zy);
scm_remember_upto_here_2 (x, y);
return take_mpz (r);
}
void
scm_integer_floor_divide_ii (scm_t_inum x, scm_t_inum y, SCM *qp, SCM *rp)
{
if (y == 0)
scm_num_overflow ("floor-divide");
scm_t_inum q = x / y;
scm_t_inum r = x % y;
int needs_adjustment = (y > 0) ? (r < 0) : (r > 0);
if (needs_adjustment)
{
r += y;
q--;
}
*qp = long_to_scm (q);
*rp = SCM_I_MAKINUM (r);
}
void
scm_integer_floor_divide_iz (scm_t_inum x, struct scm_bignum *y, SCM *qp, SCM *rp)
{
if (bignum_is_positive (y))
{
if (x < 0)
{
mpz_t zy, r;
alias_bignum_to_mpz (y, zy);
mpz_init (r);
mpz_sub_ui (r, zy, -x);
scm_remember_upto_here_1 (y);
*qp = SCM_I_MAKINUM (-1);
*rp = take_mpz (r);
}
else
{
*qp = SCM_INUM0;
*rp = SCM_I_MAKINUM (x);
}
}
else if (x <= 0)
{
*qp = SCM_INUM0;
*rp = SCM_I_MAKINUM (x);
}
else
{
mpz_t zy, r;
alias_bignum_to_mpz (y, zy);
mpz_init (r);
mpz_add_ui (r, zy, x);
scm_remember_upto_here_1 (y);
*qp = SCM_I_MAKINUM (-1);
*rp = take_mpz (r);
}
}
void
scm_integer_floor_divide_zi (struct scm_bignum *x, scm_t_inum y, SCM *qp, SCM *rp)
{
if (y == 0)
scm_num_overflow ("floor-divide");
mpz_t zx, q, r;
alias_bignum_to_mpz (x, zx);
mpz_init (q);
mpz_init (r);
if (y > 0)
mpz_fdiv_qr_ui (q, r, zx, y);
else
{
mpz_cdiv_qr_ui (q, r, zx, -y);
mpz_neg (q, q);
}
scm_remember_upto_here_1 (x);
*qp = take_mpz (q);
*rp = take_mpz (r);
}
void
scm_integer_floor_divide_zz (struct scm_bignum *x, struct scm_bignum *y, SCM *qp, SCM *rp)
{
mpz_t zx, zy, q, r;
mpz_init (q);
mpz_init (r);
alias_bignum_to_mpz (x, zx);
alias_bignum_to_mpz (y, zy);
mpz_fdiv_qr (q, r, zx, zy);
scm_remember_upto_here_2 (x, y);
*qp = take_mpz (q);
*rp = take_mpz (r);
}
SCM
scm_integer_ceiling_quotient_ii (scm_t_inum x, scm_t_inum y)
{
if (y == 0)
scm_num_overflow ("ceiling-quotient");
if (y > 0)
{
if (x >= 0)
x = x + y - 1;
}
else if (x < 0)
x = x + y + 1;
scm_t_inum q = x / y;
return long_to_scm (q);
}
SCM
scm_integer_ceiling_quotient_iz (scm_t_inum x, struct scm_bignum *y)
{
if (bignum_is_positive (y))
{
if (x > 0)
return SCM_INUM1;
else if (x == SCM_MOST_NEGATIVE_FIXNUM &&
bignum_cmp_long (y, -SCM_MOST_NEGATIVE_FIXNUM) == 0)
{
/* Special case: x == fixnum-min && y == abs (fixnum-min) */
scm_remember_upto_here_1 (y);
return SCM_I_MAKINUM (-1);
}
else
return SCM_INUM0;
}
else if (x >= 0)
return SCM_INUM0;
else
return SCM_INUM1;
}
SCM
scm_integer_ceiling_quotient_zi (struct scm_bignum *x, scm_t_inum y)
{
if (y == 0)
scm_num_overflow ("ceiling-quotient");
else if (y == 1)
return scm_from_bignum (x);
else
{
mpz_t q, zx;
mpz_init (q);
alias_bignum_to_mpz (x, zx);
if (y > 0)
mpz_cdiv_q_ui (q, zx, y);
else
{
mpz_fdiv_q_ui (q, zx, -y);
mpz_neg (q, q);
}
scm_remember_upto_here_1 (x);
return take_mpz (q);
}
}
SCM
scm_integer_ceiling_quotient_zz (struct scm_bignum *x, struct scm_bignum *y)
{
mpz_t q, zx, zy;
mpz_init (q);
alias_bignum_to_mpz (x, zx);
alias_bignum_to_mpz (y, zy);
mpz_cdiv_q (q, zx, zy);
scm_remember_upto_here_2 (x, y);
return take_mpz (q);
}
SCM
scm_integer_ceiling_remainder_ii (scm_t_inum x, scm_t_inum y)
{
if (y == 0)
scm_num_overflow ("ceiling-remainder");
scm_t_inum r = x % y;
int needs_adjustment = (y > 0) ? (r > 0) : (r < 0);
if (needs_adjustment)
r -= y;
return SCM_I_MAKINUM (r);
}
SCM
scm_integer_ceiling_remainder_iz (scm_t_inum x, struct scm_bignum *y)
{
if (bignum_is_positive (y))
{
if (x > 0)
{
mpz_t r, zy;
mpz_init (r);
alias_bignum_to_mpz (y, zy);
mpz_sub_ui (r, zy, x);
scm_remember_upto_here_1 (y);
mpz_neg (r, r);
return take_mpz (r);
}
else if (x == SCM_MOST_NEGATIVE_FIXNUM &&
bignum_cmp_long (y, -SCM_MOST_NEGATIVE_FIXNUM) == 0)
{
/* Special case: x == fixnum-min && y == abs (fixnum-min) */
scm_remember_upto_here_1 (y);
return SCM_INUM0;
}
else
return SCM_I_MAKINUM (x);
}
else if (x >= 0)
return SCM_I_MAKINUM (x);
else
{
mpz_t r, zy;
mpz_init (r);
alias_bignum_to_mpz (y, zy);
mpz_add_ui (r, zy, -x);
scm_remember_upto_here_1 (y);
mpz_neg (r, r);
return take_mpz (r);
}
}
SCM
scm_integer_ceiling_remainder_zi (struct scm_bignum *x, scm_t_inum y)
{
if (y == 0)
scm_num_overflow ("ceiling-remainder");
else
{
mpz_t zx;
alias_bignum_to_mpz (x, zx);
scm_t_inum r;
if (y > 0)
r = -mpz_cdiv_ui (zx, y);
else
r = mpz_fdiv_ui (zx, -y);
scm_remember_upto_here_1 (x);
return SCM_I_MAKINUM (r);
}
}
SCM
scm_integer_ceiling_remainder_zz (struct scm_bignum *x, struct scm_bignum *y)
{
mpz_t r, zx, zy;
mpz_init (r);
alias_bignum_to_mpz (x, zx);
alias_bignum_to_mpz (y, zy);
mpz_cdiv_r (r, zx, zy);
scm_remember_upto_here_2 (x, y);
return take_mpz (r);
}
void
scm_integer_ceiling_divide_ii (scm_t_inum x, scm_t_inum y, SCM *qp, SCM *rp)
{
if (y == 0)
scm_num_overflow ("ceiling-divide");
else
{
scm_t_inum q = x / y;
scm_t_inum r = x % y;
int needs_adjustment;
if (y > 0)
needs_adjustment = (r > 0);
else
needs_adjustment = (r < 0);
if (needs_adjustment)
{
r -= y;
q++;
}
*qp = long_to_scm (q);
*rp = SCM_I_MAKINUM (r);
}
}
void
scm_integer_ceiling_divide_iz (scm_t_inum x, struct scm_bignum *y, SCM *qp, SCM *rp)
{
if (bignum_is_positive (y))
{
if (x > 0)
{
mpz_t r, zy;
mpz_init (r);
alias_bignum_to_mpz (y, zy);
mpz_sub_ui (r, zy, x);
scm_remember_upto_here_1 (y);
mpz_neg (r, r);
*qp = SCM_INUM1;
*rp = take_mpz (r);
}
else if (x == SCM_MOST_NEGATIVE_FIXNUM &&
bignum_cmp_long (y, -SCM_MOST_NEGATIVE_FIXNUM) == 0)
{
/* Special case: x == fixnum-min && y == abs (fixnum-min) */
scm_remember_upto_here_1 (y);
*qp = SCM_I_MAKINUM (-1);
*rp = SCM_INUM0;
}
else
{
*qp = SCM_INUM0;
*rp = SCM_I_MAKINUM (x);
}
}
else if (x >= 0)
{
*qp = SCM_INUM0;
*rp = SCM_I_MAKINUM (x);
}
else
{
mpz_t r, zy;
mpz_init (r);
alias_bignum_to_mpz (y, zy);
mpz_add_ui (r, zy, -x);
scm_remember_upto_here_1 (y);
mpz_neg (r, r);
*qp = SCM_INUM1;
*rp = take_mpz (r);
}
}
void
scm_integer_ceiling_divide_zi (struct scm_bignum *x, scm_t_inum y, SCM *qp, SCM *rp)
{
if (y == 0)
scm_num_overflow ("ceiling-divide");
else
{
mpz_t q, r, zx;
mpz_init (q);
mpz_init (r);
alias_bignum_to_mpz (x, zx);
if (y > 0)
mpz_cdiv_qr_ui (q, r, zx, y);
else
{
mpz_fdiv_qr_ui (q, r, zx, -y);
mpz_neg (q, q);
}
scm_remember_upto_here_1 (x);
*qp = take_mpz (q);
*rp = take_mpz (r);
}
}
void
scm_integer_ceiling_divide_zz (struct scm_bignum *x, struct scm_bignum *y, SCM *qp, SCM *rp)
{
mpz_t q, r, zx, zy;
mpz_init (q);
mpz_init (r);
alias_bignum_to_mpz (x, zx);
alias_bignum_to_mpz (y, zy);
mpz_cdiv_qr (q, r, zx, zy);
scm_remember_upto_here_2 (x, y);
*qp = take_mpz (q);
*rp = take_mpz (r);
}
SCM
scm_integer_truncate_quotient_ii (scm_t_inum x, scm_t_inum y)
{
if (y == 0)
scm_num_overflow ("truncate-quotient");
else
{
scm_t_inum q = x / y;
return long_to_scm (q);
}
}
SCM
scm_integer_truncate_quotient_iz (scm_t_inum x, struct scm_bignum *y)
{
if (x == SCM_MOST_NEGATIVE_FIXNUM &&
bignum_cmp_long (y, -SCM_MOST_NEGATIVE_FIXNUM) == 0)
{
/* Special case: x == fixnum-min && y == abs (fixnum-min) */
scm_remember_upto_here_1 (y);
return SCM_I_MAKINUM (-1);
}
else
return SCM_INUM0;
}
SCM
scm_integer_truncate_quotient_zi (struct scm_bignum *x, scm_t_inum y)
{
if (y == 0)
scm_num_overflow ("truncate-quotient");
else if (y == 1)
return scm_from_bignum (x);
else
{
mpz_t q, zx;
mpz_init (q);
alias_bignum_to_mpz (x, zx);
if (y > 0)
mpz_tdiv_q_ui (q, zx, y);
else
{
mpz_tdiv_q_ui (q, zx, -y);
mpz_neg (q, q);
}
scm_remember_upto_here_1 (x);
return take_mpz (q);
}
}
SCM
scm_integer_truncate_quotient_zz (struct scm_bignum *x, struct scm_bignum *y)
{
mpz_t q, zx, zy;
mpz_init (q);
alias_bignum_to_mpz (x, zx);
alias_bignum_to_mpz (y, zy);
mpz_tdiv_q (q, zx, zy);
scm_remember_upto_here_2 (x, y);
return take_mpz (q);
}
SCM
scm_integer_truncate_remainder_ii (scm_t_inum x, scm_t_inum y)
{
if (y == 0)
scm_num_overflow ("truncate-remainder");
else
{
scm_t_inum q = x % y;
return long_to_scm (q);
}
}
SCM
scm_integer_truncate_remainder_iz (scm_t_inum x, struct scm_bignum *y)
{
if (x == SCM_MOST_NEGATIVE_FIXNUM &&
bignum_cmp_long (y, -SCM_MOST_NEGATIVE_FIXNUM) == 0)
{
/* Special case: x == fixnum-min && y == abs (fixnum-min) */
scm_remember_upto_here_1 (y);
return SCM_INUM0;
}
else
return SCM_I_MAKINUM (x);
}
SCM
scm_integer_truncate_remainder_zi (struct scm_bignum *x, scm_t_inum y)
{
if (y == 0)
scm_num_overflow ("truncate-remainder");
else
{
mpz_t zx;
alias_bignum_to_mpz (x, zx);
scm_t_inum r = mpz_tdiv_ui (zx, (y > 0) ? y : -y) * mpz_sgn (zx);
scm_remember_upto_here_1 (x);
return SCM_I_MAKINUM (r);
}
}
SCM
scm_integer_truncate_remainder_zz (struct scm_bignum *x, struct scm_bignum *y)
{
mpz_t r, zx, zy;
mpz_init (r);
alias_bignum_to_mpz (x, zx);
alias_bignum_to_mpz (y, zy);
mpz_tdiv_r (r, zx, zy);
scm_remember_upto_here_2 (x, y);
return take_mpz (r);
}
void
scm_integer_truncate_divide_ii (scm_t_inum x, scm_t_inum y, SCM *qp, SCM *rp)
{
if (y == 0)
scm_num_overflow ("truncate-divide");
else
{
scm_t_inum q = x / y;
scm_t_inum r = x % y;
*qp = long_to_scm (q);
*rp = SCM_I_MAKINUM (r);
}
}
void
scm_integer_truncate_divide_iz (scm_t_inum x, struct scm_bignum *y, SCM *qp, SCM *rp)
{
if (x == SCM_MOST_NEGATIVE_FIXNUM &&
bignum_cmp_long (y, -SCM_MOST_NEGATIVE_FIXNUM) == 0)
{
/* Special case: x == fixnum-min && y == abs (fixnum-min) */
scm_remember_upto_here_1 (y);
*qp = SCM_I_MAKINUM (-1);
*rp = SCM_INUM0;
}
else
{
*qp = SCM_INUM0;
*rp = SCM_I_MAKINUM (x);
}
}
void
scm_integer_truncate_divide_zi (struct scm_bignum *x, scm_t_inum y, SCM *qp, SCM *rp)
{
if (y == 0)
scm_num_overflow ("truncate-divide");
else
{
mpz_t q, zx;
mpz_init (q);
alias_bignum_to_mpz (x, zx);
scm_t_inum r;
if (y > 0)
r = mpz_tdiv_q_ui (q, zx, y);
else
{
r = mpz_tdiv_q_ui (q, zx, -y);
mpz_neg (q, q);
}
r *= mpz_sgn (zx);
scm_remember_upto_here_1 (x);
*qp = take_mpz (q);
*rp = SCM_I_MAKINUM (r);
}
}
void
scm_integer_truncate_divide_zz (struct scm_bignum *x, struct scm_bignum *y, SCM *qp, SCM *rp)
{
mpz_t q, r, zx, zy;
mpz_init (q);
mpz_init (r);
alias_bignum_to_mpz (x, zx);
alias_bignum_to_mpz (y, zy);
mpz_tdiv_qr (q, r, zx, zy);
scm_remember_upto_here_2 (x, y);
*qp = take_mpz (q);
*rp = take_mpz (r);
}
static SCM
integer_centered_quotient_zz (struct scm_bignum *x, struct scm_bignum *y)
{
mpz_t q, r, min_r, zx, zy;
mpz_init (q);
mpz_init (r);
mpz_init (min_r);
alias_bignum_to_mpz (x, zx);
alias_bignum_to_mpz (y, zy);
/* Note that x might be small enough to fit into a fixnum, so we must
not let it escape into the wild. */
/* min_r will eventually become -abs(y)/2 */
mpz_tdiv_q_2exp (min_r, zy, 1);
/* Arrange for rr to initially be non-positive, because that
simplifies the test to see if it is within the needed bounds. */
if (mpz_sgn (zy) > 0)
{
mpz_cdiv_qr (q, r, zx, zy);
scm_remember_upto_here_2 (x, y);
mpz_neg (min_r, min_r);
if (mpz_cmp (r, min_r) < 0)
mpz_sub_ui (q, q, 1);
}
else
{
mpz_fdiv_qr (q, r, zx, zy);
scm_remember_upto_here_2 (x, y);
if (mpz_cmp (r, min_r) < 0)
mpz_add_ui (q, q, 1);
}
mpz_clear (r);
mpz_clear (min_r);
return take_mpz (q);
}
SCM
scm_integer_centered_quotient_ii (scm_t_inum x, scm_t_inum y)
{
if (y == 0)
scm_num_overflow ("centered-quotient");
scm_t_inum q = x / y;
scm_t_inum r = x % y;
if (x > 0)
{
if (y > 0)
{
if (r >= (y + 1) / 2)
q++;
}
else
{
if (r >= (1 - y) / 2)
q--;
}
}
else
{
if (y > 0)
{
if (r < -y / 2)
q--;
}
else
{
if (r < y / 2)
q++;
}
}
return long_to_scm (q);
}
SCM
scm_integer_centered_quotient_iz (scm_t_inum x, struct scm_bignum *y)
{
return integer_centered_quotient_zz (long_to_bignum (x),
y);
}
SCM
scm_integer_centered_quotient_zi (struct scm_bignum *x, scm_t_inum y)
{
if (y == 0)
scm_num_overflow ("centered-quotient");
else if (y == 1)
return scm_from_bignum (x);
else
{
mpz_t q, zx;
mpz_init (q);
alias_bignum_to_mpz (x, zx);
scm_t_inum r;
/* Arrange for r to initially be non-positive, because that
simplifies the test to see if it is within the needed
bounds. */
if (y > 0)
{
r = - mpz_cdiv_q_ui (q, zx, y);
scm_remember_upto_here_1 (x);
if (r < -y / 2)
mpz_sub_ui (q, q, 1);
}
else
{
r = - mpz_cdiv_q_ui (q, zx, -y);
scm_remember_upto_here_1 (x);
mpz_neg (q, q);
if (r < y / 2)
mpz_add_ui (q, q, 1);
}
return take_mpz (q);
}
}
SCM
scm_integer_centered_quotient_zz (struct scm_bignum *x, struct scm_bignum *y)
{
return integer_centered_quotient_zz (x, y);
}
static SCM
integer_centered_remainder_zz (struct scm_bignum *x, struct scm_bignum *y)
{
mpz_t r, min_r, zx, zy;
mpz_init (r);
mpz_init (min_r);
alias_bignum_to_mpz (x, zx);
alias_bignum_to_mpz (y, zy);
/* Note that x might be small enough to fit into a
fixnum, so we must not let it escape into the wild */
/* min_r will eventually become -abs(y)/2 */
mpz_tdiv_q_2exp (min_r, zy, 1);
/* Arrange for r to initially be non-positive, because that simplifies
the test to see if it is within the needed bounds. */
if (mpz_sgn (zy) > 0)
{
mpz_cdiv_r (r, zx, zy);
mpz_neg (min_r, min_r);
if (mpz_cmp (r, min_r) < 0)
mpz_add (r, r, zy);
}
else
{
mpz_fdiv_r (r, zx, zy);
if (mpz_cmp (r, min_r) < 0)
mpz_sub (r, r, zy);
}
scm_remember_upto_here_2 (x, y);
mpz_clear (min_r);
return take_mpz (r);
}
SCM
scm_integer_centered_remainder_ii (scm_t_inum x, scm_t_inum y)
{
if (y == 0)
scm_num_overflow ("centered-remainder");
scm_t_inum r = x % y;
if (x > 0)
{
if (y > 0)
{
if (r >= (y + 1) / 2)
r -= y;
}
else
{
if (r >= (1 - y) / 2)
r += y;
}
}
else
{
if (y > 0)
{
if (r < -y / 2)
r += y;
}
else
{
if (r < y / 2)
r -= y;
}
}
return SCM_I_MAKINUM (r);
}
SCM
scm_integer_centered_remainder_iz (scm_t_inum x, struct scm_bignum *y)
{
return integer_centered_remainder_zz (long_to_bignum (x),
y);
}
SCM
scm_integer_centered_remainder_zi (struct scm_bignum *x, scm_t_inum y)
{
mpz_t zx;
alias_bignum_to_mpz (x, zx);
if (y == 0)
scm_num_overflow ("centered-remainder");
scm_t_inum r;
/* Arrange for r to initially be non-positive, because that simplifies
the test to see if it is within the needed bounds. */
if (y > 0)
{
r = - mpz_cdiv_ui (zx, y);
if (r < -y / 2)
r += y;
}
else
{
r = - mpz_cdiv_ui (zx, -y);
if (r < y / 2)
r -= y;
}
scm_remember_upto_here_1 (x);
return SCM_I_MAKINUM (r);
}
SCM
scm_integer_centered_remainder_zz (struct scm_bignum *x, struct scm_bignum *y)
{
return integer_centered_remainder_zz (x, y);
}
static void
integer_centered_divide_zz (struct scm_bignum *x, struct scm_bignum *y,
SCM *qp, SCM *rp)
{
mpz_t q, r, min_r, zx, zy;
mpz_init (q);
mpz_init (r);
mpz_init (min_r);
alias_bignum_to_mpz (x, zx);
alias_bignum_to_mpz (y, zy);
/* Note that x might be small enough to fit into a fixnum, so we must
not let it escape into the wild */
/* min_r will eventually become -abs(y/2) */
mpz_tdiv_q_2exp (min_r, zy, 1);
/* Arrange for rr to initially be non-positive, because that
simplifies the test to see if it is within the needed bounds. */
if (mpz_sgn (zy) > 0)
{
mpz_cdiv_qr (q, r, zx, zy);
mpz_neg (min_r, min_r);
if (mpz_cmp (r, min_r) < 0)
{
mpz_sub_ui (q, q, 1);
mpz_add (r, r, zy);
}
}
else
{
mpz_fdiv_qr (q, r, zx, zy);
if (mpz_cmp (r, min_r) < 0)
{
mpz_add_ui (q, q, 1);
mpz_sub (r, r, zy);
}
}
scm_remember_upto_here_2 (x, y);
mpz_clear (min_r);
*qp = take_mpz (q);
*rp = take_mpz (r);
}
void
scm_integer_centered_divide_ii (scm_t_inum x, scm_t_inum y, SCM *qp, SCM *rp)
{
if (y == 0)
scm_num_overflow ("centered-divide");
scm_t_inum q = x / y;
scm_t_inum r = x % y;
if (x > 0)
{
if (y > 0)
{
if (r >= (y + 1) / 2)
{ q++; r -= y; }
}
else
{
if (r >= (1 - y) / 2)
{ q--; r += y; }
}
}
else
{
if (y > 0)
{
if (r < -y / 2)
{ q--; r += y; }
}
else
{
if (r < y / 2)
{ q++; r -= y; }
}
}
*qp = long_to_scm (q);
*rp = SCM_I_MAKINUM (r);
}
void
scm_integer_centered_divide_iz (scm_t_inum x, struct scm_bignum *y, SCM *qp, SCM *rp)
{
integer_centered_divide_zz (long_to_bignum (x), y, qp, rp);
}
void
scm_integer_centered_divide_zi (struct scm_bignum *x, scm_t_inum y, SCM *qp, SCM *rp)
{
if (y == 0)
scm_num_overflow ("centered-divide");
mpz_t q, zx;
mpz_init (q);
alias_bignum_to_mpz (x, zx);
scm_t_inum r;
/* Arrange for r to initially be non-positive, because that
simplifies the test to see if it is within the needed bounds. */
if (y > 0)
{
r = - mpz_cdiv_q_ui (q, zx, y);
if (r < -y / 2)
{
mpz_sub_ui (q, q, 1);
r += y;
}
}
else
{
r = - mpz_cdiv_q_ui (q, zx, -y);
mpz_neg (q, q);
if (r < y / 2)
{
mpz_add_ui (q, q, 1);
r -= y;
}
}
scm_remember_upto_here_1 (x);
*qp = take_mpz (q);
*rp = SCM_I_MAKINUM (r);
}
void
scm_integer_centered_divide_zz (struct scm_bignum *x, struct scm_bignum *y, SCM *qp, SCM *rp)
{
integer_centered_divide_zz (x, y, qp, rp);
}
static SCM
integer_round_quotient_zz (struct scm_bignum *x, struct scm_bignum *y)
{
mpz_t q, r, r2, zx, zy;
int cmp, needs_adjustment;
/* Note that x might be small enough to fit into a
fixnum, so we must not let it escape into the wild */
mpz_init (q);
mpz_init (r);
mpz_init (r2);
alias_bignum_to_mpz (x, zx);
alias_bignum_to_mpz (y, zy);
mpz_fdiv_qr (q, r, zx, zy);
mpz_mul_2exp (r2, r, 1); /* r2 = 2*r */
scm_remember_upto_here_1 (x);
cmp = mpz_cmpabs (r2, zy);
if (mpz_odd_p (q))
needs_adjustment = (cmp >= 0);
else
needs_adjustment = (cmp > 0);
scm_remember_upto_here_1 (y);
if (needs_adjustment)
mpz_add_ui (q, q, 1);
mpz_clear (r);
mpz_clear (r2);
return take_mpz (q);
}
SCM
scm_integer_round_quotient_ii (scm_t_inum x, scm_t_inum y)
{
if (y == 0)
scm_num_overflow ("round-quotient");
scm_t_inum q = x / y;
scm_t_inum r = x % y;
scm_t_inum ay = y;
scm_t_inum r2 = 2 * r;
if (y < 0)
{
ay = -ay;
r2 = -r2;
}
if (q & 1L)
{
if (r2 >= ay)
q++;
else if (r2 <= -ay)
q--;
}
else
{
if (r2 > ay)
q++;
else if (r2 < -ay)
q--;
}
return long_to_scm (q);
}
SCM
scm_integer_round_quotient_iz (scm_t_inum x, struct scm_bignum *y)
{
return integer_round_quotient_zz (long_to_bignum (x), y);
}
SCM
scm_integer_round_quotient_zi (struct scm_bignum *x, scm_t_inum y)
{
if (y == 0)
scm_num_overflow ("round-quotient");
if (y == 1)
return scm_from_bignum (x);
mpz_t q, zx;
mpz_init (q);
alias_bignum_to_mpz (x, zx);
scm_t_inum r;
int needs_adjustment;
if (y > 0)
{
r = mpz_fdiv_q_ui (q, zx, y);
if (mpz_odd_p (q))
needs_adjustment = (2*r >= y);
else
needs_adjustment = (2*r > y);
}
else
{
r = - mpz_cdiv_q_ui (q, zx, -y);
mpz_neg (q, q);
if (mpz_odd_p (q))
needs_adjustment = (2*r <= y);
else
needs_adjustment = (2*r < y);
}
scm_remember_upto_here_1 (x);
if (needs_adjustment)
mpz_add_ui (q, q, 1);
return take_mpz (q);
}
SCM
scm_integer_round_quotient_zz (struct scm_bignum *x, struct scm_bignum *y)
{
mpz_t q, r, zx, zy;
int cmp, needs_adjustment;
mpz_init (q);
mpz_init (r);
alias_bignum_to_mpz (x, zx);
alias_bignum_to_mpz (y, zy);
mpz_fdiv_qr (q, r, zx, zy);
scm_remember_upto_here_1 (x);
mpz_mul_2exp (r, r, 1); /* r = 2*r */
cmp = mpz_cmpabs (r, zy);
mpz_clear (r);
scm_remember_upto_here_1 (y);
if (mpz_odd_p (q))
needs_adjustment = (cmp >= 0);
else
needs_adjustment = (cmp > 0);
if (needs_adjustment)
mpz_add_ui (q, q, 1);
return take_mpz (q);
}
static SCM
integer_round_remainder_zz (struct scm_bignum *x, struct scm_bignum *y)
{
mpz_t q, r, r2, zx, zy;
int cmp, needs_adjustment;
/* Note that x might be small enough to fit into a
fixnum, so we must not let it escape into the wild */
mpz_init (q);
mpz_init (r);
mpz_init (r2);
alias_bignum_to_mpz (x, zx);
alias_bignum_to_mpz (y, zy);
mpz_fdiv_qr (q, r, zx, zy);
scm_remember_upto_here_1 (x);
mpz_mul_2exp (r2, r, 1); /* r2 = 2*r */
cmp = mpz_cmpabs (r2, zy);
if (mpz_odd_p (q))
needs_adjustment = (cmp >= 0);
else
needs_adjustment = (cmp > 0);
if (needs_adjustment)
mpz_sub (r, r, zy);
scm_remember_upto_here_1 (y);
mpz_clear (q);
mpz_clear (r2);
return take_mpz (r);
}
SCM
scm_integer_round_remainder_ii (scm_t_inum x, scm_t_inum y)
{
if (y == 0)
scm_num_overflow ("round-remainder");
scm_t_inum q = x / y;
scm_t_inum r = x % y;
scm_t_inum ay = y;
scm_t_inum r2 = 2 * r;
if (y < 0)
{
ay = -ay;
r2 = -r2;
}
if (q & 1L)
{
if (r2 >= ay)
r -= y;
else if (r2 <= -ay)
r += y;
}
else
{
if (r2 > ay)
r -= y;
else if (r2 < -ay)
r += y;
}
return SCM_I_MAKINUM (r);
}
SCM
scm_integer_round_remainder_iz (scm_t_inum x, struct scm_bignum *y)
{
return integer_round_remainder_zz (long_to_bignum (x), y);
}
SCM
scm_integer_round_remainder_zi (struct scm_bignum *x, scm_t_inum y)
{
if (y == 0)
scm_num_overflow ("round-remainder");
mpz_t q, zx;
scm_t_inum r;
int needs_adjustment;
mpz_init (q);
alias_bignum_to_mpz (x, zx);
if (y > 0)
{
r = mpz_fdiv_q_ui (q, zx, y);
if (mpz_odd_p (q))
needs_adjustment = (2*r >= y);
else
needs_adjustment = (2*r > y);
}
else
{
r = - mpz_cdiv_q_ui (q, zx, -y);
if (mpz_odd_p (q))
needs_adjustment = (2*r <= y);
else
needs_adjustment = (2*r < y);
}
scm_remember_upto_here_1 (x);
mpz_clear (q);
if (needs_adjustment)
r -= y;
return SCM_I_MAKINUM (r);
}
SCM
scm_integer_round_remainder_zz (struct scm_bignum *x, struct scm_bignum *y)
{
return integer_round_remainder_zz (x, y);
}
static void
integer_round_divide_zz (struct scm_bignum *x, struct scm_bignum *y,
SCM *qp, SCM *rp)
{
mpz_t q, r, r2, zx, zy;
int cmp, needs_adjustment;
/* Note that x might be small enough to fit into a fixnum, so we must
not let it escape into the wild */
mpz_init (q);
mpz_init (r);
mpz_init (r2);
alias_bignum_to_mpz (x, zx);
alias_bignum_to_mpz (y, zy);
mpz_fdiv_qr (q, r, zx, zy);
scm_remember_upto_here_1 (x);
mpz_mul_2exp (r2, r, 1); /* r2 = 2*r */
cmp = mpz_cmpabs (r2, zy);
if (mpz_odd_p (q))
needs_adjustment = (cmp >= 0);
else
needs_adjustment = (cmp > 0);
if (needs_adjustment)
{
mpz_add_ui (q, q, 1);
mpz_sub (r, r, zy);
}
scm_remember_upto_here_1 (y);
mpz_clear (r2);
*qp = take_mpz (q);
*rp = take_mpz (r);
}
void
scm_integer_round_divide_ii (scm_t_inum x, scm_t_inum y, SCM *qp, SCM *rp)
{
if (y == 0)
scm_num_overflow ("round-divide");
scm_t_inum q = x / y;
scm_t_inum r = x % y;
scm_t_inum ay = y;
scm_t_inum r2 = 2 * r;
if (y < 0)
{
ay = -ay;
r2 = -r2;
}
if (q & 1L)
{
if (r2 >= ay)
{ q++; r -= y; }
else if (r2 <= -ay)
{ q--; r += y; }
}
else
{
if (r2 > ay)
{ q++; r -= y; }
else if (r2 < -ay)
{ q--; r += y; }
}
*qp = long_to_scm (q);
*rp = SCM_I_MAKINUM (r);
}
void
scm_integer_round_divide_iz (scm_t_inum x, struct scm_bignum *y, SCM *qp, SCM *rp)
{
integer_round_divide_zz (long_to_bignum (x), y, qp, rp);
}
void
scm_integer_round_divide_zi (struct scm_bignum *x, scm_t_inum y, SCM *qp, SCM *rp)
{
if (y == 0)
scm_num_overflow ("round-divide");
mpz_t q, zx;
mpz_init (q);
alias_bignum_to_mpz (x, zx);
scm_t_inum r;
int needs_adjustment;
if (y > 0)
{
r = mpz_fdiv_q_ui (q, zx, y);
if (mpz_odd_p (q))
needs_adjustment = (2*r >= y);
else
needs_adjustment = (2*r > y);
}
else
{
r = - mpz_cdiv_q_ui (q, zx, -y);
mpz_neg (q, q);
if (mpz_odd_p (q))
needs_adjustment = (2*r <= y);
else
needs_adjustment = (2*r < y);
}
scm_remember_upto_here_1 (x);
if (needs_adjustment)
{
mpz_add_ui (q, q, 1);
r -= y;
}
*qp = take_mpz (q);
*rp = SCM_I_MAKINUM (r);
}
void
scm_integer_round_divide_zz (struct scm_bignum *x, struct scm_bignum *y, SCM *qp, SCM *rp)
{
integer_round_divide_zz (x, y, qp, rp);
}
SCM
scm_integer_gcd_ii (scm_t_inum x, scm_t_inum y)
{
scm_t_inum u = x < 0 ? -x : x;
scm_t_inum v = y < 0 ? -y : y;
scm_t_inum result;
if (x == 0)
result = v;
else if (y == 0)
result = u;
else
{
int k = 0;
/* Determine a common factor 2^k */
while (((u | v) & 1) == 0)
{
k++;
u >>= 1;
v >>= 1;
}
/* Now, any factor 2^n can be eliminated */
if ((u & 1) == 0)
while ((u & 1) == 0)
u >>= 1;
else
while ((v & 1) == 0)
v >>= 1;
/* Both u and v are now odd. Subtract the smaller one
from the larger one to produce an even number, remove
more factors of two, and repeat. */
while (u != v)
{
if (u > v)
{
u -= v;
while ((u & 1) == 0)
u >>= 1;
}
else
{
v -= u;
while ((v & 1) == 0)
v >>= 1;
}
}
result = u << k;
}
return ulong_to_scm (result);
}
SCM
scm_integer_gcd_zi (struct scm_bignum *x, scm_t_inum y)
{
scm_t_bits result;
if (y == 0)
return scm_integer_abs_z (x);
if (y < 0)
y = -y;
mpz_t zx;
alias_bignum_to_mpz (x, zx);
result = mpz_gcd_ui (NULL, zx, y);
scm_remember_upto_here_1 (x);
return ulong_to_scm (result);
}
SCM
scm_integer_gcd_zz (struct scm_bignum *x, struct scm_bignum *y)
{
mpz_t result, zx, zy;
mpz_init (result);
alias_bignum_to_mpz (x, zx);
alias_bignum_to_mpz (y, zy);
mpz_gcd (result, zx, zy);
scm_remember_upto_here_2 (x, y);
return take_mpz (result);
}
SCM
scm_integer_lcm_ii (scm_t_inum x, scm_t_inum y)
{
SCM d = scm_integer_gcd_ii (x, y);
if (scm_is_eq (d, SCM_INUM0))
return d;
else
return scm_abs (scm_product (SCM_I_MAKINUM (x),
scm_quotient (SCM_I_MAKINUM (y), d)));
}
SCM
scm_integer_lcm_zi (struct scm_bignum *x, scm_t_inum y)
{
if (y == 0) return SCM_INUM0;
if (y < 0) y = - y;
mpz_t result, zx;
mpz_init (result);
alias_bignum_to_mpz (x, zx);
mpz_lcm_ui (result, zx, y);
scm_remember_upto_here_1 (x);
return take_mpz (result);
}
SCM
scm_integer_lcm_zz (struct scm_bignum *x, struct scm_bignum *y)
{
mpz_t result, zx, zy;
mpz_init (result);
alias_bignum_to_mpz (x, zx);
alias_bignum_to_mpz (y, zy);
mpz_lcm (result, zx, zy);
scm_remember_upto_here_2 (x, y);
/* shouldn't need to normalize b/c lcm of 2 bigs should be big */
return take_mpz (result);
}
/* Emulating 2's complement bignums with sign magnitude arithmetic:
Logand:
X Y Result Method:
(len)
+ + + x (map digit:logand X Y)
+ - + x (map digit:logand X (lognot (+ -1 Y)))
- + + y (map digit:logand (lognot (+ -1 X)) Y)
- - - (+ 1 (map digit:logior (+ -1 X) (+ -1 Y)))
Logior:
X Y Result Method:
+ + + (map digit:logior X Y)
+ - - y (+ 1 (map digit:logand (lognot X) (+ -1 Y)))
- + - x (+ 1 (map digit:logand (+ -1 X) (lognot Y)))
- - - x (+ 1 (map digit:logand (+ -1 X) (+ -1 Y)))
Logxor:
X Y Result Method:
+ + + (map digit:logxor X Y)
+ - - (+ 1 (map digit:logxor X (+ -1 Y)))
- + - (+ 1 (map digit:logxor (+ -1 X) Y))
- - + (map digit:logxor (+ -1 X) (+ -1 Y))
Logtest:
X Y Result
+ + (any digit:logand X Y)
+ - (any digit:logand X (lognot (+ -1 Y)))
- + (any digit:logand (lognot (+ -1 X)) Y)
- - #t
*/
SCM
scm_integer_logand_ii (scm_t_inum x, scm_t_inum y)
{
return SCM_I_MAKINUM (x & y);
}
SCM
scm_integer_logand_zi (struct scm_bignum *x, scm_t_inum y)
{
if (y == 0)
return SCM_INUM0;
if (y > 0)
{
mp_limb_t rd = bignum_limbs (x)[0];
mp_limb_t yd = y;
if (bignum_is_negative (x))
rd = ~rd + 1;
scm_remember_upto_here_1 (x);
rd &= yd;
// Result must be a positive inum.
return SCM_I_MAKINUM (rd);
}
mpz_t result, zx, zy;
mpz_init (result);
alias_bignum_to_mpz (x, zx);
mpz_init_set_si (zy, y);
mpz_and (result, zy, zx);
scm_remember_upto_here_1 (x);
mpz_clear (zy);
return take_mpz (result);
}
SCM
scm_integer_logand_zz (struct scm_bignum *x, struct scm_bignum *y)
{
mpz_t result, zx, zy;
mpz_init (result);
alias_bignum_to_mpz (x, zx);
alias_bignum_to_mpz (y, zy);
mpz_and (result, zx, zy);
scm_remember_upto_here_2 (x, y);
return take_mpz (result);
}
SCM
scm_integer_logior_ii (scm_t_inum x, scm_t_inum y)
{
return SCM_I_MAKINUM (x | y);
}
SCM
scm_integer_logior_zi (struct scm_bignum *x, scm_t_inum y)
{
if (y == 0)
return scm_from_bignum (x);
mpz_t result, zx, zy;
mpz_init (result);
alias_bignum_to_mpz (x, zx);
mpz_init_set_si (zy, y);
mpz_ior (result, zy, zx);
scm_remember_upto_here_1 (x);
mpz_clear (zy);
return take_mpz (result);
}
SCM
scm_integer_logior_zz (struct scm_bignum *x, struct scm_bignum *y)
{
mpz_t result, zx, zy;
mpz_init (result);
alias_bignum_to_mpz (x, zx);
alias_bignum_to_mpz (y, zy);
mpz_ior (result, zy, zx);
scm_remember_upto_here_2 (x, y);
return take_mpz (result);
}
SCM
scm_integer_logxor_ii (scm_t_inum x, scm_t_inum y)
{
return SCM_I_MAKINUM (x ^ y);
}
SCM
scm_integer_logxor_zi (struct scm_bignum *x, scm_t_inum y)
{
mpz_t result, zx, zy;
mpz_init (result);
alias_bignum_to_mpz (x, zx);
mpz_init_set_si (zy, y);
mpz_xor (result, zy, zx);
scm_remember_upto_here_1 (x);
mpz_clear (zy);
return take_mpz (result);
}
SCM
scm_integer_logxor_zz (struct scm_bignum *x, struct scm_bignum *y)
{
mpz_t result, zx, zy;
mpz_init (result);
alias_bignum_to_mpz (x, zx);
alias_bignum_to_mpz (y, zy);
mpz_xor (result, zy, zx);
scm_remember_upto_here_2 (x, y);
return take_mpz (result);
}
int
scm_integer_logtest_ii (scm_t_inum x, scm_t_inum y)
{
return (x & y) ? 1 : 0;
}
int
scm_integer_logtest_zi (struct scm_bignum *x, scm_t_inum y)
{
return scm_is_eq (scm_integer_logand_zi (x, y), SCM_INUM0);
}
int
scm_integer_logtest_zz (struct scm_bignum *x, struct scm_bignum *y)
{
return scm_is_eq (scm_integer_logand_zz (x, y), SCM_INUM0);
}
int
scm_integer_logbit_ui (unsigned long index, scm_t_inum n)
{
if (index < SCM_LONG_BIT)
/* Assume two's complement representation. */
return (n >> index) & 1;
else
return n < 0;
}
int
scm_integer_logbit_uz (unsigned long index, struct scm_bignum *n)
{
mpz_t zn;
alias_bignum_to_mpz (n, zn);
int val = mpz_tstbit (zn, index);
scm_remember_upto_here_1 (n);
return val;
}
SCM
scm_integer_lognot_i (scm_t_inum n)
{
return SCM_I_MAKINUM (~n);
}
SCM
scm_integer_lognot_z (struct scm_bignum *n)
{
mpz_t result, zn;
mpz_init (result);
alias_bignum_to_mpz (n, zn);
mpz_com (result, zn);
scm_remember_upto_here_1 (n);
return take_mpz (result);
}
SCM
scm_integer_expt_ii (scm_t_inum n, scm_t_inum k)
{
ASSERT (k >= 0);
if (k == 0)
return SCM_INUM1;
if (k == 1)
return SCM_I_MAKINUM (n);
if (n == -1)
return scm_is_integer_odd_i (k) ? SCM_I_MAKINUM (-1) : SCM_INUM1;
if (n == 2)
{
if (k < SCM_I_FIXNUM_BIT - 1)
return SCM_I_MAKINUM (1L << k);
if (k < 64)
return scm_integer_from_uint64 (((uint64_t) 1) << k);
size_t nlimbs = k / (sizeof (mp_limb_t)*8) + 1;
size_t high_shift = k & (sizeof (mp_limb_t)*8 - 1);
struct scm_bignum *result = allocate_bignum (nlimbs);
mp_limb_t *rd = bignum_limbs (result);
mpn_zero(rd, nlimbs - 1);
rd[nlimbs - 1] = ((mp_limb_t) 1) << high_shift;
return scm_from_bignum (result);
}
mpz_t res;
mpz_init (res);
mpz_ui_pow_ui (res, inum_magnitude (n), k);
if (n < 0 && (k & 1))
mpz_neg (res, res);
return take_mpz (res);
}
SCM
scm_integer_expt_zi (struct scm_bignum *n, scm_t_inum k)
{
ASSERT (k >= 0);
mpz_t res, zn;
mpz_init (res);
alias_bignum_to_mpz (n, zn);
mpz_pow_ui (res, zn, k);
scm_remember_upto_here_1 (n);
return take_mpz (res);
}
static void
integer_init_mpz (mpz_ptr z, SCM n)
{
if (SCM_I_INUMP (n))
mpz_init_set_si (z, SCM_I_INUM (n));
else
{
ASSERT (SCM_BIGP (n));
mpz_t zn;
alias_bignum_to_mpz (scm_bignum (n), zn);
mpz_init_set (z, zn);
scm_remember_upto_here_1 (n);
}
}
SCM
scm_integer_modulo_expt_nnn (SCM n, SCM k, SCM m)
{
if (scm_is_eq (m, SCM_INUM0))
scm_num_overflow ("modulo-expt");
mpz_t n_tmp, k_tmp, m_tmp;
integer_init_mpz (n_tmp, n);
integer_init_mpz (k_tmp, k);
integer_init_mpz (m_tmp, m);
/* if the exponent K is negative, and we simply call mpz_powm, we
will get a divide-by-zero exception when an inverse 1/n mod m
doesn't exist (or is not unique). Since exceptions are hard to
handle, we'll attempt the inversion "by hand" -- that way, we get
a simple failure code, which is easy to handle. */
if (-1 == mpz_sgn (k_tmp))
{
if (!mpz_invert (n_tmp, n_tmp, m_tmp))
{
mpz_clear (n_tmp);
mpz_clear (k_tmp);
mpz_clear (m_tmp);
scm_num_overflow ("modulo-expt");
}
mpz_neg (k_tmp, k_tmp);
}
mpz_powm (n_tmp, n_tmp, k_tmp, m_tmp);
if (mpz_sgn (m_tmp) < 0 && mpz_sgn (n_tmp) != 0)
mpz_add (n_tmp, n_tmp, m_tmp);
mpz_clear (m_tmp);
mpz_clear (k_tmp);
return take_mpz (n_tmp);
}
/* Efficiently compute (N * 2^COUNT), where N is an exact integer, and
COUNT > 0. */
SCM
scm_integer_lsh_iu (scm_t_inum n, unsigned long count)
{
ASSERT (count > 0);
/* Left shift of count >= SCM_I_FIXNUM_BIT-1 will almost[*] always
overflow a non-zero fixnum. For smaller shifts we check the
bits going into positions above SCM_I_FIXNUM_BIT-1. If they're
all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 - count)".
[*] There's one exception:
(-1) << SCM_I_FIXNUM_BIT-1 == SCM_MOST_NEGATIVE_FIXNUM */
if (n == 0)
return SCM_I_MAKINUM (n);
else if (count < SCM_I_FIXNUM_BIT-1 &&
((scm_t_bits) (SCM_SRS (n, (SCM_I_FIXNUM_BIT-1 - count)) + 1)
<= 1))
return SCM_I_MAKINUM (n < 0 ? -(-n << count) : (n << count));
else
{
mpz_t result;
mpz_init_set_si (result, n);
mpz_mul_2exp (result, result, count);
return take_mpz (result);
}
}
SCM
scm_integer_lsh_zu (struct scm_bignum *n, unsigned long count)
{
ASSERT (count > 0);
mpz_t result, zn;
mpz_init (result);
alias_bignum_to_mpz (n, zn);
mpz_mul_2exp (result, zn, count);
scm_remember_upto_here_1 (n);
return take_mpz (result);
}
/* Efficiently compute floor (N / 2^COUNT), where N is an exact integer
and COUNT > 0. */
SCM
scm_integer_floor_rsh_iu (scm_t_inum n, unsigned long count)
{
ASSERT (count > 0);
if (count >= SCM_I_FIXNUM_BIT)
return (n >= 0 ? SCM_INUM0 : SCM_I_MAKINUM (-1));
else
return SCM_I_MAKINUM (SCM_SRS (n, count));
}
SCM
scm_integer_floor_rsh_zu (struct scm_bignum *n, unsigned long count)
{
ASSERT (count > 0);
mpz_t result, zn;
mpz_init (result);
alias_bignum_to_mpz (n, zn);
mpz_fdiv_q_2exp (result, zn, count);
scm_remember_upto_here_1 (n);
return take_mpz (result);
}
/* Efficiently compute round (N / 2^COUNT), where N is an exact integer
and COUNT > 0. */
SCM
scm_integer_round_rsh_iu (scm_t_inum n, unsigned long count)
{
ASSERT (count > 0);
if (count >= SCM_I_FIXNUM_BIT)
return SCM_INUM0;
else
{
scm_t_inum q = SCM_SRS (n, count);
if (0 == (n & (1L << (count-1))))
return SCM_I_MAKINUM (q); /* round down */
else if (n & ((1L << (count-1)) - 1))
return SCM_I_MAKINUM (q + 1); /* round up */
else
return SCM_I_MAKINUM ((~1L) & (q + 1)); /* round to even */
}
}
SCM
scm_integer_round_rsh_zu (struct scm_bignum *n, unsigned long count)
{
ASSERT (count > 0);
mpz_t q, zn;
mpz_init (q);
alias_bignum_to_mpz (n, zn);
mpz_fdiv_q_2exp (q, zn, count);
if (mpz_tstbit (zn, count-1)
&& (mpz_odd_p (q) || mpz_scan1 (zn, 0) < count-1))
mpz_add_ui (q, q, 1);
scm_remember_upto_here_1 (n);
return take_mpz (q);
}
#define MIN(A, B) ((A) <= (B) ? (A) : (B))
SCM
scm_integer_bit_extract_i (scm_t_inum n, unsigned long start,
unsigned long bits)
{
/* When istart>=SCM_I_FIXNUM_BIT we can just limit the shift to
SCM_I_FIXNUM_BIT-1 to get either 0 or -1 per the sign of "n". */
n = SCM_SRS (n, MIN (start, SCM_I_FIXNUM_BIT-1));
if (n < 0 && bits >= SCM_I_FIXNUM_BIT)
{
/* Since we emulate two's complement encoded numbers, this special
case requires us to produce a result that has more bits than
can be stored in a fixnum. */
mpz_t result;
mpz_init_set_si (result, n);
mpz_fdiv_r_2exp (result, result, bits);
return take_mpz (result);
}
/* mask down to requisite bits */
bits = MIN (bits, SCM_I_FIXNUM_BIT);
return SCM_I_MAKINUM (n & ((1L << bits) - 1));
}
SCM
scm_integer_bit_extract_z (struct scm_bignum *n, unsigned long start, unsigned long bits)
{
mpz_t zn;
alias_bignum_to_mpz (n, zn);
if (bits == 1)
{
int bit = mpz_tstbit (zn, start);
scm_remember_upto_here_1 (n);
return SCM_I_MAKINUM (bit);
}
/* ENHANCE-ME: It'd be nice not to allocate a new bignum when
bits<SCM_I_FIXNUM_BIT. Would want some help from GMP to get
such bits into a ulong. */
mpz_t result;
mpz_init (result);
mpz_fdiv_q_2exp (result, zn, start);
mpz_fdiv_r_2exp (result, result, bits);
scm_remember_upto_here_1 (n);
return take_mpz (result);
}
static const char scm_logtab[] = {
0, 1, 1, 2, 1, 2, 2, 3, 1, 2, 2, 3, 2, 3, 3, 4
};
SCM
scm_integer_logcount_i (scm_t_inum n)
{
unsigned long c = 0;
if (n < 0)
n = -1 - n;
while (n)
{
c += scm_logtab[15 & n];
n >>= 4;
}
return SCM_I_MAKINUM (c);
}
SCM
scm_integer_logcount_z (struct scm_bignum *n)
{
unsigned long count;
mpz_t zn;
alias_bignum_to_mpz (n, zn);
if (mpz_sgn (zn) >= 0)
count = mpz_popcount (zn);
else
{
mpz_t z_negative_one;
mpz_init_set_si (z_negative_one, -1);
count = mpz_hamdist (zn, z_negative_one);
mpz_clear (z_negative_one);
}
scm_remember_upto_here_1 (n);
return scm_from_ulong (count);
}
static const char scm_ilentab[] = {
0, 1, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4
};
SCM
scm_integer_length_i (scm_t_inum n)
{
unsigned long c = 0;
unsigned int l = 4;
if (n < 0)
n = -1 - n;
while (n)
{
c += 4;
l = scm_ilentab [15 & n];
n >>= 4;
}
return SCM_I_MAKINUM (c - 4 + l);
}
SCM
scm_integer_length_z (struct scm_bignum *n)
{
/* mpz_sizeinbase looks at the absolute value of negatives, whereas we
want a ones-complement. If n is ...111100..00 then mpz_sizeinbase is
1 too big, so check for that and adjust. */
mpz_t zn;
alias_bignum_to_mpz (n, zn);
size_t size = mpz_sizeinbase (zn, 2);
/* If negative and no 0 bits above the lowest 1, adjust result. */
if (mpz_sgn (zn) < 0 && mpz_scan0 (zn, mpz_scan1 (zn, 0)) == ULONG_MAX)
size--;
scm_remember_upto_here_1 (n);
return scm_from_size_t (size);
}
SCM
scm_integer_to_string_i (scm_t_inum n, int base)
{
// FIXME: Use mpn_get_str instead.
char num_buf [SCM_INTBUFLEN];
size_t length = scm_iint2str (n, base, num_buf);
return scm_from_latin1_stringn (num_buf, length);
}
SCM
scm_integer_to_string_z (struct scm_bignum *n, int base)
{
mpz_t zn;
alias_bignum_to_mpz (n, zn);
char *str = mpz_get_str (NULL, base, zn);
scm_remember_upto_here_1 (n);
size_t len = strlen (str);
void (*freefunc) (void *, size_t);
mp_get_memory_functions (NULL, NULL, &freefunc);
SCM ret = scm_from_latin1_stringn (str, len);
freefunc (str, len + 1);
return ret;
}
int
scm_is_integer_equal_ir (scm_t_inum x, double y)
{
/* On a 32-bit system an inum fits a double, we can cast the inum
to a double and compare.
But on a 64-bit system an inum is bigger than a double and casting
it to a double (call that dx) will round. Although dxx will not in
general be equal to x, dx will always be an integer and within a
factor of 2 of x, so if dx==y, we know that y is an integer and
fits in scm_t_signed_bits. So we cast y to scm_t_signed_bits and
compare with plain x.
An alternative (for any size system actually) would be to check y
is an integer (with floor) and is in range of an inum (compare
against appropriate powers of 2) then test x==(scm_t_inum)y. It's
just a matter of which casts/comparisons might be fastest or
easiest for the cpu. */
return (double) x == y
&& (DBL_MANT_DIG >= SCM_I_FIXNUM_BIT-1 || x == (scm_t_inum) y);
}
int
scm_is_integer_equal_ic (scm_t_inum x, double real, double imag)
{
return imag == 0.0 && scm_is_integer_equal_ir (x, real);
}
int
scm_is_integer_equal_zz (struct scm_bignum *x, struct scm_bignum *y)
{
mpz_t zx, zy;
alias_bignum_to_mpz (x, zx);
alias_bignum_to_mpz (y, zy);
int cmp = mpz_cmp (zx, zy);
scm_remember_upto_here_2 (x, y);
return 0 == cmp;
}
int
scm_is_integer_equal_zr (struct scm_bignum *x, double y)
{
if (isnan (y))
return 0;
mpz_t zx;
alias_bignum_to_mpz (x, zx);
int cmp = mpz_cmp_d (zx, y);
scm_remember_upto_here_1 (x);
return 0 == cmp;
}
int
scm_is_integer_equal_zc (struct scm_bignum *x, double real, double imag)
{
return imag == 0.0 && scm_is_integer_equal_zr (x, real);
}
int
scm_is_integer_less_than_ir (scm_t_inum x, double y)
{
/* We can safely take the ceiling of y without changing the
result of x<y, given that x is an integer. */
y = ceil (y);
/* In the following comparisons, it's important that the right
hand side always be a power of 2, so that it can be
losslessly converted to a double even on 64-bit
machines. */
if (y >= (double) (SCM_MOST_POSITIVE_FIXNUM+1))
return 1;
else if (!(y > (double) SCM_MOST_NEGATIVE_FIXNUM))
/* The condition above is carefully written to include the
case where y==NaN. */
return 0;
else
/* y is a finite integer that fits in an inum. */
return x < (scm_t_inum) y;
}
int
scm_is_integer_less_than_ri (double x, scm_t_inum y)
{
/* We can safely take the floor of x without changing the
result of x<y, given that y is an integer. */
x = floor (x);
/* In the following comparisons, it's important that the right
hand side always be a power of 2, so that it can be
losslessly converted to a double even on 64-bit
machines. */
if (x < (double) SCM_MOST_NEGATIVE_FIXNUM)
return 1;
else if (!(x < (double) (SCM_MOST_POSITIVE_FIXNUM+1)))
/* The condition above is carefully written to include the
case where x==NaN. */
return 0;
else
/* x is a finite integer that fits in an inum. */
return (scm_t_inum) x < y;
}
int
scm_is_integer_less_than_zz (struct scm_bignum *x, struct scm_bignum *y)
{
mpz_t zx, zy;
alias_bignum_to_mpz (x, zx);
alias_bignum_to_mpz (y, zy);
int cmp = mpz_cmp (zx, zy);
scm_remember_upto_here_2 (x, y);
return cmp < 0;
}
int
scm_is_integer_less_than_zr (struct scm_bignum *x, double y)
{
if (isnan (y))
return 0;
mpz_t zx;
alias_bignum_to_mpz (x, zx);
int cmp = mpz_cmp_d (zx, y);
scm_remember_upto_here_1 (x);
return cmp < 0;
}
int
scm_is_integer_less_than_rz (double x, struct scm_bignum *y)
{
if (isnan (x))
return 0;
mpz_t zy;
alias_bignum_to_mpz (y, zy);
int cmp = mpz_cmp_d (zy, x);
scm_remember_upto_here_1 (y);
return cmp > 0;
}
int
scm_is_integer_positive_z (struct scm_bignum *x)
{
return bignum_is_positive (x);
}
int
scm_is_integer_negative_z (struct scm_bignum *x)
{
return bignum_is_negative (x);
}
#if SCM_ENABLE_MINI_GMP
static double
mpz_get_d_2exp (long *exp, mpz_srcptr z)
{
double signif = mpz_get_d (z);
int iexp;
signif = frexp (signif, &iexp);
*exp = iexp;
return signif;
}
#endif
double
scm_integer_frexp_z (struct scm_bignum *x, long *exp)
{
mpz_t zx;
alias_bignum_to_mpz (x, zx);
size_t bits = mpz_sizeinbase (zx, 2);
ASSERT (bits != 0);
size_t shift = 0;
if (bits > DBL_MANT_DIG)
{
shift = bits - DBL_MANT_DIG;
SCM xx = scm_integer_round_rsh_zu (x, shift);
if (SCM_I_INUMP (xx))
{
int expon;
double signif = frexp (SCM_I_INUM (xx), &expon);
*exp = expon + shift;
return signif;
}
x = scm_bignum (xx);
alias_bignum_to_mpz (x, zx);
}
double significand = mpz_get_d_2exp (exp, zx);
scm_remember_upto_here_1 (x);
*exp += shift;
return significand;
}
double
scm_integer_to_double_z (struct scm_bignum *x)
{
long exponent;
double significand = scm_integer_frexp_z (x, &exponent);
return ldexp (significand, exponent);
}
SCM
scm_integer_from_double (double val)
{
if (!isfinite (val))
scm_out_of_range ("inexact->exact", scm_from_double (val));
if (((double) INT64_MIN) <= val && val <= ((double) INT64_MAX))
return scm_from_int64 (val);
mpz_t result;
mpz_init_set_d (result, val);
return take_mpz (result);
}
SCM
scm_integer_add_ii (scm_t_inum x, scm_t_inum y)
{
return long_to_scm (x + y);
}
static SCM
do_add_1 (int negative, mp_limb_t *xd, size_t xn, mp_limb_t y)
{
size_t rn = xn + 1;
struct scm_bignum *result = allocate_bignum (rn);
mp_limb_t *rd = bignum_limbs (result);
if (mpn_add_1 (rd, xd, xn, y))
rd[xn] = 1;
else
result->u.z.size--;
// No need to normalize as magnitude is increasing and one operand
// already a bignum.
return scm_from_bignum (bignum_negate_if (negative, result));
}
static SCM
do_add (int negative, mp_limb_t *xd, size_t xn, mp_limb_t *yd, size_t yn)
{
size_t rn = xn + 1;
struct scm_bignum *result = allocate_bignum (rn);
mp_limb_t *rd = bignum_limbs (result);
if (mpn_add (rd, xd, xn, yd, yn))
rd[xn] = 1;
else
result->u.z.size--;
// No need to normalize as magnitude is increasing and one operand
// already a bignum.
return scm_from_bignum (bignum_negate_if (negative, result));
}
static SCM
do_sub_1 (int negative, mp_limb_t *xd, size_t xn, mp_limb_t y)
{
size_t rn = xn;
struct scm_bignum *result = allocate_bignum (rn);
mp_limb_t *rd = bignum_limbs (result);
mpn_sub_1 (rd, xd, xn, y);
return normalize_bignum
(bignum_negate_if (negative, (bignum_trim1 (result))));
}
static SCM
do_sub (int negative, mp_limb_t *xd, size_t xn, mp_limb_t *yd, size_t yn)
{
size_t rn = xn;
struct scm_bignum *result = allocate_bignum (rn);
mp_limb_t *rd = bignum_limbs (result);
mpn_sub (rd, xd, xn, yd, yn);
return normalize_bignum
(bignum_negate_if (negative, (bignum_trimn (result))));
}
static int
do_cmp (mp_limb_t *xd, size_t xn, mp_limb_t *yd, size_t yn)
{
if (xn < yn)
return -1;
if (xn > yn)
return 1;
return mpn_cmp (xd, yd, xn);
}
SCM
scm_integer_add_zi (struct scm_bignum *x, scm_t_inum y)
{
if (y == 0)
return scm_from_bignum (x);
size_t xn = bignum_limb_count (x);
if (xn == 0)
return SCM_I_MAKINUM (y);
SCM ret;
if (bignum_is_negative (x) == (y < 0))
// Magnitude increases, sign stays the same.
ret = do_add_1 (y < 0, bignum_limbs (x), xn, inum_magnitude (y));
else
// Magnitude decreases, but assuming x's magnitude is greater than
// y's, not changing sign.
ret = do_sub_1 (bignum_is_negative (x), bignum_limbs (x), xn,
inum_magnitude (y));
scm_remember_upto_here_1 (x);
return ret;
}
SCM
scm_integer_add_zz (struct scm_bignum *x, struct scm_bignum *y)
{
size_t xn = bignum_limb_count (x);
size_t yn = bignum_limb_count (y);
if (xn == 0)
return normalize_bignum (y);
if (yn == 0)
return normalize_bignum (x);
mp_limb_t *xd = bignum_limbs (x);
mp_limb_t *yd = bignum_limbs (y);
SCM ret;
if (bignum_is_negative (x) == bignum_is_negative (y))
// Magnitude increases, sign stays the same.
ret = xn < yn
? do_add (bignum_is_negative (x), yd, yn, xd, xn)
: do_add (bignum_is_negative (x), xd, xn, yd, yn);
else
// Magnitude decreases, changing sign if abs(x) < abs(y).
ret = do_cmp (xd, xn, yd, yn) < 0
? do_sub (!bignum_is_negative (x), yd, yn, xd, xn)
: do_sub (bignum_is_negative (x), xd, xn, yd, yn);
scm_remember_upto_here_2 (x, y);
return ret;
}
SCM
scm_integer_negate_i (scm_t_inum x)
{
return long_to_scm (-x);
}
SCM
scm_integer_negate_z (struct scm_bignum *x)
{
/* Must normalize here because -SCM_MOST_NEGATIVE_FIXNUM is a bignum,
but negating that gives a fixnum. */
return normalize_bignum (negate_bignum (clone_bignum (x)));
}
SCM
scm_integer_sub_ii (scm_t_inum x, scm_t_inum y)
{
// Assumes that -INUM_MIN can fit in a scm_t_inum, even if that
// scm_t_inum is not fixable, and that scm_integer_add_ii can handle
// scm_t_inum inputs outside the fixable range.
return scm_integer_add_ii (x, -y);
}
SCM
scm_integer_sub_iz (scm_t_inum x, struct scm_bignum *y)
{
if (x == 0)
return scm_integer_negate_z (y);
size_t yn = bignum_limb_count (y);
if (yn == 0)
return SCM_I_MAKINUM (x);
SCM ret;
if (bignum_is_negative (y) == (x < 0))
// Magnitude of result smaller than that of y, but assuming y's
// magnitude is greater than x's, keeping y's sign.
ret = do_sub_1 (x > 0, bignum_limbs (y), yn, inum_magnitude (x));
else
// Magnitude increases, same sign as x.
ret = do_add_1 (x < 0, bignum_limbs (y), yn, inum_magnitude (x));
scm_remember_upto_here_1 (y);
return ret;
}
SCM
scm_integer_sub_zi (struct scm_bignum *x, scm_t_inum y)
{
if (y == 0)
return scm_from_bignum (x);
size_t xn = bignum_limb_count (x);
if (xn == 0)
return SCM_I_MAKINUM (y);
SCM ret;
if (bignum_is_negative (x) == (y < 0))
// Magnitude decreases, but assuming x's magnitude is greater than
// y's, not changing sign.
ret = do_sub_1 (y < 0, bignum_limbs (x), xn, inum_magnitude (y));
else
// Magnitude increases, same sign as x.
ret = do_add_1 (bignum_is_negative (x), bignum_limbs (x), xn,
inum_magnitude (y));
scm_remember_upto_here_1 (x);
return ret;
}
SCM
scm_integer_sub_zz (struct scm_bignum *x, struct scm_bignum *y)
{
size_t xn = bignum_limb_count (x);
size_t yn = bignum_limb_count (y);
if (xn == 0)
return scm_integer_negate_z (y);
if (yn == 0)
return scm_from_bignum (x);
mp_limb_t *xd = bignum_limbs (x);
mp_limb_t *yd = bignum_limbs (y);
SCM ret;
if (bignum_is_negative (x) != bignum_is_negative (y))
// Magnitude increases, same sign as x.
ret = xn < yn
? do_add (bignum_is_negative (x), yd, yn, xd, xn)
: do_add (bignum_is_negative (x), xd, xn, yd, yn);
else
// Magnitude decreases, changing sign if abs(x) < abs(y).
ret = do_cmp (xd, xn, yd, yn) < 0
? do_sub (!bignum_is_negative (x), yd, yn, xd, xn)
: do_sub (bignum_is_negative (x), xd, xn, yd, yn);
scm_remember_upto_here_2 (x, y);
return ret;
}
SCM
scm_integer_mul_ii (scm_t_inum x, scm_t_inum y)
{
#if SCM_I_FIXNUM_BIT < 32
int64_t k = x * (int64_t) y;
if (SCM_FIXABLE (k))
return SCM_I_MAKINUM (k);
#endif
mp_limb_t xd[1] = { long_magnitude (x) };
mp_limb_t lo;
int negative = (x < 0) != (y < 0);
mp_limb_t hi = mpn_mul_1 (&lo, xd, 1, long_magnitude (y));
if (!hi)
{
if (negative)
{
if (lo <= long_magnitude (SCM_MOST_NEGATIVE_FIXNUM))
return SCM_I_MAKINUM (negative_long (lo));
}
else if (lo <= SCM_MOST_POSITIVE_FIXNUM)
return SCM_I_MAKINUM (lo);
return scm_from_bignum (make_bignum_1 (negative, lo));
}
return scm_from_bignum (make_bignum_2 (negative, lo, hi));
}
SCM
scm_integer_mul_zi (struct scm_bignum *x, scm_t_inum y)
{
switch (y)
{
case -1:
return scm_integer_negate_z (x);
case 0:
return SCM_INUM0;
case 1:
return scm_from_bignum (x);
default:
{
size_t xn = bignum_limb_count (x);
if (xn == 0)
return SCM_INUM0;
struct scm_bignum *result = allocate_bignum (xn + 1);
mp_limb_t *rd = bignum_limbs (result);
const mp_limb_t *xd = bignum_limbs (x);
mp_limb_t yd = long_magnitude (y);
int negate = bignum_is_negative (x) != (y < 0);
mp_limb_t hi = mpn_mul_1 (rd, xd, xn, yd);
if (hi)
rd[xn] = hi;
else
result->u.z.size--;
scm_remember_upto_here_1 (x);
return normalize_bignum (bignum_negate_if (negate, (result)));
}
}
}
SCM
scm_integer_mul_zz (struct scm_bignum *x, struct scm_bignum *y)
{
size_t xn = bignum_limb_count (x);
size_t yn = bignum_limb_count (y);
if (xn == 0 || yn == 0)
return SCM_INUM0;
struct scm_bignum *result = allocate_bignum (xn + yn);
mp_limb_t *rd = bignum_limbs (result);
const mp_limb_t *xd = bignum_limbs (x);
const mp_limb_t *yd = bignum_limbs (y);
int negate = bignum_is_negative (x) != bignum_is_negative (y);
if (xd == yd)
mpn_sqr (rd, xd, xn);
else if (xn <= yn)
mpn_mul (rd, yd, yn, xd, xn);
else
mpn_mul (rd, xd, xn, yd, yn);
scm_remember_upto_here_2 (x, y);
return normalize_bignum
(bignum_negate_if (negate, (bignum_trim1 (result))));
}
int
scm_is_integer_divisible_ii (scm_t_inum x, scm_t_inum y)
{
ASSERT (y != 0);
return (x % y) == 0;
}
int
scm_is_integer_divisible_zi (struct scm_bignum *x, scm_t_inum y)
{
ASSERT (y != 0);
switch (y)
{
case -1:
case 1:
return 1;
default:
{
scm_t_inum abs_y = y < 0 ? -y : y;
mpz_t zx;
alias_bignum_to_mpz (x, zx);
int divisible = mpz_divisible_ui_p (zx, abs_y);
scm_remember_upto_here_1 (x);
return divisible;
}
}
}
int
scm_is_integer_divisible_zz (struct scm_bignum *x, struct scm_bignum *y)
{
mpz_t zx, zy;
alias_bignum_to_mpz (x, zx);
alias_bignum_to_mpz (y, zy);
int divisible_p = mpz_divisible_p (zx, zy);
scm_remember_upto_here_2 (x, y);
return divisible_p;
}
SCM
scm_integer_exact_quotient_ii (scm_t_inum n, scm_t_inum d)
{
return scm_integer_truncate_quotient_ii (n, d);
}
SCM
scm_integer_exact_quotient_iz (scm_t_inum n, struct scm_bignum *d)
{
// There are only two fixnum numerators that are evenly divided by
// bignum denominators: 0, which is evenly divided 0 times by
// anything, and SCM_MOST_NEGATIVE_FIXNUM, which is evenly divided -1
// time by SCM_MOST_POSITIVE_FIXNUM+1.
if (n == 0)
return SCM_INUM0;
ASSERT (n == SCM_MOST_NEGATIVE_FIXNUM);
ASSERT (bignum_cmp_long (d, SCM_MOST_POSITIVE_FIXNUM + 1) == 0);
return SCM_I_MAKINUM (-1);
}
/* Return the exact integer q such that n = q*d, for exact integers n
and d, where d is known in advance to divide n evenly (with zero
remainder). For large integers, this can be computed more
efficiently than when the remainder is unknown. */
SCM
scm_integer_exact_quotient_zi (struct scm_bignum *n, scm_t_inum d)
{
if (SCM_UNLIKELY (d == 0))
scm_num_overflow ("quotient");
else if (SCM_UNLIKELY (d == 1))
return scm_from_bignum (n);
mpz_t q, zn;
mpz_init (q);
alias_bignum_to_mpz (n, zn);
if (d > 0)
mpz_divexact_ui (q, zn, d);
else
{
mpz_divexact_ui (q, zn, -d);
mpz_neg (q, q);
}
scm_remember_upto_here_1 (n);
return take_mpz (q);
}
SCM
scm_integer_exact_quotient_zz (struct scm_bignum *n, struct scm_bignum *d)
{
mpz_t q, zn, zd;
mpz_init (q);
alias_bignum_to_mpz (n, zn);
alias_bignum_to_mpz (d, zd);
mpz_divexact (q, zn, zd);
scm_remember_upto_here_2 (n, d);
return take_mpz (q);
}
#if SCM_SIZEOF_LONG == 4
SCM
scm_integer_from_int32 (int32_t n)
{
if (SCM_FIXABLE (n))
return SCM_I_MAKINUM (n);
return scm_from_bignum (long_to_bignum (n));
}
SCM
scm_integer_from_uint32 (uint32_t n)
{
if (SCM_POSFIXABLE (n))
return SCM_I_MAKINUM (n);
return scm_from_bignum (ulong_to_bignum (n));
}
int
scm_integer_to_int32_z (struct scm_bignum *z, int32_t *val)
{
return bignum_to_int32 (z, val);
}
int
scm_integer_to_uint32_z (struct scm_bignum *z, uint32_t *val)
{
return bignum_to_uint32 (z, val);
}
#endif
SCM
scm_integer_from_int64 (int64_t n)
{
if (SCM_FIXABLE (n))
return SCM_I_MAKINUM (n);
return scm_from_bignum (make_bignum_from_int64 (n));
}
SCM
scm_integer_from_uint64 (uint64_t n)
{
if (SCM_POSFIXABLE (n))
return SCM_I_MAKINUM (n);
return scm_from_bignum (make_bignum_from_uint64 (n));
}
int
scm_integer_to_int64_z (struct scm_bignum *z, int64_t *val)
{
return bignum_to_int64 (z, val);
}
int
scm_integer_to_uint64_z (struct scm_bignum *z, uint64_t *val)
{
return bignum_to_uint64 (z, val);
}
void
scm_integer_set_mpz_z (struct scm_bignum *z, mpz_t n)
{
mpz_t zn;
alias_bignum_to_mpz (z, zn);
mpz_set (n, zn);
scm_remember_upto_here_1 (z);
}
void
scm_integer_init_set_mpz_z (struct scm_bignum *z, mpz_t n)
{
mpz_init (n);
scm_integer_set_mpz_z (z, n);
}
void
scm_integer_exact_sqrt_i (scm_t_inum k, SCM *s, SCM *r)
{
ASSERT (k >= 0);
if (k == 0)
*s = *r = SCM_INUM0;
else
{
mp_limb_t kk = k, ss, rr;
if (mpn_sqrtrem (&ss, &rr, &kk, 1) == 0)
rr = 0;
*s = SCM_I_MAKINUM (ss);
*r = SCM_I_MAKINUM (rr);
}
}
void
scm_integer_exact_sqrt_z (struct scm_bignum *k, SCM *s, SCM *r)
{
mpz_t zk, zs, zr;
alias_bignum_to_mpz (k, zk);
mpz_init (zs);
mpz_init (zr);
mpz_sqrtrem (zs, zr, zk);
scm_remember_upto_here_1 (k);
*s = take_mpz (zs);
*r = take_mpz (zr);
}
int
scm_is_integer_perfect_square_i (scm_t_inum k)
{
if (k < 0)
return 0;
if (k == 0)
return 1;
mp_limb_t kk = k;
return mpn_perfect_square_p (&kk, 1);
}
int
scm_is_integer_perfect_square_z (struct scm_bignum *k)
{
mpz_t zk;
alias_bignum_to_mpz (k, zk);
int result = mpz_perfect_square_p (zk);
scm_remember_upto_here_1 (k);
return result;
}
SCM
scm_integer_floor_sqrt_i (scm_t_inum k)
{
if (k <= 0)
return SCM_INUM0;
mp_limb_t kk = k, ss;
mpn_sqrtrem (&ss, NULL, &kk, 1);
return SCM_I_MAKINUM (ss);
}
SCM
scm_integer_floor_sqrt_z (struct scm_bignum *k)
{
mpz_t zk, zs;
alias_bignum_to_mpz (k, zk);
mpz_init (zs);
mpz_sqrt (zs, zk);
scm_remember_upto_here_1 (k);
return take_mpz (zs);
}
double
scm_integer_inexact_sqrt_i (scm_t_inum k)
{
if (k < 0)
return -sqrt ((double) -k);
return sqrt ((double) k);
}
double
scm_integer_inexact_sqrt_z (struct scm_bignum *k)
{
long expon;
double signif = scm_integer_frexp_z (k, &expon);
int negative = signif < 0;
if (negative)
signif = -signif;
if (expon & 1)
{
signif *= 2;
expon--;
}
double result = ldexp (sqrt (signif), expon / 2);
return negative ? -result : result;
}
SCM
scm_integer_scan1_i (scm_t_inum n)
{
if (n == 0)
return SCM_I_MAKINUM (-1);
n = n ^ (n-1); /* 1 bits for each low 0 and lowest 1 */
return scm_integer_logcount_i (n >> 1);
}
SCM
scm_integer_scan1_z (struct scm_bignum *n)
{
mpz_t zn;
alias_bignum_to_mpz (n, zn);
unsigned long pos = mpz_scan1 (zn, 0L);
scm_remember_upto_here_1 (n);
return ulong_to_scm (pos);
}
View git blame Copy permalink