Fcalva/chibi-scheme
1
0
Fork 0
mirror of https://github.com/ashinn/chibi-scheme.git synced 2025-05-18 21:29:19 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions
chibi-scheme/bignum.c

1890 lines
58 KiB
C
Raw Blame History

/* bignum.c -- bignum support */
/* Copyright (c) 2009-2015 Alex Shinn. All rights reserved. */
/* BSD-style license: http://synthcode.com/license.txt */
#include "chibi/sexp.h"
#if SEXP_USE_BIGNUMS
#define SEXP_INIT_BIGNUM_SIZE 2
static int digit_value (int c) {
return (((c)<='9') ? ((c) - '0') : ((sexp_toupper(c) - 'A') + 10));
}
static int hex_digit (int n) {
return ((n<=9) ? ('0' + n) : ('A' + n - 10));
}
sexp sexp_make_bignum (sexp ctx, sexp_uint_t len) {
sexp_uint_t size = sexp_sizeof(bignum) + len*sizeof(sexp_uint_t);
sexp res = sexp_alloc_tagged(ctx, size, SEXP_BIGNUM);
if (!sexp_exceptionp(res)) {
sexp_bignum_length(res) = len;
sexp_bignum_sign(res) = 1;
}
return res;
}
sexp sexp_fixnum_to_bignum (sexp ctx, sexp a) {
sexp res = sexp_make_bignum(ctx, 1);
if (!sexp_exceptionp(res)) {
sexp_bignum_data(res)[0] = sexp_unbox_fx_abs(a);
sexp_bignum_sign(res) = sexp_fx_sign(a);
}
return res;
}
sexp sexp_make_integer_from_lsint (sexp ctx, sexp_lsint_t x) {
sexp res;
if (lsint_is_fixnum(x)) {
res = sexp_make_fixnum(lsint_to_sint(x));
} else {
res = sexp_make_bignum(ctx, 1);
if (lsint_lt_0(x)) {
sexp_bignum_sign(res) = -1;
sexp_bignum_data(res)[0] = (sexp_uint_t)-lsint_to_sint(x);
} else {
sexp_bignum_sign(res) = 1;
sexp_bignum_data(res)[0] = (sexp_uint_t)lsint_to_sint(x);
}
}
return res;
}
sexp sexp_make_unsigned_integer_from_luint (sexp ctx, sexp_luint_t x) {
sexp res;
if (luint_is_fixnum(x)) {
res = sexp_make_fixnum(luint_to_uint(x));
} else {
res = sexp_make_bignum(ctx, 1);
sexp_bignum_sign(res) = 1;
sexp_bignum_data(res)[0] = luint_to_uint(x);
}
return res;
}
#if !SEXP_USE_CUSTOM_LONG_LONGS
sexp sexp_make_integer (sexp ctx, sexp_lsint_t x) {
return sexp_make_integer_from_lsint(ctx, x);
}
sexp sexp_make_unsigned_integer (sexp ctx, sexp_luint_t x) {
return sexp_make_unsigned_integer_from_luint(ctx, x);
}
#endif /* !SEXP_USE_CUSTOM_LONG_LONGS */
#define double_trunc_10s_digit(f) (trunc((f)/10.0)*10.0)
#define double_10s_digit(f) ((f)-double_trunc_10s_digit(f))
#define double_16s_digit(f) fmod(f,16.0)
sexp sexp_double_to_bignum (sexp ctx, double f) {
int sign;
sexp_gc_var3(res, scale, tmp);
sexp_gc_preserve3(ctx, res, scale, tmp);
res = sexp_fixnum_to_bignum(ctx, SEXP_ZERO);
scale = sexp_fixnum_to_bignum(ctx, SEXP_ONE);
sign = (f < 0 ? -1 : 1);
for (f=fabs(f); f >= 1.0; f=trunc(f/16)) {
tmp = sexp_bignum_fxmul(ctx, NULL, scale, (sexp_uint_t)double_16s_digit(f), 0);
res = sexp_bignum_add(ctx, res, res, tmp);
scale = sexp_bignum_fxmul(ctx, NULL, scale, 16, 0);
}
sexp_bignum_sign(res) = sign;
sexp_gc_release3(ctx);
return res;
}
sexp sexp_copy_bignum (sexp ctx, sexp dst, sexp a, sexp_uint_t len0) {
sexp_uint_t len = (len0 > 0) ? len0 : sexp_bignum_length(a), size;
size = sexp_sizeof(bignum) + len*sizeof(sexp_uint_t);
if (! dst || sexp_bignum_length(dst) < len) {
dst = sexp_alloc_tagged(ctx, size, SEXP_BIGNUM);
sexp_bignum_length(dst) = len;
}
if (sexp_bignum_length(a) < len)
len = sexp_bignum_length(a);
sexp_bignum_sign(dst) = sexp_bignum_sign(a);
memset(sexp_bignum_data(dst), 0,
sexp_bignum_length(dst)*sizeof(sexp_uint_t));
memmove(sexp_bignum_data(dst), sexp_bignum_data(a),
len*sizeof(sexp_uint_t));
return dst;
}
int sexp_bignum_zerop (sexp a) {
int i;
sexp_uint_t *data = sexp_bignum_data(a);
for (i=sexp_bignum_length(a)-1; i>=0; i--)
if (data[i])
return 0;
return 1;
}
sexp_uint_t sexp_bignum_hi (sexp a) {
sexp_uint_t i=sexp_bignum_length(a)-1;
while ((i>0) && ! sexp_bignum_data(a)[i])
i--;
return i+1;
}
sexp_sint_t sexp_bignum_compare_abs (sexp a, sexp b) {
int ai=sexp_bignum_hi(a), bi=sexp_bignum_hi(b);
sexp_uint_t *adata=sexp_bignum_data(a), *bdata=sexp_bignum_data(b);
if (ai != bi)
return ai - bi;
for (--ai; ai >= 0; ai--) {
if (adata[ai] > bdata[ai])
return 1;
else if (adata[ai] < bdata[ai])
return -1;
}
return 0;
}
sexp_sint_t sexp_bignum_compare (sexp a, sexp b) {
if (sexp_bignum_sign(a) != sexp_bignum_sign(b))
return sexp_bignum_sign(a);
return sexp_bignum_compare_abs(a, b);
}
sexp sexp_bignum_normalize (sexp a) {
sexp_uint_t *data;
if ((! sexp_bignump(a)) || (sexp_bignum_hi(a)>1))
return a;
data = sexp_bignum_data(a);
if ((data[0] > SEXP_MAX_FIXNUM)
&& ! ((sexp_bignum_sign(a) == -1) && (data[0] == SEXP_MAX_FIXNUM+1)))
return a;
return sexp_make_fixnum((sexp_sint_t)data[0] * sexp_bignum_sign(a));
}
double sexp_bignum_to_double (sexp a) {
double res = 0;
sexp_sint_t i;
sexp_uint_t *data=sexp_bignum_data(a);
for (i=sexp_bignum_hi(a)-1; i>=0; i--)
res = res * ((double)SEXP_UINT_T_MAX+1) + data[i];
return res * sexp_bignum_sign(a);
}
sexp sexp_bignum_fxadd (sexp ctx, sexp a, sexp_uint_t b) {
sexp_uint_t len=sexp_bignum_hi(a), *data=sexp_bignum_data(a),
carry=b, i=0, n;
do { n = data[i];
data[i] += carry;
carry = (n > (SEXP_UINT_T_MAX - carry));
} while (++i<len && carry);
if (carry) {
a = sexp_copy_bignum(ctx, NULL, a, len+1);
sexp_bignum_data(a)[len] = 1;
}
return a;
}
sexp sexp_bignum_fxsub (sexp ctx, sexp a, sexp_uint_t b) {
sexp_uint_t len=sexp_bignum_hi(a), *data=sexp_bignum_data(a), borrow, i=0, n;
if ((len == 1) && (b > data[0])) {
data[0] = b - data[0];
sexp_bignum_sign(a) = -sexp_bignum_sign(a);
} else {
for (borrow=b; borrow; i++) {
n = data[i];
data[i] -= borrow;
borrow = (n < borrow);
}
}
return a;
}
sexp sexp_bignum_fxmul (sexp ctx, sexp d, sexp a, sexp_uint_t b, int offset) {
sexp_uint_t len=sexp_bignum_length(a), *data, *adata=sexp_bignum_data(a),
carry=0, i;
sexp_luint_t n;
sexp_gc_var1(tmp);
sexp_gc_preserve1(ctx, tmp);
if ((! d) || (sexp_bignum_length(d) < len+offset))
d = sexp_make_bignum(ctx, len+offset);
tmp = d;
data = sexp_bignum_data(d);
for (i=0; i<len; i++) {
n = luint_add(luint_mul_uint(luint_from_uint(adata[i]), b), luint_from_uint(carry));
data[i+offset] = luint_to_uint(n);
carry = luint_to_uint(luint_shr(n, (sizeof(sexp_uint_t)*8)));
}
if (carry) {
if (sexp_bignum_length(d) <= len+offset)
d = sexp_copy_bignum(ctx, NULL, d, len+offset+1);
sexp_bignum_data(d)[len+offset] = carry;
}
sexp_gc_release1(ctx);
return d;
}
sexp_uint_t sexp_bignum_fxdiv (sexp ctx, sexp a, sexp_uint_t b, int offset) {
sexp_uint_t len=sexp_bignum_hi(a), *data=sexp_bignum_data(a), q, r=0;
int i;
sexp_luint_t n = luint_from_uint(0);
for (i=len-1; i>=offset; i--) {
n = luint_add(luint_shl(n, sizeof(sexp_uint_t)*8), luint_from_uint(data[i]));
q = luint_to_uint(luint_div_uint(n, b));
r = luint_to_uint(luint_sub(n, luint_mul_uint(luint_from_uint(q), b)));
data[i] = q;
n = luint_from_uint(r);
}
return r;
}
sexp sexp_bignum_fxrem (sexp ctx, sexp a, sexp_sint_t b) {
sexp_uint_t len=sexp_bignum_hi(a), *data=sexp_bignum_data(a), q, b0;
int i;
sexp_luint_t n = luint_from_uint(0);
if (b > 0) {
q = b - 1;
if ((b & q) == 0)
return sexp_make_fixnum(sexp_bignum_sign(a) * (data[0] & q));
}
b0 = (b >= 0) ? b : -b;
if (b0 == 0) {
return sexp_xtype_exception(ctx, NULL, "divide by zero", a);
}
for (i=len-1; i>=0; i--) {
n = luint_add(luint_shl(n, sizeof(sexp_uint_t)*8), luint_from_uint(data[i]));
q = luint_to_uint(luint_div_uint(n, b0));
n = luint_sub(n, luint_mul_uint(luint_from_uint(q), b0));
}
return sexp_make_fixnum(sexp_bignum_sign(a) * (sexp_sint_t)luint_to_uint(n));
}
sexp sexp_read_bignum (sexp ctx, sexp in, sexp_uint_t init,
signed char sign, sexp_uint_t base) {
int c, digit;
sexp_gc_var3(res, tmp, imag);
sexp_gc_preserve3(ctx, res, tmp, imag);
res = sexp_make_bignum(ctx, SEXP_INIT_BIGNUM_SIZE);
sexp_bignum_sign(res) = sign;
sexp_bignum_data(res)[0] = init;
for (c=sexp_read_char(ctx, in); sexp_isxdigit(c); c=sexp_read_char(ctx, in)) {
digit = digit_value(c);
if ((digit < 0) || (digit >= (int)base))
break;
res = sexp_bignum_fxmul(ctx, res, res, base, 0);
res = sexp_bignum_fxadd(ctx, res, digit);
}
if (c=='.' || c=='e' || c=='E') {
if (base != 10) {
res = sexp_read_error(ctx, "found non-base 10 float", SEXP_NULL, in);
} else if (c=='.') {
res = sexp_read_float_tail(ctx, in, sexp_bignum_to_double(res), (sign==-1));
} else {
tmp = sexp_read_number(ctx, in, base, 0);
#if SEXP_USE_COMPLEX
if (sexp_complexp(tmp)) {
imag = sexp_complex_imag(tmp);
tmp = sexp_complex_real(tmp);
} else {
imag = SEXP_ZERO;
}
#endif
if (sexp_exceptionp(tmp)) {
res = tmp;
} else if (sexp_fixnump(tmp) && labs(sexp_unbox_fixnum(tmp)) < 100*1024*1024) {
tmp = sexp_expt(ctx, SEXP_TEN, tmp);
res = sexp_mul(ctx, res, tmp);
} else {
tmp = sexp_exact_to_inexact(ctx, NULL, 2, tmp);
tmp = sexp_expt(ctx, SEXP_TEN, tmp);
res = sexp_mul(ctx, res, tmp);
}
#if SEXP_USE_COMPLEX
if (imag != SEXP_ZERO && !sexp_exceptionp(res))
res = sexp_make_complex(ctx, res, imag);
#endif
}
#if SEXP_USE_RATIOS
} else if (c=='/') {
res = sexp_bignum_normalize(res);
res = sexp_make_ratio(ctx, res, SEXP_ONE);
sexp_ratio_denominator(res) = sexp_read_number(ctx, in, 10, 0);
res = sexp_ratio_normalize(ctx, res, in);
#endif
#if SEXP_USE_COMPLEX
} else if (c=='i' || c=='i' || c=='+' || c=='-') {
sexp_push_char(ctx, c, in);
res = sexp_bignum_normalize(res);
res = sexp_read_complex_tail(ctx, in, res);
#endif
} else if ((c!=EOF) && ! sexp_is_separator(c)) {
res = sexp_read_error(ctx, "invalid numeric syntax",
sexp_make_character(c), in);
} else {
sexp_push_char(ctx, c, in);
}
sexp_gc_release3(ctx);
return sexp_bignum_normalize(res);
}
static int log2i(int v) {
int i;
for (i = 0; i < sizeof(v)*8; i++)
if ((1<<(i+1)) > v)
break;
return i;
}
sexp sexp_write_bignum (sexp ctx, sexp a, sexp out, sexp_uint_t base) {
int i, str_len, lg_base = log2i(base);
char *data;
sexp_gc_var2(b, str);
sexp_gc_preserve2(ctx, b, str);
b = sexp_copy_bignum(ctx, NULL, a, 0);
sexp_bignum_sign(b) = 1;
if (lg_base < 1) {
return sexp_xtype_exception(ctx, NULL, "number base too small", a);
}
i = str_len = (sexp_bignum_length(b)*sizeof(sexp_uint_t)*8 + lg_base - 1)
/ lg_base + 1;
str = sexp_make_string(ctx, sexp_make_fixnum(str_len),
sexp_make_character(' '));
data = sexp_string_data(str);
while (! sexp_bignum_zerop(b))
data[--i] = hex_digit(sexp_bignum_fxdiv(ctx, b, base, 0));
if (i == str_len)
data[--i] = '0';
else if (sexp_bignum_sign(a) == -1)
data[--i] = '-';
sexp_write_string(ctx, data + i, out);
sexp_gc_release2(ctx);
return SEXP_VOID;
}
/****************** bignum arithmetic *************************/
sexp sexp_bignum_add_fixnum (sexp ctx, sexp a, sexp b) {
sexp_gc_var1(c);
sexp_gc_preserve1(ctx, c);
c = sexp_copy_bignum(ctx, NULL, a, 0);
if (sexp_bignum_sign(c) == sexp_fx_sign(b))
c = sexp_bignum_fxadd(ctx, c, sexp_unbox_fx_abs(b));
else
c = sexp_bignum_fxsub(ctx, c, sexp_unbox_fx_abs(b));
sexp_gc_release1(ctx);
return c;
}
sexp sexp_bignum_sub_digits (sexp ctx, sexp dst, sexp a, sexp b) {
sexp_uint_t alen=sexp_bignum_hi(a), blen=sexp_bignum_hi(b),
borrow=0, i, *adata, *bdata, *cdata;
sexp_gc_var1(c);
if ((alen < blen) || ((alen == blen) && (sexp_bignum_compare_abs(a, b) < 0)))
return sexp_bignum_sub_digits(ctx, dst, b, a);
sexp_gc_preserve1(ctx, c);
c = ((dst && sexp_bignum_hi(dst) >= alen)
? dst : sexp_copy_bignum(ctx, NULL, a, 0));
adata = sexp_bignum_data(a);
bdata = sexp_bignum_data(b);
cdata = sexp_bignum_data(c);
for (i=0; i<blen; i++) {
if (adata[i] > bdata[i] || (adata[i] == bdata[i] && !borrow)) {
cdata[i] = adata[i] - bdata[i] - borrow;
borrow = 0;
} else {
cdata[i] = (SEXP_UINT_T_MAX - bdata[i]);
cdata[i] += 1;
cdata[i] -= borrow;
cdata[i] += adata[i];
borrow = 1;
}
}
for ( ; borrow && (i<alen); i++) {
borrow = (cdata[i] == 0 ? 1 : 0);
cdata[i]--;
}
sexp_gc_release1(ctx);
return c;
}
sexp sexp_bignum_add_digits (sexp ctx, sexp dst, sexp a, sexp b) {
sexp_uint_t alen=sexp_bignum_hi(a), blen=sexp_bignum_hi(b),
carry=0, i, old_a, p_sum, *adata, *bdata, *cdata;
sexp_gc_var1(c);
if (alen < blen) return sexp_bignum_add_digits(ctx, dst, b, a);
sexp_gc_preserve1(ctx, c);
c = ((dst && sexp_bignum_hi(dst) >= alen)
? dst : sexp_copy_bignum(ctx, NULL, a, 0));
adata = sexp_bignum_data(a);
bdata = sexp_bignum_data(b);
cdata = sexp_bignum_data(c);
for (i=0; i<blen; i++) {
old_a = adata[i]; /* adata may alias cdata */
p_sum = adata[i] + bdata[i];
cdata[i] = p_sum + carry;
carry = (old_a > (SEXP_UINT_T_MAX - bdata[i]) ? 1 : 0)
+ (p_sum > (SEXP_UINT_T_MAX - carry) ? 1 : 0);
}
for ( ; carry && (i<alen); i++) {
carry = (cdata[i] == SEXP_UINT_T_MAX ? 1 : 0);
cdata[i]++;
}
if (carry) {
c = sexp_copy_bignum(ctx, NULL, c, alen+1);
sexp_bignum_data(c)[alen] = 1;
}
sexp_gc_release1(ctx);
return c;
}
sexp sexp_bignum_add (sexp ctx, sexp dst, sexp a, sexp b) {
sexp res;
if (sexp_bignum_sign(a) == sexp_bignum_sign(b)) {
res = sexp_bignum_add_digits(ctx, dst, a, b);
sexp_bignum_sign(res) = sexp_bignum_sign(a);
} else {
res = sexp_bignum_sub_digits(ctx, dst, a, b);
sexp_bignum_sign(res)
= sexp_bignum_sign(sexp_bignum_compare_abs(a, b) >= 0 ? a : b);
}
return res;
}
sexp sexp_bignum_sub (sexp ctx, sexp dst, sexp a, sexp b) {
sexp res;
int sign;
if (sexp_bignum_sign(a) == sexp_bignum_sign(b)) {
sign = (sexp_bignum_compare_abs(a, b) >= 0 ? sexp_bignum_sign(a)
: -sexp_bignum_sign(a));
res = sexp_bignum_sub_digits(ctx, dst, a, b);
} else {
sign = sexp_bignum_sign(a);
res = sexp_bignum_add_digits(ctx, dst, a, b);
}
sexp_bignum_sign(res) = sign;
return res;
}
static void sexp_bignum_split (sexp ctx, sexp a, sexp_uint_t k, sexp* lo, sexp* hi) {
sexp_uint_t alen=sexp_bignum_hi(a), i, *adata=sexp_bignum_data(a),
*lodata, *hidata;
*lo = sexp_make_bignum(ctx, k); /* must be gc protected by caller */
*hi = sexp_make_bignum(ctx, alen-k+1);
lodata = sexp_bignum_data(*lo);
hidata = sexp_bignum_data(*hi);
for (i=0; i<k; i++) /* split into a[0..k-1], a[k..] */
lodata[i] = adata[i];
for (i=k; i<alen; i++)
hidata[i-k] = adata[i];
}
static sexp sexp_bignum_shift (sexp ctx, sexp a, sexp_uint_t k) {
sexp res;
sexp_uint_t alen = sexp_bignum_hi(a), i;
res = sexp_make_bignum(ctx, alen + k + 1);
for (i=0; i<alen; i++)
sexp_bignum_data(res)[i+k] = sexp_bignum_data(a)[i];
return res;
}
sexp sexp_bignum_mul (sexp ctx, sexp dst, sexp a, sexp b) {
sexp_uint_t alen=sexp_bignum_hi(a), blen=sexp_bignum_hi(b), k,
*bdata=sexp_bignum_data(b);
sexp_gc_var7(a0, a1, b0, b1, z0, z1, z2);
if (alen < blen) return sexp_bignum_mul(ctx, dst, b, a);
if (blen == 1) {
z1 = sexp_bignum_fxmul(ctx, dst, a, bdata[0], 0);
} else {
/* karatsuba: */
/* ab = (a1B^k + a0) * (b1B^k + b0) */
/* = z2B^2k + z1B^k + z0 */
/* where: */
/* z2 = a1b1 */
/* z1 = a1b0 + a0b1 */
/* z0 = a0b0 */
/* then optimize further: */
/* z1 = (a1 + a0)(b1 + b0) - z2 - z0 */
sexp_gc_preserve7(ctx, a0, a1, b0, b1, z0, z1, z2);
k = blen / 2;
sexp_bignum_split(ctx, a, k, &a0, &a1);
sexp_bignum_split(ctx, b, k, &b0, &b1);
z0 = sexp_bignum_add(ctx, NULL, a1, a0); /* temp */
z1 = sexp_bignum_add(ctx, NULL, b1, b0);
z1 = sexp_bignum_mul(ctx, NULL, z1, z0); /* 1 */
z0 = sexp_bignum_mul(ctx, NULL, a0, b0); /* 2 */
z2 = sexp_bignum_mul(ctx, NULL, a1, b1); /* 3 */
z1 = sexp_bignum_sub(ctx, NULL, z1, z0);
z1 = sexp_bignum_sub(ctx, NULL, z1, z2);
z2 = sexp_bignum_shift(ctx, z2, 2*k);
z1 = sexp_bignum_shift(ctx, z1, k);
z1 = sexp_bignum_add(ctx, z1, z1, z0);
z1 = sexp_bignum_add(ctx, z1, z1, z2);
sexp_gc_release7(ctx);
}
sexp_bignum_sign(z1) = sexp_bignum_sign(a) * sexp_bignum_sign(b);
return z1;
}
sexp sexp_bignum_quot_rem (sexp ctx, sexp *rem, sexp a, sexp b) {
sexp_luint_t d, dn, dd;
sexp_uint_t dlo, dhi;
sexp_sint_t alen, blen=sexp_bignum_hi(b), sign=1, off=0;
sexp_gc_var5(q, x, y, a1, b1);
if (blen == 1 && sexp_bignum_data(b)[0] == 0)
return sexp_xtype_exception(ctx, NULL, "divide by zero", a);
sexp_gc_preserve5(ctx, q, x, y, a1, b1);
a1 = sexp_copy_bignum(ctx, NULL, a, 0);
sexp_bignum_sign(a1) = 1;
/* fast path for single bigit divisor */
if (blen == 1) {
b1 = sexp_make_bignum(ctx, 1);
sexp_bignum_data(b1)[0] = sexp_bignum_fxdiv(ctx, a1, sexp_bignum_data(b)[0], 0);
*rem = sexp_bignum_normalize(b1);
if (sexp_bignum_sign(a) * sexp_bignum_sign(b) < 0) {
sexp_negate_exact(a1);
}
if (sexp_bignum_sign(a) < 0) {
sexp_negate_exact(*rem);
}
sexp_gc_release5(ctx);
return a1;
}
/* general case */
b1 = sexp_copy_bignum(ctx, NULL, b, 0);
sexp_bignum_sign(b1) = 1;
q = SEXP_ZERO;
x = sexp_make_bignum(ctx, sexp_bignum_length(a));
while (sexp_bignum_compare_abs(a1, b1) >= 0) { /* a1, b1 at least 2 bigits */
/* guess divisor x */
alen = sexp_bignum_hi(a1);
sexp_bignum_data(x)[off] = 0;
if (off > 0) sexp_bignum_data(x)[off-1] = 0;
off = alen - blen + 1;
dn = luint_add_uint(luint_shl(luint_from_uint(sexp_bignum_data(a1)[alen-1])
, (sizeof(sexp_uint_t)*8))
, sexp_bignum_data(a1)[alen-2]);
dd = luint_add_uint(luint_shl(luint_from_uint(sexp_bignum_data(b1)[blen-1])
, (sizeof(sexp_uint_t)*8))
, sexp_bignum_data(b1)[blen-2]);
if (alen > 2 && blen > 2 &&
luint_lt(luint_from_uint(sexp_bignum_data(a1)[alen-1]), luint_shl(luint_from_uint(1), (sizeof(sexp_uint_t)*4))) &&
luint_lt(luint_from_uint(sexp_bignum_data(b1)[blen-1]), luint_shl(luint_from_uint(1), (sizeof(sexp_uint_t)*4)))) {
dn = luint_add_uint(luint_shl(dn, (sizeof(sexp_uint_t)*4))
, (sexp_bignum_data(a1)[alen-3] >> (sizeof(sexp_uint_t)*4)));
dd = luint_add_uint(luint_shl(dd, (sizeof(sexp_uint_t)*4))
, (sexp_bignum_data(b1)[blen-3] >> (sizeof(sexp_uint_t)*4)));
}
d = luint_div(dn, dd);
if (luint_eq(d, luint_from_uint(0))) {
dn = luint_add_uint(luint_shl(luint_from_uint(sexp_bignum_data(a1)[alen-1])
, (sizeof(sexp_uint_t)*8))
, sexp_bignum_data(a1)[alen-2]);
dd = luint_from_uint(sexp_bignum_data(b1)[blen-1]);
if (luint_lt(luint_from_uint(sexp_bignum_data(a1)[alen-1]), (luint_shl(luint_from_uint(1), (sizeof(sexp_uint_t)*4)))) &&
luint_lt(luint_from_uint(sexp_bignum_data(b1)[blen-1]), (luint_shl(luint_from_uint(1), (sizeof(sexp_uint_t)*4))))) {
dn = luint_add_uint(luint_shl(dn, (sizeof(sexp_uint_t)*4))
, (sexp_bignum_data(a1)[alen-3] >> (sizeof(sexp_uint_t)*4)));
dd = luint_add_uint(luint_shl(dd, (sizeof(sexp_uint_t)*4))
, (sexp_bignum_data(b1)[blen-2] >> (sizeof(sexp_uint_t)*4)));
}
d = luint_div(dn, dd);
off--;
}
dhi = luint_to_uint(luint_shr(d, (sizeof(sexp_uint_t)*8)));
dlo = luint_to_uint(luint_and(d, luint_sub(luint_shl(luint_from_uint(1), (sizeof(sexp_uint_t)*8)), luint_from_uint(1))));
sexp_bignum_data(x)[off] = dhi;
if (off > 0) sexp_bignum_data(x)[off-1] = dlo;
/* update quotient q and remainder a1 estimates */
y = sexp_bignum_mul(ctx, NULL, b1, x);
if (sign < 0) {
a1 = sexp_bignum_add(ctx, NULL, a1, y);
q = sexp_sub(ctx, q, x);
} else {
a1 = sexp_bignum_sub(ctx, NULL, a1, y);
q = sexp_add(ctx, q, x);
}
/* flip the sign if we overshot in our estimate */
if (sexp_bignum_sign(a1) != sign) {
sexp_bignum_sign(a1) = (char)(-sign);
sign *= -1;
}
}
/* adjust signs */
a1 = sexp_bignum_normalize(a1);
if (sign < 0 && a1 != SEXP_ZERO) {
q = sexp_sub(ctx, q, SEXP_ONE);
a1 = sexp_add(ctx, a1, b1);
}
*rem = a1;
if (sexp_bignum_sign(a) * sexp_bignum_sign(b) < 0) {
sexp_negate_exact(q);
}
if (sexp_bignum_sign(a) < 0) {
sexp_negate_exact(*rem);
}
sexp_gc_release5(ctx);
return q;
}
sexp sexp_bignum_quotient (sexp ctx, sexp a, sexp b) {
sexp res;
sexp_gc_var1(rem);
sexp_gc_preserve1(ctx, rem);
res = sexp_bignum_quot_rem(ctx, &rem, a, b);
sexp_gc_release1(ctx);
return res;
}
sexp sexp_bignum_remainder (sexp ctx, sexp a, sexp b) {
sexp_gc_var1(rem);
sexp_gc_preserve1(ctx, rem);
sexp_bignum_quot_rem(ctx, &rem, a, b); /* discard quotient */
sexp_gc_release1(ctx);
return rem;
}
sexp sexp_bignum_expt (sexp ctx, sexp a, sexp b) {
sexp_sint_t e = sexp_unbox_fixnum(b);
sexp_sint_t abs_e;
if (e < 0)
abs_e = -e;
else
abs_e = e;
sexp_gc_var2(res, acc);
sexp_gc_preserve2(ctx, res, acc);
res = sexp_fixnum_to_bignum(ctx, SEXP_ONE);
acc = sexp_copy_bignum(ctx, NULL, a, 0);
for (; abs_e; abs_e>>=1, acc=sexp_bignum_mul(ctx, NULL, acc, acc))
if (abs_e & 1)
res = sexp_bignum_mul(ctx, NULL, res, acc);
if (e < 0)
res = sexp_div(ctx, sexp_fixnum_to_bignum(ctx, SEXP_ONE), res);
sexp_gc_release2(ctx);
return sexp_bignum_normalize(res);
}
#if SEXP_USE_MATH
/*
* a = x * 2^2n, with 0.1 <= x < 1.0 (base 2) => sqrt(a) ~ 2^n
*/
sexp sexp_bignum_sqrt_estimate (sexp ctx, sexp a) {
sexp_uint_t alen=sexp_bignum_hi(a), adata_hi;
int nbits, i;
sexp_gc_var1(res);
adata_hi = sexp_bignum_data(a)[alen - 1];
for (i = sizeof(sexp_uint_t)*8-1; i > 0; i--)
if (adata_hi & (1ul << i))
break;
nbits = sizeof(sexp_uint_t) * 8 * (alen - 1) + i + 1;
sexp_gc_preserve1(ctx, res);
res = sexp_fixnum_to_bignum(ctx, SEXP_TWO);
res = sexp_bignum_expt(ctx, res, sexp_make_fixnum(nbits / 2));
sexp_gc_release1(ctx);
return res;
}
#define SEXP_MAX_ACCURATE_FLONUM_SQRT 1073741824.0 /* 2^30 */
/* Babylonian method */
sexp sexp_bignum_sqrt (sexp ctx, sexp a, sexp* rem_out) {
sexp_gc_var4(res, rem, tmp, tmpa);
if (! sexp_bignump(a)) return sexp_type_exception(ctx, NULL, SEXP_BIGNUM, a);
sexp_gc_preserve4(ctx, res, rem, tmp, tmpa);
/* initial estimate via flonum, ignoring signs */
if (sexp_exact_negativep(a)) {
tmpa = sexp_copy_bignum(ctx, NULL, a, 0);
a = tmpa;
sexp_negate(a);
}
res = sexp_make_flonum(ctx, sexp_bignum_to_double(a));
res = sexp_inexact_sqrt(ctx, NULL, 1, res);
if (sexp_flonump(res) &&
sexp_flonum_value(res) > SEXP_MAX_ACCURATE_FLONUM_SQRT) {
if (isinf(sexp_flonum_value(res)))
res = sexp_bignum_sqrt_estimate(ctx, a);
else
res = sexp_double_to_bignum(ctx, sexp_flonum_value(res));
loop: /* until 0 <= a - res*res < 2*res + 1 */
rem = sexp_mul(ctx, res, res);
tmp = rem = sexp_sub(ctx, a, rem);
if (!sexp_positivep(tmp)) goto adjust;
tmp = sexp_sub(ctx, tmp, SEXP_ONE);
tmp = sexp_quotient(ctx, tmp, SEXP_TWO);
tmp = sexp_compare(ctx, tmp, res);
if (sexp_positivep(tmp)) {
adjust:
tmp = sexp_quotient(ctx, a, res);
res = sexp_add(ctx, res, tmp);
res = sexp_quotient(ctx, res, SEXP_TWO);
goto loop;
}
} else {
if (sexp_flonump(res))
res = sexp_bignum_normalize(sexp_double_to_bignum(ctx, sexp_flonum_value(res)));
tmp = sexp_mul(ctx, res, res);
rem = sexp_sub(ctx, a, tmp);
}
*rem_out = sexp_bignum_normalize(rem);
sexp_gc_release4(ctx);
return sexp_bignum_normalize(res);
}
#endif /* SEXP_USE_MATH */
/************************ ratios ******************************/
#if SEXP_USE_RATIOS
double sexp_ratio_to_double (sexp rat) {
sexp num = sexp_ratio_numerator(rat), den = sexp_ratio_denominator(rat);
return (sexp_bignump(num) ? sexp_bignum_to_double(num)
: sexp_fixnum_to_double(num))
/ (sexp_bignump(den) ? sexp_bignum_to_double(den)
: sexp_fixnum_to_double(den));
}
sexp sexp_double_to_ratio (sexp ctx, double f) {
int sign, i;
sexp_gc_var3(res, whole, scale);
if (f == trunc(f))
return sexp_bignum_normalize(sexp_double_to_bignum(ctx, f));
sexp_gc_preserve3(ctx, res, whole, scale);
whole = sexp_double_to_bignum(ctx, trunc(f));
res = sexp_fixnum_to_bignum(ctx, SEXP_ZERO);
scale = SEXP_ONE;
sign = (f < 0 ? -1 : 1);
for (i=0, f=fabs(f-trunc(f)); f != trunc(f) && i < 15; i++) {
res = sexp_bignum_fxmul(ctx, NULL, res, 10, 0);
f = f * 10;
res = sexp_bignum_fxadd(ctx, res, (sexp_uint_t)double_10s_digit(f));
f = f - trunc(f);
scale = sexp_mul(ctx, scale, SEXP_TEN);
}
sexp_bignum_sign(res) = sign;
res = sexp_bignum_normalize(res);
scale = sexp_bignum_normalize(scale);
res = sexp_make_ratio(ctx, res, scale);
res = sexp_ratio_normalize(ctx, res, SEXP_FALSE);
res = sexp_add(ctx, res, whole);
sexp_gc_release3(ctx);
return res;
}
//
// For conversion that does not introduce round-off error,
// no matter what FLT_RADIX is.
//
sexp sexp_double_to_ratio_2 (sexp ctx, double f) {
int sign,i;
sexp_gc_var3(res, whole, scale);
if (f == trunc(f))
return sexp_bignum_normalize(sexp_double_to_bignum(ctx, f));
sexp_gc_preserve3(ctx, res, whole, scale);
whole = sexp_double_to_bignum(ctx, trunc(f));
res = sexp_fixnum_to_bignum(ctx, SEXP_ZERO);
scale = SEXP_ONE;
sign = (f < 0 ? -1 : 1);
f = fabs(f-trunc(f));
while(f) {
res = sexp_bignum_fxmul(ctx, NULL, res, FLT_RADIX, 0);
scale = sexp_mul(ctx, scale, sexp_make_fixnum(FLT_RADIX));
f *= FLT_RADIX;
i = trunc(f);
if (i) {
f -= i;
res = sexp_bignum_fxadd(ctx, res, i);
}
}
sexp_bignum_sign(res) = sign;
res = sexp_bignum_normalize(res);
scale = sexp_bignum_normalize(scale);
res = sexp_make_ratio(ctx, res, scale);
res = sexp_ratio_normalize(ctx, res, SEXP_FALSE);
res = sexp_add(ctx, res, whole);
sexp_gc_release3(ctx);
return res;
}
sexp sexp_ratio_add (sexp ctx, sexp a, sexp b) {
sexp_gc_var3(res, num, den);
sexp_gc_preserve3(ctx, res, num, den);
num = sexp_mul(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(b));
den = sexp_mul(ctx, sexp_ratio_numerator(b), sexp_ratio_denominator(a));
num = sexp_add(ctx, num, den);
den = sexp_mul(ctx, sexp_ratio_denominator(a), sexp_ratio_denominator(b));
res = sexp_make_ratio(ctx, num, den);
res = sexp_ratio_normalize(ctx, res, SEXP_FALSE);
sexp_gc_release3(ctx);
return res;
}
sexp sexp_ratio_mul (sexp ctx, sexp a, sexp b) {
sexp_gc_var3(res, num, den);
sexp_gc_preserve3(ctx, res, num, den);
num = sexp_mul(ctx, sexp_ratio_numerator(a), sexp_ratio_numerator(b));
den = sexp_mul(ctx, sexp_ratio_denominator(a), sexp_ratio_denominator(b));
res = sexp_make_ratio(ctx, num, den);
res = sexp_ratio_normalize(ctx, res, SEXP_FALSE);
sexp_gc_release3(ctx);
return res;
}
sexp sexp_ratio_div (sexp ctx, sexp a, sexp b) {
sexp_gc_var3(res, num, den);
sexp_gc_preserve3(ctx, res, num, den);
num = sexp_mul(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(b));
den = sexp_mul(ctx, sexp_ratio_denominator(a), sexp_ratio_numerator(b));
res = sexp_make_ratio(ctx, num, den);
res = sexp_ratio_normalize(ctx, res, SEXP_FALSE);
sexp_gc_release3(ctx);
return res;
}
sexp sexp_ratio_compare (sexp ctx, sexp a, sexp b) {
sexp_gc_var2(a2, b2);
sexp_gc_preserve2(ctx, a2, b2);
a2 = sexp_mul(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(b));
b2 = sexp_mul(ctx, sexp_ratio_numerator(b), sexp_ratio_denominator(a));
a2 = sexp_compare(ctx, a2, b2);
sexp_gc_release2(ctx);
return a2;
}
sexp sexp_ratio_round (sexp ctx, sexp a) {
sexp_gc_var2(q, r);
sexp_gc_preserve2(ctx, q, r);
q = sexp_quotient(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(a));
if ((sexp_ratio_denominator(a) == SEXP_TWO) && sexp_oddp(q)) {
q = sexp_add(ctx, q, (sexp_exact_positivep(q) ? SEXP_ONE : SEXP_NEG_ONE));
} else {
r = sexp_remainder(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(a));
r = sexp_mul(ctx, r, SEXP_TWO);
if (sexp_exact_negativep(r)) {sexp_negate(r);}
if (sexp_unbox_fixnum(sexp_compare(ctx, r, sexp_ratio_denominator(a))) > 0)
q = sexp_add(ctx, q, (sexp_exact_positivep(q) ? SEXP_ONE : SEXP_NEG_ONE));
}
sexp_gc_release2(ctx);
return q;
}
sexp sexp_ratio_trunc (sexp ctx, sexp a) {
return sexp_quotient(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(a));
}
sexp sexp_ratio_floor (sexp ctx, sexp a) {
sexp_gc_var1(q);
sexp_gc_preserve1(ctx, q);
q = sexp_quotient(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(a));
if (sexp_exact_negativep(sexp_ratio_numerator(a)))
q = sexp_add(ctx, q, SEXP_NEG_ONE);
sexp_gc_release1(ctx);
return q;
}
sexp sexp_ratio_ceiling (sexp ctx, sexp a) {
sexp_gc_var1(q);
sexp_gc_preserve1(ctx, q);
q = sexp_quotient(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(a));
if (sexp_exact_positivep(sexp_ratio_numerator(a)))
q = sexp_add(ctx, q, SEXP_ONE);
sexp_gc_release1(ctx);
return q;
}
#endif
double sexp_to_double (sexp x) {
if (sexp_flonump(x))
return sexp_flonum_value(x);
else if (sexp_fixnump(x))
return sexp_fixnum_to_double(x);
else if (sexp_bignump(x))
return sexp_bignum_to_double(x);
#if SEXP_USE_RATIOS
else if (sexp_ratiop(x))
return sexp_ratio_to_double(x);
#endif
else
return 0.0;
}
/************************ complex numbers ****************************/
#if SEXP_USE_COMPLEX
static sexp sexp_complex_copy (sexp ctx, sexp a) {
sexp_gc_var1(res);
sexp_gc_preserve1(ctx, res);
res = sexp_make_complex(ctx, sexp_complex_real(a), sexp_complex_imag(a));
if (sexp_flonump(sexp_complex_real(a)))
sexp_complex_real(a) = sexp_make_flonum(ctx, sexp_flonum_value(sexp_complex_real(a)));
else if (sexp_bignump(sexp_complex_real(a)))
sexp_complex_real(a) = sexp_copy_bignum(ctx, NULL, sexp_complex_real(a), 0);
if (sexp_flonump(sexp_complex_imag(a)))
sexp_complex_imag(a) = sexp_make_flonum(ctx, sexp_flonum_value(sexp_complex_imag(a)));
else if (sexp_bignump(sexp_complex_imag(a)))
sexp_complex_imag(a) = sexp_copy_bignum(ctx, NULL, sexp_complex_imag(a), 0);
sexp_gc_release1(ctx);
return res;
}
sexp sexp_complex_add (sexp ctx, sexp a, sexp b) {
sexp_gc_var3(res, real, imag);
sexp_gc_preserve3(ctx, res, real, imag);
real = sexp_add(ctx, sexp_complex_real(a), sexp_complex_real(b));
imag = sexp_add(ctx, sexp_complex_imag(a), sexp_complex_imag(b));
res = sexp_make_complex(ctx, real, imag);
sexp_gc_release3(ctx);
return sexp_complex_normalize(res);
}
sexp sexp_complex_sub (sexp ctx, sexp a, sexp b) {
sexp_gc_var2(res, tmp);
sexp_gc_preserve2(ctx, res, tmp);
tmp = sexp_complex_copy(ctx, b);
sexp_negate(sexp_complex_real(tmp));
sexp_negate(sexp_complex_imag(tmp));
res = sexp_complex_add(ctx, a, tmp);
sexp_gc_release2(ctx);
return res;
}
sexp sexp_complex_mul (sexp ctx, sexp a, sexp b) {
sexp_gc_var3(res, real, imag);
sexp_gc_preserve3(ctx, res, real, imag);
real = sexp_mul(ctx, sexp_complex_real(a), sexp_complex_real(b));
res = sexp_mul(ctx, sexp_complex_imag(a), sexp_complex_imag(b));
real = sexp_sub(ctx, real, res);
imag = sexp_mul(ctx, sexp_complex_real(a), sexp_complex_imag(b));
res = sexp_mul(ctx, sexp_complex_imag(a), sexp_complex_real(b));
imag = sexp_add(ctx, imag, res);
res = sexp_make_complex(ctx, real, imag);
sexp_gc_release3(ctx);
return sexp_complex_normalize(res);
}
/* (a + bi) (ac + bd) (bc - ad) */
/* -------- = ----------- + ----------- i */
/* (c + di) (c^2 + d^2) (c^2 + d^2) */
sexp sexp_complex_div (sexp ctx, sexp a, sexp b) {
sexp_gc_var4(res, real, imag, denom);
sexp_gc_preserve4(ctx, res, real, imag, denom);
/* c^2 + d^2 */
denom = sexp_mul(ctx, sexp_complex_real(b), sexp_complex_real(b));
res = sexp_mul(ctx, sexp_complex_imag(b), sexp_complex_imag(b));
denom = sexp_add(ctx, denom, res);
/* ac + bd */
real = sexp_mul(ctx, sexp_complex_real(a), sexp_complex_real(b));
res = sexp_mul(ctx, sexp_complex_imag(a), sexp_complex_imag(b));
real = sexp_add(ctx, real, res);
real = sexp_div(ctx, real, denom);
/* bc - ad */
imag = sexp_mul(ctx, sexp_complex_imag(a), sexp_complex_real(b));
res = sexp_mul(ctx, sexp_complex_real(a), sexp_complex_imag(b));
imag = sexp_sub(ctx, imag, res);
imag = sexp_div(ctx, imag, denom);
res = sexp_make_complex(ctx, real, imag);
sexp_gc_release4(ctx);
return sexp_complex_normalize(res);
}
static sexp sexp_to_complex (sexp ctx, sexp x) {
#if SEXP_USE_RATIOS
sexp_gc_var1(tmp);
#endif
if (sexp_flonump(x) || sexp_fixnump(x) || sexp_bignump(x)) {
return sexp_make_complex(ctx, x, SEXP_ZERO);
#if SEXP_USE_RATIOS
} else if (sexp_ratiop(x)) {
sexp_gc_preserve1(ctx, tmp);
tmp = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
sexp_complex_real(tmp) = sexp_make_flonum(ctx, sexp_to_double(x));
sexp_gc_release1(ctx);
return tmp;
#endif
} else {
return x;
}
}
sexp sexp_complex_exp (sexp ctx, sexp z) {
double e2x = exp(sexp_to_double(sexp_complex_real(z))),
y = sexp_to_double(sexp_complex_imag(z));
sexp_gc_var1(res);
sexp_gc_preserve1(ctx, res);
res = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
sexp_complex_real(res) = sexp_make_flonum(ctx, e2x*cos(y));
sexp_complex_imag(res) = sexp_make_flonum(ctx, e2x*sin(y));
sexp_gc_release1(ctx);
return res;
}
sexp sexp_complex_log (sexp ctx, sexp z) {
double x = sexp_to_double(sexp_complex_real(z)),
y = sexp_to_double(sexp_complex_imag(z));
sexp_gc_var1(res);
sexp_gc_preserve1(ctx, res);
res = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
sexp_complex_real(res) = sexp_make_flonum(ctx, log(sqrt(x*x + y*y)));
sexp_complex_imag(res) = sexp_make_flonum(ctx, atan2(y, x));
sexp_gc_release1(ctx);
return res;
}
sexp sexp_complex_expt (sexp ctx, sexp a, sexp b) {
sexp_gc_var1(res);
sexp_gc_preserve1(ctx, res);
res = sexp_to_complex(ctx, a);
res = sexp_complex_log(ctx, res);
res = sexp_mul(ctx, b, res);
res = sexp_complex_exp(ctx, res);
sexp_gc_release1(ctx);
return res;
}
#if SEXP_USE_MATH
sexp sexp_complex_sqrt (sexp ctx, sexp z) {
double x = sexp_to_double(sexp_complex_real(z)),
y = sexp_to_double(sexp_complex_imag(z)), r;
sexp_gc_var1(res);
sexp_gc_preserve1(ctx, res);
r = sqrt(x*x + y*y);
res = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
sexp_complex_real(res) = sexp_make_flonum(ctx, sqrt((x+r)/2));
sexp_complex_imag(res) = sexp_make_flonum(ctx, ((y<0||(y==0&&1/y<0))?-1:1)*sqrt((-x+r)/2));
sexp_gc_release1(ctx);
return res;
}
sexp sexp_complex_sin (sexp ctx, sexp z) {
double x = sexp_to_double(sexp_complex_real(z)),
y = sexp_to_double(sexp_complex_imag(z));
sexp_gc_var1(res);
sexp_gc_preserve1(ctx, res);
res = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
sexp_complex_real(res) = sexp_make_flonum(ctx, sin(x)*cosh(y));
sexp_complex_imag(res) = sexp_make_flonum(ctx, cos(x)*sinh(y));
sexp_gc_release1(ctx);
return res;
}
sexp sexp_complex_cos (sexp ctx, sexp z) {
double x = sexp_to_double(sexp_complex_real(z)),
y = sexp_to_double(sexp_complex_imag(z));
sexp_gc_var1(res);
sexp_gc_preserve1(ctx, res);
res = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
sexp_complex_real(res) = sexp_make_flonum(ctx, cos(x)*cosh(y));
sexp_complex_imag(res) = sexp_make_flonum(ctx, -sin(x)*sinh(y));
sexp_gc_release1(ctx);
return res;
}
sexp sexp_complex_tan (sexp ctx, sexp z) {
sexp res;
sexp_gc_var2(sin, cos);
sexp_gc_preserve2(ctx, sin, cos);
sin = sexp_complex_sin(ctx, z);
cos = sexp_complex_cos(ctx, z);
res = sexp_complex_div(ctx, sin, cos);
sexp_gc_release2(ctx);
return res;
}
sexp sexp_complex_asin (sexp ctx, sexp z) {
sexp_gc_var3(res, tmp, tmp2);
sexp_gc_preserve3(ctx, res, tmp, tmp2);
res = sexp_complex_mul(ctx, z, z);
res = sexp_sub(ctx, SEXP_ONE, res);
res = sexp_sqrt(ctx, NULL, 1, res);
tmp = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
sexp_complex_real(tmp) = sexp_mul(ctx, SEXP_NEG_ONE, sexp_complex_imag(z));
sexp_complex_imag(tmp) = sexp_complex_real(z);
res = sexp_add(ctx, tmp, res);
res = sexp_log(ctx, NULL, 1, res);
tmp = sexp_make_complex(ctx, SEXP_ZERO, SEXP_NEG_ONE);
res = sexp_mul(ctx, res, tmp);
sexp_gc_release3(ctx);
return res;
}
sexp sexp_complex_acos (sexp ctx, sexp z) {
sexp_gc_var2(res, tmp);
sexp_gc_preserve2(ctx, res, tmp);
res = sexp_complex_asin(ctx, z);
tmp = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
sexp_complex_real(tmp) = sexp_make_flonum(ctx, acos(-1)/2);
res = sexp_sub(ctx, tmp, res);
sexp_gc_release2(ctx);
return res;
}
sexp sexp_complex_atan (sexp ctx, sexp z) {
sexp_gc_var3(res, tmp1, tmp2);
sexp_gc_preserve3(ctx, res, tmp1, tmp2);
tmp1 = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ONE);
tmp1 = sexp_complex_mul(ctx, z, tmp1);
res = sexp_make_complex(ctx, SEXP_ONE, SEXP_ZERO);
res = sexp_complex_sub(ctx, res, tmp1);
res = sexp_complex_log(ctx, res);
tmp2 = sexp_make_complex(ctx, SEXP_ONE, SEXP_ZERO);
tmp2 = sexp_complex_add(ctx, tmp2, tmp1);
tmp2 = sexp_complex_log(ctx, tmp2);
res = sexp_complex_sub(ctx, res, tmp2);
tmp1 = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ONE);
sexp_complex_imag(tmp1) = sexp_make_flonum(ctx, 0.5);
res = sexp_complex_mul(ctx, res, tmp1);
sexp_gc_release3(ctx);
return res;
}
#endif
#endif
/****************** generic arithmetic ************************/
enum sexp_number_types {
SEXP_NUM_NOT = 0,
SEXP_NUM_FIX,
SEXP_NUM_FLO,
SEXP_NUM_BIG,
#if SEXP_USE_RATIOS
SEXP_NUM_RAT,
#endif
#if SEXP_USE_COMPLEX
SEXP_NUM_CPX,
#endif
};
enum sexp_number_combs {
SEXP_NUM_NOT_NOT = 0,
SEXP_NUM_NOT_FIX,
SEXP_NUM_NOT_FLO,
SEXP_NUM_NOT_BIG,
#if SEXP_USE_RATIOS
SEXP_NUM_NOT_RAT,
#endif
#if SEXP_USE_COMPLEX
SEXP_NUM_NOT_CPX,
#endif
SEXP_NUM_FIX_NOT,
SEXP_NUM_FIX_FIX,
SEXP_NUM_FIX_FLO,
SEXP_NUM_FIX_BIG,
#if SEXP_USE_RATIOS
SEXP_NUM_FIX_RAT,
#endif
#if SEXP_USE_COMPLEX
SEXP_NUM_FIX_CPX,
#endif
SEXP_NUM_FLO_NOT,
SEXP_NUM_FLO_FIX,
SEXP_NUM_FLO_FLO,
SEXP_NUM_FLO_BIG,
#if SEXP_USE_RATIOS
SEXP_NUM_FLO_RAT,
#endif
#if SEXP_USE_COMPLEX
SEXP_NUM_FLO_CPX,
#endif
SEXP_NUM_BIG_NOT,
SEXP_NUM_BIG_FIX,
SEXP_NUM_BIG_FLO,
SEXP_NUM_BIG_BIG,
#if SEXP_USE_RATIOS
SEXP_NUM_BIG_RAT,
#endif
#if SEXP_USE_COMPLEX
SEXP_NUM_BIG_CPX,
#endif
#if SEXP_USE_RATIOS
SEXP_NUM_RAT_NOT,
SEXP_NUM_RAT_FIX,
SEXP_NUM_RAT_FLO,
SEXP_NUM_RAT_BIG,
SEXP_NUM_RAT_RAT,
#if SEXP_USE_COMPLEX
SEXP_NUM_RAT_CPX,
#endif
#endif
#if SEXP_USE_COMPLEX
SEXP_NUM_CPX_NOT,
SEXP_NUM_CPX_FIX,
SEXP_NUM_CPX_FLO,
SEXP_NUM_CPX_BIG,
#if SEXP_USE_RATIOS
SEXP_NUM_CPX_RAT,
#endif
SEXP_NUM_CPX_CPX,
#endif
};
static int sexp_number_types[] =
#if SEXP_USE_RATIOS && SEXP_USE_COMPLEX
{0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 5, 0, 0};
#else
#if SEXP_USE_RATIOS || SEXP_USE_COMPLEX
{0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 0, 0, 0};
#else
{0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 0, 0, 0, 0};
#endif
#endif
#define SEXP_NUM_NUMBER_TYPES (4 + SEXP_USE_RATIOS + SEXP_USE_COMPLEX)
static int sexp_number_type (sexp a) {
return sexp_pointerp(a) ?
(sexp_pointer_tag(a)<(sizeof(sexp_number_types)/sizeof(sexp_number_types[0]))
? sexp_number_types[sexp_pointer_tag(a)] : 0)
#if SEXP_USE_IMMEDIATE_FLONUMS
: sexp_flonump(a) ? 2
#endif
: sexp_fixnump(a);
}
sexp sexp_add (sexp ctx, sexp a, sexp b) {
sexp_sint_t sum;
int at=sexp_number_type(a), bt=sexp_number_type(b), t;
sexp r=SEXP_VOID;
sexp_gc_var1(tmp);
sexp_gc_preserve1(ctx, tmp);
if (at > bt) {r=a; a=b; b=r; t=at; at=bt; bt=t;}
switch ((at * SEXP_NUM_NUMBER_TYPES) + bt) {
case SEXP_NUM_NOT_NOT: case SEXP_NUM_NOT_FIX:
case SEXP_NUM_NOT_FLO: case SEXP_NUM_NOT_BIG:
#if SEXP_USE_RATIOS
case SEXP_NUM_NOT_RAT:
#endif
#if SEXP_USE_COMPLEX
case SEXP_NUM_NOT_CPX:
#endif
r = sexp_type_exception(ctx, NULL, SEXP_NUMBER, a);
break;
case SEXP_NUM_FIX_FIX:
sum = sexp_unbox_fixnum(a) + sexp_unbox_fixnum(b);
if ((sum < SEXP_MIN_FIXNUM) || (sum > SEXP_MAX_FIXNUM))
r = sexp_add(ctx, tmp=sexp_fixnum_to_bignum(ctx, a), b);
else
r = sexp_make_fixnum(sum);
break;
case SEXP_NUM_FIX_FLO:
r = a == SEXP_ZERO ? b : sexp_make_flonum(ctx, sexp_fixnum_to_double(a)+sexp_flonum_value(b));
break;
case SEXP_NUM_FIX_BIG:
r = sexp_bignum_normalize(sexp_bignum_add_fixnum(ctx, b, a));
break;
case SEXP_NUM_FLO_FLO:
r = sexp_fp_add(ctx, a, b);
break;
case SEXP_NUM_FLO_BIG:
r = sexp_make_flonum(ctx, sexp_flonum_value(a) + sexp_bignum_to_double(b));
break;
case SEXP_NUM_BIG_BIG:
r = sexp_bignum_normalize(sexp_bignum_add(ctx, NULL, b, a));
break;
#if SEXP_USE_RATIOS
case SEXP_NUM_FLO_RAT:
r = sexp_make_flonum(ctx, sexp_flonum_value(a) + sexp_ratio_to_double(b));
break;
case SEXP_NUM_FIX_RAT:
case SEXP_NUM_BIG_RAT:
a = tmp = sexp_make_ratio(ctx, a, SEXP_ONE);
/* ... FALLTHROUGH ... */
case SEXP_NUM_RAT_RAT:
r = sexp_ratio_add(ctx, a, b);
break;
#endif
#if SEXP_USE_COMPLEX
#if SEXP_USE_RATIOS
case SEXP_NUM_RAT_CPX:
#endif
case SEXP_NUM_FLO_CPX:
case SEXP_NUM_FIX_CPX:
case SEXP_NUM_BIG_CPX:
a = tmp = sexp_make_complex(ctx, a, SEXP_ZERO);
/* ... FALLTHROUGH ... */
case SEXP_NUM_CPX_CPX:
r = sexp_complex_add(ctx, a, b);
break;
#endif
}
sexp_gc_release1(ctx);
return r;
}
sexp sexp_sub (sexp ctx, sexp a, sexp b) {
#if SEXP_USE_FLONUMS
int negatep=0;
#endif
int at=sexp_number_type(a), bt=sexp_number_type(b);
sexp r=SEXP_VOID;
sexp_gc_var2(tmp1, tmp2);
sexp_gc_preserve2(ctx, tmp1, tmp2);
switch ((at * SEXP_NUM_NUMBER_TYPES) + bt) {
case SEXP_NUM_NOT_NOT: case SEXP_NUM_NOT_FIX:
case SEXP_NUM_NOT_FLO: case SEXP_NUM_NOT_BIG:
#if SEXP_USE_RATIOS
case SEXP_NUM_NOT_RAT:
#endif
#if SEXP_USE_COMPLEX
case SEXP_NUM_NOT_CPX:
#endif
r = sexp_type_exception(ctx, NULL, SEXP_NUMBER, a);
break;
case SEXP_NUM_FIX_NOT: case SEXP_NUM_FLO_NOT: case SEXP_NUM_BIG_NOT:
#if SEXP_USE_RATIOS
case SEXP_NUM_RAT_NOT:
#endif
#if SEXP_USE_COMPLEX
case SEXP_NUM_CPX_NOT:
#endif
r = sexp_type_exception(ctx, NULL, SEXP_NUMBER, b);
break;
case SEXP_NUM_FIX_FIX:
r = sexp_fx_sub(a, b); /* VM catches this case */
break;
case SEXP_NUM_FIX_FLO:
r = sexp_make_flonum(ctx, sexp_fixnum_to_double(a)-sexp_flonum_value(b));
break;
case SEXP_NUM_FIX_BIG:
tmp1 = sexp_fixnum_to_bignum(ctx, a);
r = sexp_bignum_sub(ctx, NULL, b, tmp1);
sexp_negate_exact(r);
r = sexp_bignum_normalize(r);
break;
case SEXP_NUM_FLO_FIX:
r = sexp_make_flonum(ctx, sexp_flonum_value(a) - sexp_fixnum_to_double(b));
break;
case SEXP_NUM_FLO_FLO:
r = sexp_fp_sub(ctx, a, b);
break;
case SEXP_NUM_FLO_BIG:
r = sexp_make_flonum(ctx, sexp_flonum_value(a) - sexp_bignum_to_double(b));
break;
case SEXP_NUM_BIG_FIX:
tmp1 = sexp_fixnum_to_bignum(ctx, b);
r = sexp_bignum_normalize(sexp_bignum_sub(ctx, NULL, a, tmp1));
break;
case SEXP_NUM_BIG_FLO:
r = sexp_make_flonum(ctx, sexp_bignum_to_double(a) - sexp_flonum_value(b));
break;
case SEXP_NUM_BIG_BIG:
r = sexp_bignum_normalize(sexp_bignum_sub(ctx, NULL, a, b));
break;
#if SEXP_USE_RATIOS
case SEXP_NUM_FLO_RAT:
r = sexp_make_flonum(ctx, sexp_flonum_value(a) - sexp_ratio_to_double(b));
break;
case SEXP_NUM_RAT_FLO:
r = sexp_make_flonum(ctx, sexp_ratio_to_double(a) - sexp_flonum_value(b));
break;
case SEXP_NUM_RAT_FIX:
case SEXP_NUM_RAT_BIG:
tmp1 = a; a = b; b = tmp1;
negatep = 1;
/* ... FALLTHROUGH ... */
case SEXP_NUM_FIX_RAT:
case SEXP_NUM_BIG_RAT:
a = tmp1 = sexp_make_ratio(ctx, a, SEXP_ONE);
/* ... FALLTHROUGH ... */
case SEXP_NUM_RAT_RAT:
tmp2 = sexp_make_ratio(ctx, sexp_ratio_numerator(b), sexp_ratio_denominator(b));
if (sexp_bignump(sexp_ratio_numerator(tmp2)))
sexp_ratio_numerator(tmp2) = sexp_copy_bignum(ctx, NULL, sexp_ratio_numerator(tmp2), 0);
sexp_negate_exact(sexp_ratio_numerator(tmp2));
r = sexp_ratio_add(ctx, a, tmp2);
if (negatep) {
if (sexp_ratiop(r)) {
sexp_negate_exact(sexp_ratio_numerator(r));
} else {
sexp_negate_exact(r);
}
}
break;
#endif
#if SEXP_USE_COMPLEX
#if SEXP_USE_RATIOS
case SEXP_NUM_RAT_CPX:
a = tmp1 = sexp_make_flonum(ctx, sexp_ratio_to_double(a));
goto complex_sub;
case SEXP_NUM_CPX_RAT:
b = tmp1 = sexp_make_flonum(ctx, sexp_ratio_to_double(b));
/* ... FALLTHROUGH ... */
#endif
case SEXP_NUM_CPX_FLO:
case SEXP_NUM_CPX_FIX:
case SEXP_NUM_CPX_BIG:
tmp1 = a; a = b; b = tmp1;
negatep = 1;
/* ... FALLTHROUGH ... */
case SEXP_NUM_FLO_CPX:
case SEXP_NUM_FIX_CPX:
case SEXP_NUM_BIG_CPX:
a = tmp1 = sexp_make_complex(ctx, a, SEXP_ZERO);
/* ... FALLTHROUGH ... */
case SEXP_NUM_CPX_CPX:
#if SEXP_USE_RATIOS
complex_sub:
#endif
r = sexp_complex_sub(ctx, a, b);
if (negatep) {
if (sexp_complexp(r)) {
r = sexp_complex_copy(ctx, r);
sexp_negate(sexp_complex_real(r));
sexp_negate(sexp_complex_imag(r));
} else {
sexp_negate(r);
}
}
break;
#endif
}
sexp_gc_release2(ctx);
return r;
}
sexp sexp_mul (sexp ctx, sexp a, sexp b) {
sexp_lsint_t prod;
int at=sexp_number_type(a), bt=sexp_number_type(b), t;
sexp r=SEXP_VOID;
sexp_gc_var1(tmp);
sexp_gc_preserve1(ctx, tmp);
if (at > bt) {r=a; a=b; b=r; t=at; at=bt; bt=t;}
switch ((at * SEXP_NUM_NUMBER_TYPES) + bt) {
case SEXP_NUM_NOT_NOT: case SEXP_NUM_NOT_FIX:
case SEXP_NUM_NOT_FLO: case SEXP_NUM_NOT_BIG:
#if SEXP_USE_RATIOS
case SEXP_NUM_NOT_RAT:
#endif
r = sexp_type_exception(ctx, NULL, SEXP_NUMBER, a);
break;
case SEXP_NUM_FIX_FIX:
prod = lsint_mul_sint(lsint_from_sint(sexp_unbox_fixnum(a)), sexp_unbox_fixnum(b));
if (!lsint_is_fixnum(prod))
r = sexp_mul(ctx, tmp=sexp_fixnum_to_bignum(ctx, a), b);
else
r = sexp_make_fixnum(lsint_to_sint(prod));
break;
case SEXP_NUM_FIX_FLO:
r = (a==SEXP_ZERO ? a : sexp_make_flonum(ctx, sexp_fixnum_to_double(a)*sexp_flonum_value(b)));
break;
case SEXP_NUM_FIX_BIG:
r = sexp_bignum_fxmul(ctx, NULL, b, sexp_unbox_fx_abs(a), 0);
sexp_bignum_sign(r) = sexp_fx_sign(a) * sexp_bignum_sign(b);
r = sexp_bignum_normalize(r);
break;
case SEXP_NUM_FLO_FLO:
r = sexp_fp_mul(ctx, a, b);
break;
case SEXP_NUM_FLO_BIG:
r = sexp_make_flonum(ctx, sexp_flonum_value(a) * sexp_bignum_to_double(b));
break;
case SEXP_NUM_BIG_BIG:
r = sexp_bignum_normalize(sexp_bignum_mul(ctx, NULL, a, b));
break;
#if SEXP_USE_RATIOS
case SEXP_NUM_FLO_RAT:
r = sexp_make_flonum(ctx, sexp_flonum_value(a) * sexp_ratio_to_double(b));
break;
case SEXP_NUM_FIX_RAT:
case SEXP_NUM_BIG_RAT:
a = tmp = sexp_make_ratio(ctx, a, SEXP_ONE);
/* ... FALLTHROUGH ... */
case SEXP_NUM_RAT_RAT:
r = sexp_ratio_mul(ctx, a, b);
break;
#endif
#if SEXP_USE_COMPLEX
#if SEXP_USE_RATIOS
case SEXP_NUM_RAT_CPX:
#endif
case SEXP_NUM_FLO_CPX:
case SEXP_NUM_FIX_CPX:
case SEXP_NUM_BIG_CPX:
a = tmp = sexp_make_complex(ctx, a, SEXP_ZERO);
/* ... FALLTHROUGH ... */
case SEXP_NUM_CPX_CPX:
r = sexp_complex_mul(ctx, a, b);
break;
#endif
}
sexp_gc_release1(ctx);
return r;
}
sexp sexp_div (sexp ctx, sexp a, sexp b) {
int at=sexp_number_type(a), bt=sexp_number_type(b);
#if ! SEXP_USE_RATIOS
double f;
#endif
sexp r=SEXP_VOID;
sexp_gc_var1(tmp);
sexp_gc_preserve1(ctx, tmp);
switch ((at * SEXP_NUM_NUMBER_TYPES) + bt) {
case SEXP_NUM_NOT_NOT: case SEXP_NUM_NOT_FIX:
case SEXP_NUM_NOT_FLO: case SEXP_NUM_NOT_BIG:
#if SEXP_USE_RATIOS
case SEXP_NUM_NOT_RAT:
#endif
r = sexp_type_exception(ctx, NULL, SEXP_NUMBER, a);
break;
case SEXP_NUM_FIX_NOT: case SEXP_NUM_FLO_NOT: case SEXP_NUM_BIG_NOT:
#if SEXP_USE_RATIOS
case SEXP_NUM_RAT_NOT:
#endif
r = sexp_type_exception(ctx, NULL, SEXP_NUMBER, b);
break;
case SEXP_NUM_FIX_FIX:
#if SEXP_USE_RATIOS
tmp = sexp_make_ratio(ctx, a, b);
r = sexp_ratio_normalize(ctx, tmp, SEXP_FALSE);
#else
f = sexp_fixnum_to_double(a) / sexp_fixnum_to_double(b);
r = ((f == trunc(f)) ? sexp_make_fixnum((sexp_sint_t)f)
: sexp_make_flonum(ctx, f));
#endif
break;
case SEXP_NUM_FIX_FLO:
r = sexp_make_flonum(ctx, sexp_fixnum_to_double(a)/sexp_flonum_value(b));
break;
case SEXP_NUM_FIX_BIG:
#if SEXP_USE_RATIOS
tmp = sexp_make_ratio(ctx, a, b);
r = sexp_ratio_normalize(ctx, tmp, SEXP_FALSE);
#else
r = sexp_make_flonum(ctx, sexp_fixnum_to_double(a)/sexp_bignum_to_double(b));
#endif
break;
case SEXP_NUM_FLO_FIX:
r = sexp_make_flonum(ctx, sexp_flonum_value(a)/sexp_fixnum_to_double(b));
break;
case SEXP_NUM_FLO_FLO:
r = sexp_fp_div(ctx, a, b);
break;
case SEXP_NUM_FLO_BIG:
r = sexp_make_flonum(ctx, sexp_flonum_value(a) / sexp_bignum_to_double(b));
break;
case SEXP_NUM_BIG_FIX:
#if SEXP_USE_RATIOS
tmp = sexp_make_ratio(ctx, a, b);
r = sexp_ratio_normalize(ctx, tmp, SEXP_FALSE);
break;
#else
b = tmp = sexp_fixnum_to_bignum(ctx, b);
#endif
/* ... FALLTHROUGH if ! SEXP_USE_RATIOS ... */
case SEXP_NUM_BIG_BIG:
#if SEXP_USE_RATIOS
tmp = sexp_make_ratio(ctx, a, b);
r = sexp_ratio_normalize(ctx, tmp, SEXP_FALSE);
#else
r = sexp_bignum_quot_rem(ctx, &tmp, a, b);
if (sexp_bignum_normalize(tmp) != SEXP_ZERO)
r = sexp_make_flonum(ctx, sexp_bignum_to_double(a)
/ sexp_bignum_to_double(b));
else
r = sexp_bignum_normalize(r);
#endif
break;
case SEXP_NUM_BIG_FLO:
r = sexp_make_flonum(ctx, sexp_bignum_to_double(a) / sexp_flonum_value(b));
break;
#if SEXP_USE_RATIOS
case SEXP_NUM_FLO_RAT:
r = sexp_make_flonum(ctx, sexp_flonum_value(a) / sexp_ratio_to_double(b));
break;
case SEXP_NUM_RAT_FLO:
r = sexp_make_flonum(ctx, sexp_ratio_to_double(a) / sexp_flonum_value(b));
break;
case SEXP_NUM_RAT_FIX:
case SEXP_NUM_RAT_BIG:
b = tmp = sexp_make_ratio(ctx, b, SEXP_ONE);
/* ... FALLTHROUGH ... */
case SEXP_NUM_FIX_RAT:
case SEXP_NUM_BIG_RAT:
if (!sexp_ratiop(a))
a = tmp = sexp_make_ratio(ctx, a, SEXP_ONE);
/* ... FALLTHROUGH ... */
case SEXP_NUM_RAT_RAT:
r = sexp_ratio_div(ctx, a, b);
break;
#endif
#if SEXP_USE_COMPLEX
#if SEXP_USE_RATIOS
case SEXP_NUM_CPX_RAT:
b = tmp = sexp_make_flonum(ctx, sexp_ratio_to_double(b));
/* ... FALLTHROUGH ... */
#endif
case SEXP_NUM_CPX_FLO:
case SEXP_NUM_CPX_FIX:
case SEXP_NUM_CPX_BIG:
b = tmp = sexp_make_complex(ctx, b, SEXP_ZERO);
/* ... FALLTHROUGH ... */
#if SEXP_USE_RATIOS
case SEXP_NUM_RAT_CPX:
if (sexp_ratiop(a))
a = tmp = sexp_make_flonum(ctx, sexp_ratio_to_double(a));
/* ... FALLTHROUGH ... */
#endif
case SEXP_NUM_FLO_CPX:
case SEXP_NUM_FIX_CPX:
case SEXP_NUM_BIG_CPX:
if (!sexp_complexp(a))
a = tmp = sexp_make_complex(ctx, a, SEXP_ZERO);
/* ... FALLTHROUGH ... */
case SEXP_NUM_CPX_CPX:
r = sexp_complex_div(ctx, a, b);
break;
#endif
}
sexp_gc_release1(ctx);
return r;
}
sexp sexp_to_inexact (sexp ctx, sexp a) {
if (sexp_fixnump(a)) return sexp_fixnum_to_flonum(ctx, a);
if (sexp_bignump(a)) return sexp_make_flonum(ctx, sexp_bignum_to_double(a));
return a;
}
sexp sexp_quotient (sexp ctx, sexp a, sexp b) {
int at=sexp_number_type(a), bt=sexp_number_type(b);
sexp r=SEXP_VOID;
sexp_gc_var1(tmp);
if (b == SEXP_ONE) return a;
sexp_gc_preserve1(ctx, tmp);
switch ((at * SEXP_NUM_NUMBER_TYPES) + bt) {
case SEXP_NUM_NOT_NOT: case SEXP_NUM_NOT_FIX:
case SEXP_NUM_NOT_FLO: case SEXP_NUM_NOT_BIG:
r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, a);
break;
case SEXP_NUM_FIX_NOT: case SEXP_NUM_FLO_NOT: case SEXP_NUM_BIG_NOT:
r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, b);
break;
case SEXP_NUM_FLO_FIX: case SEXP_NUM_FLO_FLO: case SEXP_NUM_FLO_BIG:
#if SEXP_USE_RATIOS
case SEXP_NUM_FLO_RAT:
#endif
if (sexp_flonum_value(a) != trunc(sexp_flonum_value(a))) {
r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, a);
} else {
tmp = sexp_bignum_normalize(sexp_double_to_bignum(ctx, sexp_flonum_value(a)));
tmp = sexp_quotient(ctx, tmp, b);
r = sexp_to_inexact(ctx, tmp);
}
break;
#if SEXP_USE_RATIOS
case SEXP_NUM_RAT_FIX: case SEXP_NUM_RAT_BIG: case SEXP_NUM_RAT_RAT:
#endif
#if SEXP_USE_COMPLEX
case SEXP_NUM_FLO_CPX: case SEXP_NUM_CPX_FIX: case SEXP_NUM_CPX_FLO:
case SEXP_NUM_CPX_BIG: case SEXP_NUM_CPX_CPX:
#if SEXP_USE_RATIOS
case SEXP_NUM_CPX_RAT:
#endif
#endif
r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, a);
break;
case SEXP_NUM_FIX_FLO: case SEXP_NUM_BIG_FLO:
#if SEXP_USE_RATIOS
case SEXP_NUM_RAT_FLO:
#endif
if (sexp_flonum_value(b) != trunc(sexp_flonum_value(b))) {
r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, b);
} else {
tmp = sexp_bignum_normalize(sexp_double_to_bignum(ctx, sexp_flonum_value(b)));
tmp = sexp_quotient(ctx, a, tmp);
r = sexp_to_inexact(ctx, tmp);
}
break;
#if SEXP_USE_RATIOS
case SEXP_NUM_FIX_RAT: case SEXP_NUM_BIG_RAT:
#endif
#if SEXP_USE_COMPLEX
case SEXP_NUM_FIX_CPX: case SEXP_NUM_BIG_CPX:
#endif
r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, b);
break;
case SEXP_NUM_FIX_FIX:
r = sexp_fx_div(a, b);
break;
case SEXP_NUM_FIX_BIG:
r = SEXP_ZERO;
break;
case SEXP_NUM_BIG_FIX:
b = tmp = sexp_fixnum_to_bignum(ctx, b);
/* ... FALLTHROUGH ... */
case SEXP_NUM_BIG_BIG:
r = sexp_bignum_normalize(sexp_bignum_quotient(ctx, a, b));
break;
}
sexp_gc_release1(ctx);
return r;
}
sexp sexp_remainder (sexp ctx, sexp a, sexp b) {
int at=sexp_number_type(a), bt=sexp_number_type(b);
sexp r=SEXP_VOID;
sexp_gc_var1(tmp);
if (b == SEXP_ONE) return SEXP_ZERO;
sexp_gc_preserve1(ctx, tmp);
switch ((at * SEXP_NUM_NUMBER_TYPES) + bt) {
case SEXP_NUM_NOT_NOT: case SEXP_NUM_NOT_FIX:
case SEXP_NUM_NOT_FLO: case SEXP_NUM_NOT_BIG:
r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, a);
break;
case SEXP_NUM_FIX_NOT: case SEXP_NUM_FLO_NOT: case SEXP_NUM_BIG_NOT:
r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, b);
break;
case SEXP_NUM_FLO_FIX: case SEXP_NUM_FLO_FLO: case SEXP_NUM_FLO_BIG:
#if SEXP_USE_RATIOS
case SEXP_NUM_FLO_RAT:
#endif
if (isinf(sexp_flonum_value(a)) ||
sexp_flonum_value(a) != trunc(sexp_flonum_value(a))) {
r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, a);
} else if (bt == SEXP_NUM_FLO && isinf(sexp_flonum_value(b))) {
r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, b);
} else {
tmp = sexp_bignum_normalize(sexp_double_to_bignum(ctx, sexp_flonum_value(a)));
tmp = sexp_remainder(ctx, tmp, b);
r = sexp_to_inexact(ctx, tmp);
}
break;
#if SEXP_USE_RATIOS
case SEXP_NUM_RAT_FIX: case SEXP_NUM_RAT_BIG: case SEXP_NUM_RAT_RAT:
#endif
#if SEXP_USE_COMPLEX
case SEXP_NUM_FLO_CPX: case SEXP_NUM_CPX_FIX: case SEXP_NUM_CPX_FLO:
case SEXP_NUM_CPX_BIG: case SEXP_NUM_CPX_CPX:
#if SEXP_USE_RATIOS
case SEXP_NUM_CPX_RAT:
#endif
#endif
r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, a);
break;
case SEXP_NUM_FIX_FLO: case SEXP_NUM_BIG_FLO:
#if SEXP_USE_RATIOS
case SEXP_NUM_RAT_FLO:
#endif
if (isinf(sexp_flonum_value(b)) ||
sexp_flonum_value(b) != trunc(sexp_flonum_value(b))) {
r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, b);
} else {
tmp = sexp_bignum_normalize(sexp_double_to_bignum(ctx, sexp_flonum_value(b)));
tmp = sexp_remainder(ctx, a, tmp);
r = sexp_to_inexact(ctx, tmp);
}
break;
#if SEXP_USE_RATIOS
case SEXP_NUM_FIX_RAT: case SEXP_NUM_BIG_RAT:
#endif
#if SEXP_USE_COMPLEX
case SEXP_NUM_FIX_CPX: case SEXP_NUM_BIG_CPX:
#endif
r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, b);
break;
case SEXP_NUM_FIX_FIX:
r = sexp_fx_rem(a, b);
break;
case SEXP_NUM_FIX_BIG:
r = a;
break;
case SEXP_NUM_BIG_FIX:
r = sexp_bignum_fxrem(ctx, a, sexp_unbox_fixnum(b));
break;
case SEXP_NUM_BIG_BIG:
r = sexp_bignum_normalize(sexp_bignum_remainder(ctx, a, b));
break;
}
sexp_gc_release1(ctx);
return r;
}
sexp sexp_compare (sexp ctx, sexp a, sexp b) {
int at=sexp_number_type(a), bt=sexp_number_type(b);
sexp r=SEXP_VOID;
double f, g;
sexp_gc_var1(tmp);
sexp_gc_preserve1(ctx, tmp);
if (at > bt) {
r = sexp_compare(ctx, b, a);
sexp_negate(r);
} else {
switch ((at * SEXP_NUM_NUMBER_TYPES) + bt) {
case SEXP_NUM_NOT_NOT: case SEXP_NUM_NOT_FIX:
case SEXP_NUM_NOT_FLO: case SEXP_NUM_NOT_BIG:
#if SEXP_USE_COMPLEX
case SEXP_NUM_CPX_CPX: case SEXP_NUM_CPX_FIX:
case SEXP_NUM_CPX_FLO: case SEXP_NUM_CPX_BIG:
#if SEXP_USE_RATIOS
case SEXP_NUM_CPX_RAT:
#endif
#endif
r = sexp_type_exception(ctx, NULL, SEXP_NUMBER, a);
break;
case SEXP_NUM_FIX_FIX:
r = sexp_make_fixnum(sexp_unbox_fixnum(a) - sexp_unbox_fixnum(b));
break;
case SEXP_NUM_FIX_FLO:
f = sexp_fixnum_to_double(a);
g = sexp_flonum_value(b);
r = sexp_make_fixnum(f < g ? -1 : f == g ? 0 : 1);
break;
case SEXP_NUM_FIX_BIG:
if ((sexp_bignum_hi(b) > 1) ||
(sexp_bignum_data(b)[0] > SEXP_MAX_FIXNUM))
r = sexp_make_fixnum(sexp_bignum_sign(b) < 0 ? 1 : -1);
else
r = sexp_make_fixnum(sexp_unbox_fixnum(a) - (sexp_sint_t)sexp_bignum_data(b)[0]);
break;
case SEXP_NUM_FLO_FLO:
f = sexp_flonum_value(a);
g = sexp_flonum_value(b);
r = sexp_make_fixnum(f < g ? -1 : f == g ? 0 : 1);
break;
case SEXP_NUM_FLO_BIG:
f = sexp_flonum_value(a);
if (isinf(f)) {
r = f > 0 ? SEXP_ONE : SEXP_NEG_ONE;
break;
} else if (isnan(f)) {
r = sexp_xtype_exception(ctx, NULL, "can't compare NaN", a);
break;
} else {
a = tmp = sexp_double_to_bignum(ctx, f);
}
/* ... FALLTHROUGH ... */
case SEXP_NUM_BIG_BIG:
r = sexp_make_fixnum(sexp_bignum_compare(a, b));
break;
#if SEXP_USE_RATIOS
case SEXP_NUM_FLO_RAT:
f = sexp_flonum_value(a);
if (isinf(f)) {
r = f > 0 ? SEXP_ONE : SEXP_NEG_ONE;
} else if (isnan(f)) {
r = sexp_xtype_exception(ctx, NULL, "can't compare NaN", a);
} else {
g = sexp_ratio_to_double(b);
r = sexp_make_fixnum(f < g ? -1 : f == g ? 0 : 1);
}
break;
case SEXP_NUM_FIX_RAT:
case SEXP_NUM_BIG_RAT:
a = tmp = sexp_make_ratio(ctx, a, SEXP_ONE);
/* ... FALLTHROUGH ... */
case SEXP_NUM_RAT_RAT:
r = sexp_ratio_compare(ctx, a, b);
break;
#endif
}
}
sexp_gc_release1(ctx);
return r;
}
#endif
Reference in a new issue View git blame Copy permalink