1 /*  bignum.c -- bignum support                                */
2 /*  Copyright (c) 2009-2015 Alex Shinn.  All rights reserved. */
3 /*  BSD-style license: http://synthcode.com/license.txt       */
4 
5 #include "chibi/sexp.h"
6 
7 #if SEXP_USE_BIGNUMS
8 
9 #define SEXP_INIT_BIGNUM_SIZE 2
10 
digit_value(int c)11 static int digit_value (int c) {
12   return (((c)<='9') ? ((c) - '0') : ((sexp_toupper(c) - 'A') + 10));
13 }
14 
hex_digit(int n)15 static int hex_digit (int n) {
16   return ((n<=9) ? ('0' + n) : ('A' + n - 10));
17 }
18 
sexp_make_bignum(sexp ctx,sexp_uint_t len)19 sexp sexp_make_bignum (sexp ctx, sexp_uint_t len) {
20   sexp_uint_t size = sexp_sizeof(bignum) + len*sizeof(sexp_uint_t);
21   sexp res = sexp_alloc_tagged(ctx, size, SEXP_BIGNUM);
22   if (!sexp_exceptionp(res)) {
23     sexp_bignum_length(res) = len;
24     sexp_bignum_sign(res) = 1;
25   }
26   return res;
27 }
28 
sexp_fixnum_to_bignum(sexp ctx,sexp a)29 sexp sexp_fixnum_to_bignum (sexp ctx, sexp a) {
30   sexp res = sexp_make_bignum(ctx, 1);
31   if (!sexp_exceptionp(res)) {
32     sexp_bignum_data(res)[0] = sexp_unbox_fx_abs(a);
33     sexp_bignum_sign(res) = sexp_fx_sign(a);
34   }
35   return res;
36 }
37 
sexp_make_integer_from_lsint(sexp ctx,sexp_lsint_t x)38 sexp sexp_make_integer_from_lsint (sexp ctx, sexp_lsint_t x) {
39   sexp res;
40   if (lsint_is_fixnum(x)) {
41     res = sexp_make_fixnum(lsint_to_sint(x));
42   } else if (sexp_lsint_fits_sint(x)) {
43     res = sexp_make_bignum(ctx, 1);
44     if (lsint_lt_0(x)) {
45       sexp_bignum_sign(res) = -1;
46       sexp_bignum_data(res)[0] = (sexp_uint_t)-lsint_to_sint(x);
47     } else {
48       sexp_bignum_sign(res) = 1;
49       sexp_bignum_data(res)[0] = (sexp_uint_t)lsint_to_sint(x);
50     }
51   } else {
52     res = sexp_make_bignum(ctx, 2);
53     if (lsint_lt_0(x)) {
54       sexp_bignum_sign(res) = -1;
55       sexp_bignum_data(res)[0] = (sexp_uint_t)-lsint_to_sint(x);
56       sexp_bignum_data(res)[1] = (sexp_uint_t)~lsint_to_sint_hi(x);
57     } else {
58       sexp_bignum_sign(res) = 1;
59       sexp_bignum_data(res)[0] = (sexp_uint_t)lsint_to_sint(x);
60       sexp_bignum_data(res)[1] = (sexp_uint_t)lsint_to_sint_hi(x);
61     }
62   }
63   return res;
64 }
65 
sexp_make_unsigned_integer_from_luint(sexp ctx,sexp_luint_t x)66 sexp sexp_make_unsigned_integer_from_luint (sexp ctx, sexp_luint_t x) {
67   sexp res;
68   if (luint_is_fixnum(x)) {
69     res = sexp_make_fixnum(luint_to_uint(x));
70   } else if (sexp_luint_fits_uint(x)) {
71     res = sexp_make_bignum(ctx, 1);
72     sexp_bignum_sign(res) = 1;
73     sexp_bignum_data(res)[0] = luint_to_uint(x);
74   } else {
75     res = sexp_make_bignum(ctx, 2);
76     sexp_bignum_sign(res) = 1;
77     sexp_bignum_data(res)[0] = luint_to_uint(x);
78     sexp_bignum_data(res)[1] = luint_to_uint_hi(x);
79   }
80   return res;
81 }
82 
83 #if SEXP_USE_CUSTOM_LONG_LONGS
sexp_make_integer(sexp ctx,long long x)84 sexp sexp_make_integer(sexp ctx, long long x) {
85   return sexp_make_integer_from_lsint(ctx, lsint_from_sint(x));
86 }
sexp_make_unsigned_integer(sexp ctx,unsigned long long x)87 sexp sexp_make_unsigned_integer(sexp ctx, unsigned long long x) {
88   return sexp_make_unsigned_integer_from_luint(ctx, luint_from_uint(x));
89 }
90 #else
sexp_make_integer(sexp ctx,sexp_lsint_t x)91 sexp sexp_make_integer (sexp ctx, sexp_lsint_t x) {
92   return sexp_make_integer_from_lsint(ctx, x);
93 }
sexp_make_unsigned_integer(sexp ctx,sexp_luint_t x)94 sexp sexp_make_unsigned_integer (sexp ctx, sexp_luint_t x) {
95   return sexp_make_unsigned_integer_from_luint(ctx, x);
96 }
97 #endif
98 
99 #if !SEXP_64_BIT
sexp_bignum_to_sint(sexp x)100 long long sexp_bignum_to_sint(sexp x) {
101   if (!sexp_bignump(x))
102     return 0;
103   if (sexp_bignum_length(x) > 1)
104     return sexp_bignum_sign(x) * (
105       (((long long)(sexp_bignum_data(x)[1]))<<(8*sizeof(sexp_bignum_data(x)[0]))) + sexp_bignum_data(x)[0]);
106   return sexp_bignum_sign(x) * sexp_bignum_data(x)[0];
107 }
108 
sexp_bignum_to_uint(sexp x)109 unsigned long long sexp_bignum_to_uint(sexp x) {
110   if (!sexp_bignump(x))
111     return 0;
112   if (sexp_bignum_length(x) > 1)
113     return (((unsigned long long)(sexp_bignum_data(x)[1]))<<(8*sizeof(sexp_bignum_data(x)[0]))) + sexp_bignum_data(x)[0];
114   return sexp_bignum_data(x)[0];
115 }
116 #endif
117 
118 #define double_trunc_10s_digit(f) (trunc((f)/10.0)*10.0)
119 #define double_10s_digit(f) ((f)-double_trunc_10s_digit(f))
120 
121 #define double_16s_digit(f) fmod(f,16.0)
122 
sexp_double_to_bignum(sexp ctx,double f)123 sexp sexp_double_to_bignum (sexp ctx, double f) {
124   int sign;
125   sexp_gc_var3(res, scale, tmp);
126   sexp_gc_preserve3(ctx, res, scale, tmp);
127   res = sexp_fixnum_to_bignum(ctx, SEXP_ZERO);
128   scale = sexp_fixnum_to_bignum(ctx, SEXP_ONE);
129   sign = (f < 0 ? -1 : 1);
130   for (f=fabs(f); f >= 1.0; f=trunc(f/16)) {
131     tmp = sexp_bignum_fxmul(ctx, NULL, scale, (sexp_uint_t)double_16s_digit(f), 0);
132     res = sexp_bignum_add(ctx, res, res, tmp);
133     scale = sexp_bignum_fxmul(ctx, NULL, scale, 16, 0);
134   }
135   sexp_bignum_sign(res) = sign;
136   sexp_gc_release3(ctx);
137   return res;
138 }
139 
sexp_copy_bignum(sexp ctx,sexp dst,sexp a,sexp_uint_t len0)140 sexp sexp_copy_bignum (sexp ctx, sexp dst, sexp a, sexp_uint_t len0) {
141   sexp_uint_t len = (len0 > 0) ? len0 : sexp_bignum_length(a), size;
142   size = sexp_sizeof(bignum) + len*sizeof(sexp_uint_t);
143   if (! dst || sexp_bignum_length(dst) < len) {
144     dst = sexp_alloc_tagged(ctx, size, SEXP_BIGNUM);
145     sexp_bignum_length(dst) = len;
146   }
147   if (sexp_bignum_length(a) < len)
148     len = sexp_bignum_length(a);
149 
150   sexp_bignum_sign(dst) = sexp_bignum_sign(a);
151   memset(sexp_bignum_data(dst), 0,
152          sexp_bignum_length(dst)*sizeof(sexp_uint_t));
153   memmove(sexp_bignum_data(dst), sexp_bignum_data(a),
154           len*sizeof(sexp_uint_t));
155   return dst;
156 }
157 
sexp_bignum_zerop(sexp a)158 int sexp_bignum_zerop (sexp a) {
159   int i;
160   sexp_uint_t *data = sexp_bignum_data(a);
161   for (i=sexp_bignum_length(a)-1; i>=0; i--)
162     if (data[i])
163       return 0;
164   return 1;
165 }
166 
sexp_bignum_hi(sexp a)167 sexp_uint_t sexp_bignum_hi (sexp a) {
168   sexp_uint_t i=sexp_bignum_length(a)-1;
169   while ((i>0) && ! sexp_bignum_data(a)[i])
170     i--;
171   return i+1;
172 }
173 
sexp_bignum_compare_abs(sexp a,sexp b)174 sexp_sint_t sexp_bignum_compare_abs (sexp a, sexp b) {
175   int ai=sexp_bignum_hi(a), bi=sexp_bignum_hi(b);
176   sexp_uint_t *adata=sexp_bignum_data(a), *bdata=sexp_bignum_data(b);
177   if (ai != bi)
178     return ai - bi;
179   for (--ai; ai >= 0; ai--) {
180     if (adata[ai] > bdata[ai])
181       return 1;
182     else if (adata[ai] < bdata[ai])
183       return -1;
184   }
185   return 0;
186 }
187 
sexp_bignum_compare(sexp a,sexp b)188 sexp_sint_t sexp_bignum_compare (sexp a, sexp b) {
189   if (sexp_bignum_sign(a) != sexp_bignum_sign(b))
190     return sexp_bignum_sign(a);
191   sexp_sint_t cmp = sexp_bignum_compare_abs(a, b);
192   return sexp_bignum_sign(a) < 0 ? -cmp : cmp;
193 }
194 
sexp_bignum_normalize(sexp a)195 sexp sexp_bignum_normalize (sexp a) {
196   sexp_uint_t *data;
197   if ((! sexp_bignump(a)) || (sexp_bignum_hi(a)>1))
198     return a;
199   data = sexp_bignum_data(a);
200   if ((data[0] > SEXP_MAX_FIXNUM)
201       && ! ((sexp_bignum_sign(a) == -1) && (data[0] == SEXP_MAX_FIXNUM+1)))
202     return a;
203   return sexp_make_fixnum((sexp_sint_t)data[0] * sexp_bignum_sign(a));
204 }
205 
sexp_bignum_to_double(sexp a)206 double sexp_bignum_to_double (sexp a) {
207   double res = 0;
208   sexp_sint_t i;
209   sexp_uint_t *data=sexp_bignum_data(a);
210   for (i=sexp_bignum_hi(a)-1; i>=0; i--)
211     res = res * ((double)SEXP_UINT_T_MAX+1) + data[i];
212   return res * sexp_bignum_sign(a);
213 }
214 
sexp_bignum_fxadd(sexp ctx,sexp a,sexp_uint_t b)215 sexp sexp_bignum_fxadd (sexp ctx, sexp a, sexp_uint_t b) {
216   sexp_uint_t len=sexp_bignum_hi(a), *data=sexp_bignum_data(a),
217     carry=b, i=0, n;
218   do { n = data[i];
219        data[i] += carry;
220        carry = (n > (SEXP_UINT_T_MAX - carry));
221   } while (++i<len && carry);
222   if (carry) {
223     a = sexp_copy_bignum(ctx, NULL, a, len+1);
224     sexp_bignum_data(a)[len] = 1;
225   }
226   return a;
227 }
228 
sexp_bignum_fxsub(sexp ctx,sexp a,sexp_uint_t b)229 sexp sexp_bignum_fxsub (sexp ctx, sexp a, sexp_uint_t b) {
230   sexp_uint_t len=sexp_bignum_hi(a), *data=sexp_bignum_data(a), borrow, i=0, n;
231   if ((len == 1) && (b > data[0])) {
232     data[0] = b - data[0];
233     sexp_bignum_sign(a) = -sexp_bignum_sign(a);
234   } else {
235     for (borrow=b; borrow; i++) {
236       n = data[i];
237       data[i] -= borrow;
238       borrow = (n < borrow);
239     }
240   }
241   return a;
242 }
243 
sexp_bignum_fxmul(sexp ctx,sexp d,sexp a,sexp_uint_t b,int offset)244 sexp sexp_bignum_fxmul (sexp ctx, sexp d, sexp a, sexp_uint_t b, int offset) {
245   sexp_uint_t len=sexp_bignum_length(a), *data, *adata=sexp_bignum_data(a),
246     carry=0, i;
247   sexp_luint_t n;
248   sexp_gc_var1(tmp);
249   sexp_gc_preserve1(ctx, tmp);
250   if ((! d) || (sexp_bignum_length(d) < len+offset))
251     d = sexp_make_bignum(ctx, len+offset);
252   tmp = d;
253   data = sexp_bignum_data(d);
254   for (i=0; i<len; i++) {
255     n = luint_add(luint_mul_uint(luint_from_uint(adata[i]), b), luint_from_uint(carry));
256     data[i+offset] = luint_to_uint(n);
257     carry = luint_to_uint(luint_shr(n, (sizeof(sexp_uint_t)*8)));
258   }
259   if (carry) {
260     if (sexp_bignum_length(d) <= len+offset)
261       d = sexp_copy_bignum(ctx, NULL, d, len+offset+1);
262     sexp_bignum_data(d)[len+offset] = carry;
263   }
264   sexp_gc_release1(ctx);
265   return d;
266 }
267 
sexp_bignum_fxdiv(sexp ctx,sexp a,sexp_uint_t b,int offset)268 sexp_uint_t sexp_bignum_fxdiv (sexp ctx, sexp a, sexp_uint_t b, int offset) {
269   sexp_uint_t len=sexp_bignum_hi(a), *data=sexp_bignum_data(a), q, r=0;
270   int i;
271   sexp_luint_t n = luint_from_uint(0);
272   for (i=len-1; i>=offset; i--) {
273     n = luint_add(luint_shl(n,  sizeof(sexp_uint_t)*8), luint_from_uint(data[i]));
274     q = luint_to_uint(luint_div_uint(n, b));
275     r = luint_to_uint(luint_sub(n, luint_mul_uint(luint_from_uint(q),  b)));
276     data[i] = q;
277     n = luint_from_uint(r);
278   }
279   return r;
280 }
281 
sexp_bignum_fxrem(sexp ctx,sexp a,sexp_sint_t b)282 sexp sexp_bignum_fxrem (sexp ctx, sexp a, sexp_sint_t b) {
283   sexp_uint_t len=sexp_bignum_hi(a), *data=sexp_bignum_data(a), q, b0;
284   int i;
285   sexp_luint_t n = luint_from_uint(0);
286   if (b > 0) {
287     q = b - 1;
288     if ((b & q) == 0)
289       return sexp_make_fixnum(sexp_bignum_sign(a) * (data[0] & q));
290   }
291   b0 = (b >= 0) ? b : -b;
292   if (b0 == 0) {
293     return sexp_xtype_exception(ctx, NULL, "divide by zero", a);
294   }
295   for (i=len-1; i>=0; i--) {
296     n = luint_add(luint_shl(n, sizeof(sexp_uint_t)*8), luint_from_uint(data[i]));
297     q = luint_to_uint(luint_div_uint(n, b0));
298     n = luint_sub(n, luint_mul_uint(luint_from_uint(q), b0));
299   }
300   return sexp_make_fixnum(sexp_bignum_sign(a) * (sexp_sint_t)luint_to_uint(n));
301 }
302 
sexp_read_bignum(sexp ctx,sexp in,sexp_uint_t init,signed char sign,sexp_uint_t base)303 sexp sexp_read_bignum (sexp ctx, sexp in, sexp_uint_t init,
304                        signed char sign, sexp_uint_t base) {
305   int c, digit;
306   sexp_gc_var3(res, tmp, imag);
307   sexp_gc_preserve3(ctx, res, tmp, imag);
308   res = sexp_make_bignum(ctx, SEXP_INIT_BIGNUM_SIZE);
309   sexp_bignum_sign(res) = sign;
310   sexp_bignum_data(res)[0] = init;
311   for (c=sexp_read_char(ctx, in); sexp_isxdigit(c); c=sexp_read_char(ctx, in)) {
312     digit = digit_value(c);
313     if ((digit < 0) || (digit >= (int)base))
314       break;
315     res = sexp_bignum_fxmul(ctx, res, res, base, 0);
316     res = sexp_bignum_fxadd(ctx, res, digit);
317   }
318   if (c=='.' || c=='e' || c=='E') {
319     if (base != 10) {
320       res = sexp_read_error(ctx, "found non-base 10 float", SEXP_NULL, in);
321     } else if (c=='.') {
322       res = sexp_read_float_tail(ctx, in, sexp_bignum_to_double(res), (sign==-1));
323     } else {
324       tmp = sexp_read_number(ctx, in, base, 0);
325 #if SEXP_USE_COMPLEX
326       if (sexp_complexp(tmp)) {
327         imag = sexp_complex_imag(tmp);
328         tmp = sexp_complex_real(tmp);
329       } else {
330         imag = SEXP_ZERO;
331       }
332 #endif
333       if (sexp_exceptionp(tmp)) {
334         res = tmp;
335       } else if (sexp_fixnump(tmp) && labs(sexp_unbox_fixnum(tmp)) < 100*1024*1024) {
336         tmp = sexp_expt(ctx, SEXP_TEN, tmp);
337         res = sexp_mul(ctx, res, tmp);
338       } else {
339         tmp = sexp_exact_to_inexact(ctx, NULL, 2, tmp);
340         tmp = sexp_expt(ctx, SEXP_TEN, tmp);
341         res = sexp_mul(ctx, res, tmp);
342       }
343 #if SEXP_USE_COMPLEX
344       if (imag != SEXP_ZERO && !sexp_exceptionp(res))
345         res = sexp_make_complex(ctx, res, imag);
346 #endif
347     }
348 #if SEXP_USE_RATIOS
349   } else if (c=='/') {
350     res = sexp_bignum_normalize(res);
351     res = sexp_make_ratio(ctx, res, SEXP_ONE);
352     sexp_ratio_denominator(res) = sexp_read_number(ctx, in, 10, 0);
353     res = sexp_ratio_normalize(ctx, res, in);
354 #endif
355 #if SEXP_USE_COMPLEX
356   } else if (c=='i' || c=='i' || c=='+' || c=='-') {
357     sexp_push_char(ctx, c, in);
358     res = sexp_bignum_normalize(res);
359     res = sexp_read_complex_tail(ctx, in, res);
360 #endif
361   } else if ((c!=EOF) && ! sexp_is_separator(c)) {
362     res = sexp_read_error(ctx, "invalid numeric syntax",
363                           sexp_make_character(c), in);
364   } else {
365     sexp_push_char(ctx, c, in);
366   }
367   sexp_gc_release3(ctx);
368   return sexp_bignum_normalize(res);
369 }
370 
log2i(int v)371 static int log2i(int v) {
372   int i;
373   for (i = 0; i < sizeof(v)*8; i++)
374     if ((1<<(i+1)) > v)
375       break;
376   return i;
377 }
378 
sexp_write_bignum(sexp ctx,sexp a,sexp out,sexp_uint_t base)379 sexp sexp_write_bignum (sexp ctx, sexp a, sexp out, sexp_uint_t base) {
380   int i, str_len, lg_base = log2i(base);
381   char *data;
382   sexp_gc_var2(b, str);
383   sexp_gc_preserve2(ctx, b, str);
384   b = sexp_copy_bignum(ctx, NULL, a, 0);
385   sexp_bignum_sign(b) = 1;
386   if (lg_base < 1) {
387     return sexp_xtype_exception(ctx, NULL, "number base too small", a);
388   }
389   i = str_len = (sexp_bignum_length(b)*sizeof(sexp_uint_t)*8 + lg_base - 1)
390     / lg_base + 1;
391   str = sexp_make_string(ctx, sexp_make_fixnum(str_len),
392                          sexp_make_character(' '));
393   data = sexp_string_data(str);
394   while (! sexp_bignum_zerop(b))
395     data[--i] = hex_digit(sexp_bignum_fxdiv(ctx, b, base, 0));
396   if (i == str_len)
397     data[--i] = '0';
398   else if (sexp_bignum_sign(a) == -1)
399     data[--i] = '-';
400   sexp_write_string(ctx, data + i, out);
401   sexp_gc_release2(ctx);
402   return SEXP_VOID;
403 }
404 
405 /****************** bignum arithmetic *************************/
406 
sexp_bignum_add_fixnum(sexp ctx,sexp a,sexp b)407 sexp sexp_bignum_add_fixnum (sexp ctx, sexp a, sexp b) {
408   sexp_gc_var1(c);
409   sexp_gc_preserve1(ctx, c);
410   c = sexp_copy_bignum(ctx, NULL, a, 0);
411   if (sexp_bignum_sign(c) == sexp_fx_sign(b))
412     c = sexp_bignum_fxadd(ctx, c, sexp_unbox_fx_abs(b));
413   else
414     c = sexp_bignum_fxsub(ctx, c, sexp_unbox_fx_abs(b));
415   sexp_gc_release1(ctx);
416   return c;
417 }
418 
sexp_bignum_sub_digits(sexp ctx,sexp dst,sexp a,sexp b)419 sexp sexp_bignum_sub_digits (sexp ctx, sexp dst, sexp a, sexp b) {
420   sexp_uint_t alen=sexp_bignum_hi(a), blen=sexp_bignum_hi(b),
421     borrow=0, i, *adata, *bdata, *cdata;
422   sexp_gc_var1(c);
423   if ((alen < blen) || ((alen == blen) && (sexp_bignum_compare_abs(a, b) < 0)))
424     return sexp_bignum_sub_digits(ctx, dst, b, a);
425   sexp_gc_preserve1(ctx, c);
426   c = ((dst && sexp_bignum_hi(dst) >= alen)
427        ? dst : sexp_copy_bignum(ctx, NULL, a, 0));
428   adata = sexp_bignum_data(a);
429   bdata = sexp_bignum_data(b);
430   cdata = sexp_bignum_data(c);
431   for (i=0; i<blen; i++) {
432     if (adata[i] > bdata[i] || (adata[i] == bdata[i] && !borrow)) {
433       cdata[i] = adata[i] - bdata[i] - borrow;
434       borrow = 0;
435     } else {
436       cdata[i] = (SEXP_UINT_T_MAX - bdata[i]);
437       cdata[i] += 1;
438       cdata[i] -= borrow;
439       cdata[i] += adata[i];
440       borrow = 1;
441     }
442   }
443   for ( ; borrow && (i<alen); i++) {
444     borrow = (cdata[i] == 0 ? 1 : 0);
445     cdata[i]--;
446   }
447   sexp_gc_release1(ctx);
448   return c;
449 }
450 
sexp_bignum_add_digits(sexp ctx,sexp dst,sexp a,sexp b)451 sexp sexp_bignum_add_digits (sexp ctx, sexp dst, sexp a, sexp b) {
452   sexp_uint_t alen=sexp_bignum_hi(a), blen=sexp_bignum_hi(b),
453     carry=0, i, old_a, p_sum, *adata, *bdata, *cdata;
454   sexp_gc_var1(c);
455   if (alen < blen) return sexp_bignum_add_digits(ctx, dst, b, a);
456   sexp_gc_preserve1(ctx, c);
457   c = ((dst && sexp_bignum_hi(dst) >= alen)
458        ? dst : sexp_copy_bignum(ctx, NULL, a, 0));
459   adata = sexp_bignum_data(a);
460   bdata = sexp_bignum_data(b);
461   cdata = sexp_bignum_data(c);
462   for (i=0; i<blen; i++) {
463     old_a = adata[i]; /* adata may alias cdata */
464     p_sum = adata[i] + bdata[i];
465     cdata[i] = p_sum + carry;
466     carry = (old_a > (SEXP_UINT_T_MAX - bdata[i]) ? 1 : 0)
467           + (p_sum > (SEXP_UINT_T_MAX - carry) ? 1 : 0);
468   }
469   for ( ; carry && (i<alen); i++) {
470     carry = (cdata[i] == SEXP_UINT_T_MAX ? 1 : 0);
471     cdata[i]++;
472   }
473   if (carry) {
474     c = sexp_copy_bignum(ctx, NULL, c, alen+1);
475     sexp_bignum_data(c)[alen] = 1;
476   }
477   sexp_gc_release1(ctx);
478   return c;
479 }
480 
sexp_bignum_add(sexp ctx,sexp dst,sexp a,sexp b)481 sexp sexp_bignum_add (sexp ctx, sexp dst, sexp a, sexp b) {
482   sexp res;
483   if (sexp_bignum_sign(a) == sexp_bignum_sign(b)) {
484     res = sexp_bignum_add_digits(ctx, dst, a, b);
485     sexp_bignum_sign(res) = sexp_bignum_sign(a);
486   } else {
487     res = sexp_bignum_sub_digits(ctx, dst, a, b);
488     sexp_bignum_sign(res)
489       = sexp_bignum_sign(sexp_bignum_compare_abs(a, b) >= 0 ? a : b);
490   }
491   return res;
492 }
493 
sexp_bignum_sub(sexp ctx,sexp dst,sexp a,sexp b)494 sexp sexp_bignum_sub (sexp ctx, sexp dst, sexp a, sexp b) {
495   sexp res;
496   int sign;
497   if (sexp_bignum_sign(a) == sexp_bignum_sign(b)) {
498     sign = (sexp_bignum_compare_abs(a, b) >= 0 ? sexp_bignum_sign(a)
499             : -sexp_bignum_sign(a));
500     res = sexp_bignum_sub_digits(ctx, dst, a, b);
501   } else {
502     sign = sexp_bignum_sign(a);
503     res = sexp_bignum_add_digits(ctx, dst, a, b);
504   }
505   sexp_bignum_sign(res) = sign;
506   return res;
507 }
508 
sexp_bignum_split(sexp ctx,sexp a,sexp_uint_t k,sexp * lo,sexp * hi)509 static void sexp_bignum_split (sexp ctx, sexp a, sexp_uint_t k, sexp* lo, sexp* hi) {
510   sexp_uint_t alen=sexp_bignum_hi(a), i, *adata=sexp_bignum_data(a),
511     *lodata, *hidata;
512   *lo = sexp_make_bignum(ctx, k);        /* must be gc protected by caller */
513   *hi = sexp_make_bignum(ctx, alen-k+1);
514   lodata = sexp_bignum_data(*lo);
515   hidata = sexp_bignum_data(*hi);
516   for (i=0; i<k; i++)                    /* split into a[0..k-1], a[k..] */
517     lodata[i] = adata[i];
518   for (i=k; i<alen; i++)
519     hidata[i-k] = adata[i];
520 }
521 
sexp_bignum_shift(sexp ctx,sexp a,sexp_uint_t k)522 static sexp sexp_bignum_shift (sexp ctx, sexp a, sexp_uint_t k) {
523   sexp res;
524   sexp_uint_t alen = sexp_bignum_hi(a), i;
525   res = sexp_make_bignum(ctx, alen + k + 1);
526   for (i=0; i<alen; i++)
527     sexp_bignum_data(res)[i+k] = sexp_bignum_data(a)[i];
528   return res;
529 }
530 
sexp_bignum_mul(sexp ctx,sexp dst,sexp a,sexp b)531 sexp sexp_bignum_mul (sexp ctx, sexp dst, sexp a, sexp b) {
532   sexp_uint_t alen=sexp_bignum_hi(a), blen=sexp_bignum_hi(b), k,
533     *bdata=sexp_bignum_data(b);
534   sexp_gc_var7(a0, a1, b0, b1, z0, z1, z2);
535   if (alen < blen) return sexp_bignum_mul(ctx, dst, b, a);
536   if (blen == 1) {
537     z1 = sexp_bignum_fxmul(ctx, dst, a, bdata[0], 0);
538   } else {
539     /* karatsuba:                             */
540     /*   ab = (a1B^k + a0) * (b1B^k + b0)     */
541     /*      = z2B^2k + z1B^k + z0             */
542     /* where:                                 */
543     /*   z2 = a1b1                            */
544     /*   z1 = a1b0 + a0b1                     */
545     /*   z0 = a0b0                            */
546     /* then optimize further:                 */
547     /*   z1 = (a1 + a0)(b1 + b0) - z2 - z0    */
548     sexp_gc_preserve7(ctx, a0, a1, b0, b1, z0, z1, z2);
549     k = blen / 2;
550     sexp_bignum_split(ctx, a, k, &a0, &a1);
551     sexp_bignum_split(ctx, b, k, &b0, &b1);
552     z0 = sexp_bignum_add(ctx, NULL, a1, a0);  /* temp */
553     z1 = sexp_bignum_add(ctx, NULL, b1, b0);
554     z1 = sexp_bignum_mul(ctx, NULL, z1, z0);  /* 1 */
555     z0 = sexp_bignum_mul(ctx, NULL, a0, b0);  /* 2 */
556     z2 = sexp_bignum_mul(ctx, NULL, a1, b1);  /* 3 */
557     z1 = sexp_bignum_sub(ctx, NULL, z1, z0);
558     z1 = sexp_bignum_sub(ctx, NULL, z1, z2);
559     z2 = sexp_bignum_shift(ctx, z2, 2*k);
560     z1 = sexp_bignum_shift(ctx, z1, k);
561     z1 = sexp_bignum_add(ctx, z1, z1, z0);
562     z1 = sexp_bignum_add(ctx, z1, z1, z2);
563     sexp_gc_release7(ctx);
564   }
565   sexp_bignum_sign(z1) = sexp_bignum_sign(a) * sexp_bignum_sign(b);
566   return z1;
567 }
568 
sexp_bignum_quot_rem(sexp ctx,sexp * rem,sexp a,sexp b)569 sexp sexp_bignum_quot_rem (sexp ctx, sexp *rem, sexp a, sexp b) {
570   sexp_luint_t d, dn, dd;
571   sexp_uint_t dlo, dhi;
572   sexp_sint_t alen, blen=sexp_bignum_hi(b), sign=1, off=0;
573   sexp_gc_var5(q, x, y, a1, b1);
574   if (blen == 1 && sexp_bignum_data(b)[0] == 0)
575     return sexp_xtype_exception(ctx, NULL, "divide by zero", a);
576   sexp_gc_preserve5(ctx, q, x, y, a1, b1);
577   a1 = sexp_copy_bignum(ctx, NULL, a, 0);
578   sexp_bignum_sign(a1) = 1;
579   /* fast path for single bigit divisor */
580   if (blen == 1) {
581     b1 = sexp_make_bignum(ctx, 1);
582     sexp_bignum_data(b1)[0] = sexp_bignum_fxdiv(ctx, a1, sexp_bignum_data(b)[0], 0);
583     *rem = sexp_bignum_normalize(b1);
584     if (sexp_bignum_sign(a) * sexp_bignum_sign(b) < 0) {
585       sexp_negate_exact(a1);
586     }
587     if (sexp_bignum_sign(a) < 0) {
588       sexp_negate_exact(*rem);
589     }
590     sexp_gc_release5(ctx);
591     return a1;
592   }
593   /* general case */
594   b1 = sexp_copy_bignum(ctx, NULL, b, 0);
595   sexp_bignum_sign(b1) = 1;
596   q = SEXP_ZERO;
597   x = sexp_make_bignum(ctx, sexp_bignum_length(a));
598   while (sexp_bignum_compare_abs(a1, b1) >= 0) { /* a1, b1 at least 2 bigits */
599     /* guess divisor x */
600     alen = sexp_bignum_hi(a1);
601     sexp_bignum_data(x)[off] = 0;
602     if (off > 0) sexp_bignum_data(x)[off-1] = 0;
603     off = alen - blen + 1;
604     dn = luint_add_uint(luint_shl(luint_from_uint(sexp_bignum_data(a1)[alen-1])
605            , (sizeof(sexp_uint_t)*8))
606           , sexp_bignum_data(a1)[alen-2]);
607     dd = luint_add_uint(luint_shl(luint_from_uint(sexp_bignum_data(b1)[blen-1])
608            , (sizeof(sexp_uint_t)*8))
609           , sexp_bignum_data(b1)[blen-2]);
610     if (alen > 2 && blen > 2 &&
611         luint_lt(luint_from_uint(sexp_bignum_data(a1)[alen-1]), luint_shl(luint_from_uint(1), (sizeof(sexp_uint_t)*4))) &&
612         luint_lt(luint_from_uint(sexp_bignum_data(b1)[blen-1]), luint_shl(luint_from_uint(1), (sizeof(sexp_uint_t)*4)))) {
613       dn = luint_add_uint(luint_shl(dn, (sizeof(sexp_uint_t)*4))
614         , (sexp_bignum_data(a1)[alen-3] >> (sizeof(sexp_uint_t)*4)));
615       dd = luint_add_uint(luint_shl(dd, (sizeof(sexp_uint_t)*4))
616         , (sexp_bignum_data(b1)[blen-3] >> (sizeof(sexp_uint_t)*4)));
617     }
618     d = luint_div(dn, dd);
619     if (luint_eq(d, luint_from_uint(0))) {
620       dn = luint_add_uint(luint_shl(luint_from_uint(sexp_bignum_data(a1)[alen-1])
621              , (sizeof(sexp_uint_t)*8))
622             , sexp_bignum_data(a1)[alen-2]);
623       dd = luint_from_uint(sexp_bignum_data(b1)[blen-1]);
624       if (luint_lt(luint_from_uint(sexp_bignum_data(a1)[alen-1]), (luint_shl(luint_from_uint(1), (sizeof(sexp_uint_t)*4)))) &&
625           luint_lt(luint_from_uint(sexp_bignum_data(b1)[blen-1]), (luint_shl(luint_from_uint(1), (sizeof(sexp_uint_t)*4))))) {
626         dn = luint_add_uint(luint_shl(dn, (sizeof(sexp_uint_t)*4))
627           , (sexp_bignum_data(a1)[alen-3] >> (sizeof(sexp_uint_t)*4)));
628         dd = luint_add_uint(luint_shl(dd, (sizeof(sexp_uint_t)*4))
629           , (sexp_bignum_data(b1)[blen-2] >> (sizeof(sexp_uint_t)*4)));
630       }
631       d = luint_div(dn, dd);
632       off--;
633     }
634     dhi = luint_to_uint(luint_shr(d, (sizeof(sexp_uint_t)*8)));
635     dlo = luint_to_uint(luint_and(d, luint_sub(luint_shl(luint_from_uint(1), (sizeof(sexp_uint_t)*8)), luint_from_uint(1))));
636     sexp_bignum_data(x)[off] = dhi;
637     if (off > 0) sexp_bignum_data(x)[off-1] = dlo;
638     /* update quotient q and remainder a1 estimates */
639     y = sexp_bignum_mul(ctx, NULL, b1, x);
640     if (sign < 0) {
641       a1 = sexp_bignum_add(ctx, NULL, a1, y);
642       q = sexp_sub(ctx, q, x);
643     } else {
644       a1 = sexp_bignum_sub(ctx, NULL, a1, y);
645       q = sexp_add(ctx, q, x);
646     }
647     /* flip the sign if we overshot in our estimate */
648     if (sexp_bignum_sign(a1) != sign) {
649       sexp_bignum_sign(a1) = (char)(-sign);
650       sign *= -1;
651     }
652   }
653   /* adjust signs */
654   a1 = sexp_bignum_normalize(a1);
655   if (sign < 0 && a1 != SEXP_ZERO) {
656     q = sexp_sub(ctx, q, SEXP_ONE);
657     a1 = sexp_add(ctx, a1, b1);
658   }
659   *rem = a1;
660   if (sexp_bignum_sign(a) * sexp_bignum_sign(b) < 0) {
661     sexp_negate_exact(q);
662   }
663   if (sexp_bignum_sign(a) < 0) {
664     sexp_negate_exact(*rem);
665   }
666   sexp_gc_release5(ctx);
667   return q;
668 }
669 
sexp_bignum_quotient(sexp ctx,sexp a,sexp b)670 sexp sexp_bignum_quotient (sexp ctx, sexp a, sexp b) {
671   sexp res;
672   sexp_gc_var1(rem);
673   sexp_gc_preserve1(ctx, rem);
674   res = sexp_bignum_quot_rem(ctx, &rem, a, b);
675   sexp_gc_release1(ctx);
676   return res;
677 }
678 
sexp_bignum_remainder(sexp ctx,sexp a,sexp b)679 sexp sexp_bignum_remainder (sexp ctx, sexp a, sexp b) {
680   sexp_gc_var1(rem);
681   sexp_gc_preserve1(ctx, rem);
682   sexp_bignum_quot_rem(ctx, &rem, a, b); /* discard quotient */
683   sexp_gc_release1(ctx);
684   return rem;
685 }
686 
sexp_bignum_expt(sexp ctx,sexp a,sexp b)687 sexp sexp_bignum_expt (sexp ctx, sexp a, sexp b) {
688   sexp_sint_t e = sexp_unbox_fixnum(b);
689   sexp_sint_t abs_e;
690   if (e < 0)
691     abs_e = -e;
692   else
693     abs_e = e;
694   sexp_gc_var2(res, acc);
695   sexp_gc_preserve2(ctx, res, acc);
696   res = sexp_fixnum_to_bignum(ctx, SEXP_ONE);
697   acc = sexp_copy_bignum(ctx, NULL, a, 0);
698   for (; abs_e; abs_e>>=1, acc=sexp_bignum_mul(ctx, NULL, acc, acc))
699     if (abs_e & 1)
700       res = sexp_bignum_mul(ctx, NULL, res, acc);
701   if (e < 0)
702     res = sexp_div(ctx, sexp_fixnum_to_bignum(ctx, SEXP_ONE), res);
703   sexp_gc_release2(ctx);
704   return sexp_bignum_normalize(res);
705 }
706 
707 #if SEXP_USE_MATH
708 
709 /*
710  * a = x * 2^2n, with 0.1 <= x < 1.0 (base 2)  =>  sqrt(a) ~ 2^n
711  */
sexp_bignum_sqrt_estimate(sexp ctx,sexp a)712 sexp sexp_bignum_sqrt_estimate (sexp ctx, sexp a) {
713   sexp_uint_t alen=sexp_bignum_hi(a), adata_hi;
714   int nbits, i;
715   sexp_gc_var1(res);
716 
717   adata_hi = sexp_bignum_data(a)[alen - 1];
718   for (i = sizeof(sexp_uint_t)*8-1; i > 0; i--)
719     if (adata_hi & (1ul << i))
720       break;
721   nbits = sizeof(sexp_uint_t) * 8 * (alen - 1) + i + 1;
722 
723   sexp_gc_preserve1(ctx, res);
724   res = sexp_fixnum_to_bignum(ctx, SEXP_TWO);
725   res = sexp_bignum_expt(ctx, res, sexp_make_fixnum(nbits / 2));
726   sexp_gc_release1(ctx);
727   return res;
728 }
729 
730 #define SEXP_MAX_ACCURATE_FLONUM_SQRT 1073741824.0 /* 2^30 */
731 
732 /* Babylonian method */
sexp_bignum_sqrt(sexp ctx,sexp a,sexp * rem_out)733 sexp sexp_bignum_sqrt (sexp ctx, sexp a, sexp* rem_out) {
734   sexp_gc_var4(res, rem, tmp, tmpa);
735   if (! sexp_bignump(a)) return sexp_type_exception(ctx, NULL, SEXP_BIGNUM, a);
736   sexp_gc_preserve4(ctx, res, rem, tmp, tmpa);
737   /* initial estimate via flonum, ignoring signs */
738   if (sexp_exact_negativep(a)) {
739     tmpa = sexp_copy_bignum(ctx, NULL, a, 0);
740     a = tmpa;
741     sexp_negate(a);
742   }
743   res = sexp_make_flonum(ctx, sexp_bignum_to_double(a));
744   res = sexp_inexact_sqrt(ctx, NULL, 1, res);
745   if (sexp_flonump(res) &&
746       sexp_flonum_value(res) > SEXP_MAX_ACCURATE_FLONUM_SQRT) {
747     if (isinf(sexp_flonum_value(res)))
748       res = sexp_bignum_sqrt_estimate(ctx, a);
749     else
750       res = sexp_double_to_bignum(ctx, sexp_flonum_value(res));
751   loop:                           /* until 0 <= a - res*res < 2*res + 1 */
752     rem = sexp_mul(ctx, res, res);
753     tmp = rem = sexp_sub(ctx, a, rem);
754     if (!sexp_positivep(tmp)) goto adjust;
755     tmp = sexp_sub(ctx, tmp, SEXP_ONE);
756     tmp = sexp_quotient(ctx, tmp, SEXP_TWO);
757     tmp = sexp_compare(ctx, tmp, res);
758     if (sexp_positivep(tmp)) {
759     adjust:
760       tmp = sexp_quotient(ctx, a, res);
761       res = sexp_add(ctx, res, tmp);
762       res = sexp_quotient(ctx, res, SEXP_TWO);
763       goto loop;
764     }
765   } else {
766     if (sexp_flonump(res))
767       res = sexp_bignum_normalize(sexp_double_to_bignum(ctx, sexp_flonum_value(res)));
768     tmp = sexp_mul(ctx, res, res);
769     rem = sexp_sub(ctx, a, tmp);
770   }
771   *rem_out = sexp_bignum_normalize(rem);
772   sexp_gc_release4(ctx);
773   return sexp_bignum_normalize(res);
774 }
775 
776 #endif  /* SEXP_USE_MATH */
777 
778 /************************ ratios ******************************/
779 
780 #if SEXP_USE_RATIOS
781 
sexp_ratio_to_double(sexp ctx,sexp rat)782 double sexp_ratio_to_double (sexp ctx, sexp rat) {
783   sexp_gc_var1(quot);
784   sexp num = sexp_ratio_numerator(rat), den = sexp_ratio_denominator(rat);
785   double res = (sexp_bignump(num) ? sexp_bignum_to_double(num)
786           : sexp_fixnum_to_double(num))
787     / (sexp_bignump(den) ? sexp_bignum_to_double(den)
788        : sexp_fixnum_to_double(den));
789   if (!isfinite(res)) {
790     sexp_gc_preserve1(ctx, quot);
791     if (sexp_unbox_fixnum(sexp_compare(ctx, sexp_ratio_numerator(rat),  sexp_ratio_denominator(rat))) < 0) {
792       quot = sexp_quotient(ctx, sexp_ratio_denominator(rat),  sexp_ratio_numerator(rat));
793       res = 1 / sexp_to_double(ctx, quot);
794     } else {
795       quot = sexp_quotient(ctx, sexp_ratio_numerator(rat),  sexp_ratio_denominator(rat));
796       res = sexp_to_double(ctx, quot);
797     }
798     sexp_gc_release1(ctx);
799   }
800   return res;
801 }
802 
sexp_double_to_ratio(sexp ctx,double f)803 sexp sexp_double_to_ratio (sexp ctx, double f) {
804   int sign, i;
805   sexp_gc_var3(res, whole, scale);
806   if (f == trunc(f))
807     return sexp_bignum_normalize(sexp_double_to_bignum(ctx, f));
808   sexp_gc_preserve3(ctx, res, whole, scale);
809   whole = sexp_double_to_bignum(ctx, trunc(f));
810   res = sexp_fixnum_to_bignum(ctx, SEXP_ZERO);
811   scale = SEXP_ONE;
812   sign = (f < 0 ? -1 : 1);
813   for (i=0, f=fabs(f-trunc(f)); f != trunc(f) && i < 15; i++) {
814     res = sexp_bignum_fxmul(ctx, NULL, res, 10, 0);
815     f = f * 10;
816     res = sexp_bignum_fxadd(ctx, res, (sexp_uint_t)double_10s_digit(f));
817     f = f - trunc(f);
818     scale = sexp_mul(ctx, scale, SEXP_TEN);
819   }
820   sexp_bignum_sign(res) = sign;
821   res = sexp_bignum_normalize(res);
822   scale = sexp_bignum_normalize(scale);
823   res = sexp_make_ratio(ctx, res, scale);
824   res = sexp_ratio_normalize(ctx, res, SEXP_FALSE);
825   res = sexp_add(ctx, res, whole);
826   sexp_gc_release3(ctx);
827   return res;
828 }
829 
830 /*
831  * For conversion that does not introduce round-off error,
832  * no matter what FLT_RADIX is.
833  */
sexp_double_to_ratio_2(sexp ctx,double f)834 sexp sexp_double_to_ratio_2 (sexp ctx, double f) {
835   int sign,i;
836   sexp_gc_var3(res, whole, scale);
837   if (f == trunc(f))
838     return sexp_bignum_normalize(sexp_double_to_bignum(ctx, f));
839   sexp_gc_preserve3(ctx, res, whole, scale);
840   whole = sexp_double_to_bignum(ctx, trunc(f));
841   res = sexp_fixnum_to_bignum(ctx, SEXP_ZERO);
842   scale = SEXP_ONE;
843   sign = (f < 0 ? -1 : 1);
844   f = fabs(f-trunc(f));
845   while(f) {
846     res = sexp_bignum_fxmul(ctx, NULL, res, FLT_RADIX, 0);
847     scale = sexp_mul(ctx, scale, sexp_make_fixnum(FLT_RADIX));
848     f *= FLT_RADIX;
849     i = trunc(f);
850     if (i) {
851       f -= i;
852       res = sexp_bignum_fxadd(ctx, res, i);
853     }
854   }
855   sexp_bignum_sign(res) = sign;
856   res = sexp_bignum_normalize(res);
857   scale = sexp_bignum_normalize(scale);
858   res = sexp_make_ratio(ctx, res, scale);
859   res = sexp_ratio_normalize(ctx, res, SEXP_FALSE);
860   res = sexp_add(ctx, res, whole);
861   sexp_gc_release3(ctx);
862   return res;
863 }
864 
sexp_ratio_add(sexp ctx,sexp a,sexp b)865 sexp sexp_ratio_add (sexp ctx, sexp a, sexp b) {
866   sexp_gc_var3(res, num, den);
867   sexp_gc_preserve3(ctx, res, num, den);
868   num = sexp_mul(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(b));
869   den = sexp_mul(ctx, sexp_ratio_numerator(b), sexp_ratio_denominator(a));
870   num = sexp_add(ctx, num, den);
871   den = sexp_mul(ctx, sexp_ratio_denominator(a), sexp_ratio_denominator(b));
872   res = sexp_make_ratio(ctx, num, den);
873   res = sexp_ratio_normalize(ctx, res, SEXP_FALSE);
874   sexp_gc_release3(ctx);
875   return res;
876 }
877 
sexp_ratio_mul(sexp ctx,sexp a,sexp b)878 sexp sexp_ratio_mul (sexp ctx, sexp a, sexp b) {
879   sexp_gc_var3(res, num, den);
880   sexp_gc_preserve3(ctx, res, num, den);
881   num = sexp_mul(ctx, sexp_ratio_numerator(a), sexp_ratio_numerator(b));
882   den = sexp_mul(ctx, sexp_ratio_denominator(a), sexp_ratio_denominator(b));
883   res = sexp_make_ratio(ctx, num, den);
884   res = sexp_ratio_normalize(ctx, res, SEXP_FALSE);
885   sexp_gc_release3(ctx);
886   return res;
887 }
888 
sexp_ratio_div(sexp ctx,sexp a,sexp b)889 sexp sexp_ratio_div (sexp ctx, sexp a, sexp b) {
890   sexp_gc_var3(res, num, den);
891   sexp_gc_preserve3(ctx, res, num, den);
892   num = sexp_mul(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(b));
893   den = sexp_mul(ctx, sexp_ratio_denominator(a), sexp_ratio_numerator(b));
894   res = sexp_make_ratio(ctx, num, den);
895   res = sexp_ratio_normalize(ctx, res, SEXP_FALSE);
896   sexp_gc_release3(ctx);
897   return res;
898 }
899 
sexp_ratio_compare(sexp ctx,sexp a,sexp b)900 sexp sexp_ratio_compare (sexp ctx, sexp a, sexp b) {
901   sexp_gc_var2(a2, b2);
902   sexp_gc_preserve2(ctx, a2, b2);
903   a2 = sexp_mul(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(b));
904   b2 = sexp_mul(ctx, sexp_ratio_numerator(b), sexp_ratio_denominator(a));
905   a2 = sexp_compare(ctx, a2, b2);
906   sexp_gc_release2(ctx);
907   return a2;
908 }
909 
sexp_ratio_round(sexp ctx,sexp a)910 sexp sexp_ratio_round (sexp ctx, sexp a) {
911   sexp_gc_var2(q, r);
912   sexp_gc_preserve2(ctx, q, r);
913   q = sexp_quotient(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(a));
914   if ((sexp_ratio_denominator(a) == SEXP_TWO) && sexp_oddp(q)) {
915     q = sexp_add(ctx, q, (sexp_exact_positivep(q) ? SEXP_ONE : SEXP_NEG_ONE));
916   } else {
917     r = sexp_remainder(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(a));
918     r = sexp_mul(ctx, r, SEXP_TWO);
919     if (sexp_exact_negativep(r)) {sexp_negate(r);}
920     if (sexp_unbox_fixnum(sexp_compare(ctx, r, sexp_ratio_denominator(a))) > 0)
921       q = sexp_add(ctx, q, (sexp_exact_negativep(sexp_ratio_numerator(a)) ? SEXP_NEG_ONE : SEXP_ONE));
922   }
923   sexp_gc_release2(ctx);
924   return q;
925 }
926 
sexp_ratio_trunc(sexp ctx,sexp a)927 sexp sexp_ratio_trunc (sexp ctx, sexp a) {
928   return sexp_quotient(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(a));
929 }
930 
sexp_ratio_floor(sexp ctx,sexp a)931 sexp sexp_ratio_floor (sexp ctx, sexp a) {
932   sexp_gc_var1(q);
933   sexp_gc_preserve1(ctx, q);
934   q = sexp_quotient(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(a));
935   if (sexp_exact_negativep(sexp_ratio_numerator(a)))
936     q = sexp_add(ctx, q, SEXP_NEG_ONE);
937   sexp_gc_release1(ctx);
938   return q;
939 }
940 
sexp_ratio_ceiling(sexp ctx,sexp a)941 sexp sexp_ratio_ceiling (sexp ctx, sexp a) {
942   sexp_gc_var1(q);
943   sexp_gc_preserve1(ctx, q);
944   q = sexp_quotient(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(a));
945   if (sexp_exact_positivep(sexp_ratio_numerator(a)))
946     q = sexp_add(ctx, q, SEXP_ONE);
947   sexp_gc_release1(ctx);
948   return q;
949 }
950 
951 #endif
952 
sexp_to_double(sexp ctx,sexp x)953 double sexp_to_double (sexp ctx, sexp x) {
954   if (sexp_flonump(x))
955     return sexp_flonum_value(x);
956   else if (sexp_fixnump(x))
957     return sexp_fixnum_to_double(x);
958   else if (sexp_bignump(x))
959     return sexp_bignum_to_double(x);
960 #if SEXP_USE_RATIOS
961   else if (sexp_ratiop(x))
962     return sexp_ratio_to_double(ctx, x);
963 #endif
964   else
965     return 0.0;
966 }
967 
968 /************************ complex numbers ****************************/
969 
970 #if SEXP_USE_COMPLEX
971 
sexp_complex_copy(sexp ctx,sexp a)972 static sexp sexp_complex_copy (sexp ctx, sexp a) {
973   sexp_gc_var1(res);
974   sexp_gc_preserve1(ctx, res);
975   res = sexp_make_complex(ctx, sexp_complex_real(a), sexp_complex_imag(a));
976   if (sexp_flonump(sexp_complex_real(a)))
977     sexp_complex_real(a) = sexp_make_flonum(ctx, sexp_flonum_value(sexp_complex_real(a)));
978   else if (sexp_bignump(sexp_complex_real(a)))
979     sexp_complex_real(a) = sexp_copy_bignum(ctx, NULL, sexp_complex_real(a), 0);
980   if (sexp_flonump(sexp_complex_imag(a)))
981     sexp_complex_imag(a) = sexp_make_flonum(ctx, sexp_flonum_value(sexp_complex_imag(a)));
982   else if (sexp_bignump(sexp_complex_imag(a)))
983     sexp_complex_imag(a) = sexp_copy_bignum(ctx, NULL, sexp_complex_imag(a), 0);
984   sexp_gc_release1(ctx);
985   return res;
986 }
987 
sexp_complex_add(sexp ctx,sexp a,sexp b)988 sexp sexp_complex_add (sexp ctx, sexp a, sexp b) {
989   sexp_gc_var3(res, real, imag);
990   sexp_gc_preserve3(ctx, res, real, imag);
991   real = sexp_add(ctx, sexp_complex_real(a), sexp_complex_real(b));
992   imag = sexp_add(ctx, sexp_complex_imag(a), sexp_complex_imag(b));
993   res = sexp_make_complex(ctx, real, imag);
994   sexp_gc_release3(ctx);
995   return sexp_complex_normalize(res);
996 }
997 
sexp_complex_sub(sexp ctx,sexp a,sexp b)998 sexp sexp_complex_sub (sexp ctx, sexp a, sexp b) {
999   sexp_gc_var2(res, tmp);
1000   sexp_gc_preserve2(ctx, res, tmp);
1001   tmp = sexp_complex_copy(ctx, b);
1002   sexp_negate(sexp_complex_real(tmp));
1003   sexp_negate(sexp_complex_imag(tmp));
1004   res = sexp_complex_add(ctx, a, tmp);
1005   sexp_gc_release2(ctx);
1006   return res;
1007 }
1008 
sexp_complex_mul(sexp ctx,sexp a,sexp b)1009 sexp sexp_complex_mul (sexp ctx, sexp a, sexp b) {
1010   sexp_gc_var3(res, real, imag);
1011   sexp_gc_preserve3(ctx, res, real, imag);
1012   real = sexp_mul(ctx, sexp_complex_real(a), sexp_complex_real(b));
1013   res = sexp_mul(ctx, sexp_complex_imag(a), sexp_complex_imag(b));
1014   real = sexp_sub(ctx, real, res);
1015   imag = sexp_mul(ctx, sexp_complex_real(a), sexp_complex_imag(b));
1016   res = sexp_mul(ctx, sexp_complex_imag(a), sexp_complex_real(b));
1017   imag = sexp_add(ctx, imag, res);
1018   res = sexp_make_complex(ctx, real, imag);
1019   sexp_gc_release3(ctx);
1020   return sexp_complex_normalize(res);
1021 }
1022 
1023 /* (a + bi)    (ac + bd)      (bc - ad)    */
1024 /* -------- = -----------  + ----------- i */
1025 /* (c + di)   (c^2 + d^2)    (c^2 + d^2)   */
1026 
sexp_complex_div(sexp ctx,sexp a,sexp b)1027 sexp sexp_complex_div (sexp ctx, sexp a, sexp b) {
1028   sexp_gc_var4(res, real, imag, denom);
1029   sexp_gc_preserve4(ctx, res, real, imag, denom);
1030   /* c^2 + d^2 */
1031   denom = sexp_mul(ctx, sexp_complex_real(b), sexp_complex_real(b));
1032   res = sexp_mul(ctx, sexp_complex_imag(b), sexp_complex_imag(b));
1033   denom = sexp_add(ctx, denom, res);
1034   /* ac + bd */
1035   real = sexp_mul(ctx, sexp_complex_real(a), sexp_complex_real(b));
1036   res = sexp_mul(ctx, sexp_complex_imag(a), sexp_complex_imag(b));
1037   real = sexp_add(ctx, real, res);
1038   real = sexp_div(ctx, real, denom);
1039   /* bc - ad */
1040   imag = sexp_mul(ctx, sexp_complex_imag(a), sexp_complex_real(b));
1041   res = sexp_mul(ctx, sexp_complex_real(a), sexp_complex_imag(b));
1042   imag = sexp_sub(ctx, imag, res);
1043   imag = sexp_div(ctx, imag, denom);
1044   res = sexp_make_complex(ctx, real, imag);
1045   sexp_gc_release4(ctx);
1046   return sexp_complex_normalize(res);
1047 }
1048 
sexp_to_complex(sexp ctx,sexp x)1049 static sexp sexp_to_complex (sexp ctx, sexp x) {
1050 #if SEXP_USE_RATIOS
1051   sexp_gc_var1(tmp);
1052 #endif
1053   if (sexp_flonump(x) || sexp_fixnump(x) || sexp_bignump(x)) {
1054     return sexp_make_complex(ctx, x, SEXP_ZERO);
1055 #if SEXP_USE_RATIOS
1056   } else if (sexp_ratiop(x)) {
1057     sexp_gc_preserve1(ctx, tmp);
1058     tmp = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
1059     sexp_complex_real(tmp) = sexp_make_flonum(ctx, sexp_to_double(ctx, x));
1060     sexp_gc_release1(ctx);
1061     return tmp;
1062 #endif
1063   } else {
1064     return x;
1065   }
1066 }
1067 
sexp_complex_exp(sexp ctx,sexp z)1068 sexp sexp_complex_exp (sexp ctx, sexp z) {
1069   double e2x = exp(sexp_to_double(ctx, sexp_complex_real(z))),
1070     y = sexp_to_double(ctx, sexp_complex_imag(z));
1071   sexp_gc_var1(res);
1072   sexp_gc_preserve1(ctx, res);
1073   res = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
1074   sexp_complex_real(res) = sexp_make_flonum(ctx, e2x*cos(y));
1075   sexp_complex_imag(res) = sexp_make_flonum(ctx, e2x*sin(y));
1076   sexp_gc_release1(ctx);
1077   return res;
1078 }
1079 
sexp_complex_log(sexp ctx,sexp z)1080 sexp sexp_complex_log (sexp ctx, sexp z) {
1081   double x = sexp_to_double(ctx, sexp_complex_real(z)),
1082     y = sexp_to_double(ctx, sexp_complex_imag(z));
1083   sexp_gc_var1(res);
1084   sexp_gc_preserve1(ctx, res);
1085   res = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
1086   sexp_complex_real(res) = sexp_make_flonum(ctx, log(sqrt(x*x + y*y)));
1087   sexp_complex_imag(res) = sexp_make_flonum(ctx, atan2(y, x));
1088   sexp_gc_release1(ctx);
1089   return res;
1090 }
1091 
sexp_complex_expt(sexp ctx,sexp a,sexp b)1092 sexp sexp_complex_expt (sexp ctx, sexp a, sexp b) {
1093   sexp_gc_var1(res);
1094   sexp_gc_preserve1(ctx, res);
1095   res = sexp_to_complex(ctx, a);
1096   res = sexp_complex_log(ctx, res);
1097   res = sexp_mul(ctx, b, res);
1098   res = sexp_complex_exp(ctx, res);
1099   sexp_gc_release1(ctx);
1100   return res;
1101 }
1102 
1103 #if SEXP_USE_MATH
1104 
sexp_complex_sqrt(sexp ctx,sexp z)1105 sexp sexp_complex_sqrt (sexp ctx, sexp z) {
1106   double x = sexp_to_double(ctx, sexp_complex_real(z)),
1107     y = sexp_to_double(ctx, sexp_complex_imag(z)), r;
1108   sexp_gc_var1(res);
1109   sexp_gc_preserve1(ctx, res);
1110   r = sqrt(x*x + y*y);
1111   res = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
1112   sexp_complex_real(res) = sexp_make_flonum(ctx, sqrt((x+r)/2));
1113   sexp_complex_imag(res) = sexp_make_flonum(ctx, ((y<0||(y==0&&1/y<0))?-1:1)*sqrt((-x+r)/2));
1114   sexp_gc_release1(ctx);
1115   return res;
1116 }
1117 
sexp_complex_sin(sexp ctx,sexp z)1118 sexp sexp_complex_sin (sexp ctx, sexp z) {
1119   double x = sexp_to_double(ctx, sexp_complex_real(z)),
1120     y = sexp_to_double(ctx, sexp_complex_imag(z));
1121   sexp_gc_var1(res);
1122   sexp_gc_preserve1(ctx, res);
1123   res = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
1124   sexp_complex_real(res) = sexp_make_flonum(ctx, sin(x)*cosh(y));
1125   sexp_complex_imag(res) = sexp_make_flonum(ctx, cos(x)*sinh(y));
1126   sexp_gc_release1(ctx);
1127   return res;
1128 }
1129 
sexp_complex_cos(sexp ctx,sexp z)1130 sexp sexp_complex_cos (sexp ctx, sexp z) {
1131   double x = sexp_to_double(ctx, sexp_complex_real(z)),
1132     y = sexp_to_double(ctx, sexp_complex_imag(z));
1133   sexp_gc_var1(res);
1134   sexp_gc_preserve1(ctx, res);
1135   res = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
1136   sexp_complex_real(res) = sexp_make_flonum(ctx, cos(x)*cosh(y));
1137   sexp_complex_imag(res) = sexp_make_flonum(ctx, -sin(x)*sinh(y));
1138   sexp_gc_release1(ctx);
1139   return res;
1140 }
1141 
sexp_complex_tan(sexp ctx,sexp z)1142 sexp sexp_complex_tan (sexp ctx, sexp z) {
1143   sexp res;
1144   sexp_gc_var2(sin, cos);
1145   sexp_gc_preserve2(ctx, sin, cos);
1146   sin = sexp_complex_sin(ctx, z);
1147   cos = sexp_complex_cos(ctx, z);
1148   res = sexp_complex_div(ctx, sin, cos);
1149   sexp_gc_release2(ctx);
1150   return res;
1151 }
1152 
sexp_complex_asin(sexp ctx,sexp z)1153 sexp sexp_complex_asin (sexp ctx, sexp z) {
1154   sexp_gc_var3(res, tmp, tmp2);
1155   sexp_gc_preserve3(ctx, res, tmp, tmp2);
1156   res = sexp_complex_mul(ctx, z, z);
1157   res = sexp_sub(ctx, SEXP_ONE, res);
1158   res = sexp_sqrt(ctx, NULL, 1, res);
1159   tmp = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
1160   sexp_complex_real(tmp) = sexp_mul(ctx, SEXP_NEG_ONE, sexp_complex_imag(z));
1161   sexp_complex_imag(tmp) = sexp_complex_real(z);
1162   res = sexp_add(ctx, tmp, res);
1163   res = sexp_log(ctx, NULL, 1, res);
1164   tmp = sexp_make_complex(ctx, SEXP_ZERO, SEXP_NEG_ONE);
1165   res = sexp_mul(ctx, res, tmp);
1166   sexp_gc_release3(ctx);
1167   return res;
1168 }
1169 
sexp_complex_acos(sexp ctx,sexp z)1170 sexp sexp_complex_acos (sexp ctx, sexp z) {
1171   sexp_gc_var2(res, tmp);
1172   sexp_gc_preserve2(ctx, res, tmp);
1173   res = sexp_complex_asin(ctx, z);
1174   tmp = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
1175   sexp_complex_real(tmp) = sexp_make_flonum(ctx, acos(-1)/2);
1176   res = sexp_sub(ctx, tmp, res);
1177   sexp_gc_release2(ctx);
1178   return res;
1179 }
1180 
sexp_complex_atan(sexp ctx,sexp z)1181 sexp sexp_complex_atan (sexp ctx, sexp z) {
1182   sexp_gc_var3(res, tmp1, tmp2);
1183   sexp_gc_preserve3(ctx, res, tmp1, tmp2);
1184   tmp1 = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ONE);
1185   tmp1 = sexp_complex_mul(ctx, z, tmp1);
1186   res = sexp_make_complex(ctx, SEXP_ONE, SEXP_ZERO);
1187   res = sexp_complex_sub(ctx, res, tmp1);
1188   res = sexp_complex_log(ctx, res);
1189   tmp2 = sexp_make_complex(ctx, SEXP_ONE, SEXP_ZERO);
1190   tmp2 = sexp_complex_add(ctx, tmp2, tmp1);
1191   tmp2 = sexp_complex_log(ctx, tmp2);
1192   res = sexp_complex_sub(ctx, res, tmp2);
1193   tmp1 = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ONE);
1194   sexp_complex_imag(tmp1) = sexp_make_flonum(ctx, 0.5);
1195   res = sexp_complex_mul(ctx, res, tmp1);
1196   sexp_gc_release3(ctx);
1197   return res;
1198 }
1199 
1200 #endif
1201 #endif
1202 
1203 /****************** generic arithmetic ************************/
1204 
1205 enum sexp_number_types {
1206   SEXP_NUM_NOT = 0,
1207   SEXP_NUM_FIX,
1208   SEXP_NUM_FLO,
1209   SEXP_NUM_BIG,
1210 #if SEXP_USE_RATIOS
1211   SEXP_NUM_RAT,
1212 #endif
1213 #if SEXP_USE_COMPLEX
1214   SEXP_NUM_CPX,
1215 #endif
1216 };
1217 
1218 enum sexp_number_combs {
1219   SEXP_NUM_NOT_NOT = 0,
1220   SEXP_NUM_NOT_FIX,
1221   SEXP_NUM_NOT_FLO,
1222   SEXP_NUM_NOT_BIG,
1223 #if SEXP_USE_RATIOS
1224   SEXP_NUM_NOT_RAT,
1225 #endif
1226 #if SEXP_USE_COMPLEX
1227   SEXP_NUM_NOT_CPX,
1228 #endif
1229   SEXP_NUM_FIX_NOT,
1230   SEXP_NUM_FIX_FIX,
1231   SEXP_NUM_FIX_FLO,
1232   SEXP_NUM_FIX_BIG,
1233 #if SEXP_USE_RATIOS
1234   SEXP_NUM_FIX_RAT,
1235 #endif
1236 #if SEXP_USE_COMPLEX
1237   SEXP_NUM_FIX_CPX,
1238 #endif
1239   SEXP_NUM_FLO_NOT,
1240   SEXP_NUM_FLO_FIX,
1241   SEXP_NUM_FLO_FLO,
1242   SEXP_NUM_FLO_BIG,
1243 #if SEXP_USE_RATIOS
1244   SEXP_NUM_FLO_RAT,
1245 #endif
1246 #if SEXP_USE_COMPLEX
1247   SEXP_NUM_FLO_CPX,
1248 #endif
1249   SEXP_NUM_BIG_NOT,
1250   SEXP_NUM_BIG_FIX,
1251   SEXP_NUM_BIG_FLO,
1252   SEXP_NUM_BIG_BIG,
1253 #if SEXP_USE_RATIOS
1254   SEXP_NUM_BIG_RAT,
1255 #endif
1256 #if SEXP_USE_COMPLEX
1257   SEXP_NUM_BIG_CPX,
1258 #endif
1259 #if SEXP_USE_RATIOS
1260   SEXP_NUM_RAT_NOT,
1261   SEXP_NUM_RAT_FIX,
1262   SEXP_NUM_RAT_FLO,
1263   SEXP_NUM_RAT_BIG,
1264   SEXP_NUM_RAT_RAT,
1265 #if SEXP_USE_COMPLEX
1266   SEXP_NUM_RAT_CPX,
1267 #endif
1268 #endif
1269 #if SEXP_USE_COMPLEX
1270   SEXP_NUM_CPX_NOT,
1271   SEXP_NUM_CPX_FIX,
1272   SEXP_NUM_CPX_FLO,
1273   SEXP_NUM_CPX_BIG,
1274 #if SEXP_USE_RATIOS
1275   SEXP_NUM_CPX_RAT,
1276 #endif
1277   SEXP_NUM_CPX_CPX,
1278 #endif
1279 };
1280 
1281 static int sexp_number_types[] =
1282 #if SEXP_USE_RATIOS && SEXP_USE_COMPLEX
1283   {0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 5, 0, 0};
1284 #else
1285 #if SEXP_USE_RATIOS || SEXP_USE_COMPLEX
1286   {0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 0, 0, 0};
1287 #else
1288   {0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 0, 0, 0, 0};
1289 #endif
1290 #endif
1291 
1292 #define SEXP_NUM_NUMBER_TYPES (4 + SEXP_USE_RATIOS + SEXP_USE_COMPLEX)
1293 
sexp_number_type(sexp a)1294 static int sexp_number_type (sexp a) {
1295   return sexp_pointerp(a) ?
1296     (sexp_pointer_tag(a)<(sizeof(sexp_number_types)/sizeof(sexp_number_types[0]))
1297      ? sexp_number_types[sexp_pointer_tag(a)] : 0)
1298 #if SEXP_USE_IMMEDIATE_FLONUMS
1299     : sexp_flonump(a) ? 2
1300 #endif
1301     : sexp_fixnump(a);
1302 }
1303 
sexp_add(sexp ctx,sexp a,sexp b)1304 sexp sexp_add (sexp ctx, sexp a, sexp b) {
1305   sexp_sint_t sum;
1306   int at=sexp_number_type(a), bt=sexp_number_type(b), t;
1307   sexp r=SEXP_VOID;
1308   sexp_gc_var1(tmp);
1309   sexp_gc_preserve1(ctx, tmp);
1310   if (at > bt) {r=a; a=b; b=r; t=at; at=bt; bt=t;}
1311   switch ((at * SEXP_NUM_NUMBER_TYPES) + bt) {
1312   case SEXP_NUM_NOT_NOT: case SEXP_NUM_NOT_FIX:
1313   case SEXP_NUM_NOT_FLO: case SEXP_NUM_NOT_BIG:
1314 #if SEXP_USE_RATIOS
1315   case SEXP_NUM_NOT_RAT:
1316 #endif
1317 #if SEXP_USE_COMPLEX
1318   case SEXP_NUM_NOT_CPX:
1319 #endif
1320     r = sexp_type_exception(ctx, NULL, SEXP_NUMBER, a);
1321     break;
1322   case SEXP_NUM_FIX_FIX:
1323     sum = sexp_unbox_fixnum(a) + sexp_unbox_fixnum(b);
1324     if ((sum < SEXP_MIN_FIXNUM) || (sum > SEXP_MAX_FIXNUM))
1325       r = sexp_add(ctx, tmp=sexp_fixnum_to_bignum(ctx, a), b);
1326     else
1327       r = sexp_make_fixnum(sum);
1328     break;
1329   case SEXP_NUM_FIX_FLO:
1330     r = a == SEXP_ZERO ? b : sexp_make_flonum(ctx, sexp_fixnum_to_double(a)+sexp_flonum_value(b));
1331     break;
1332   case SEXP_NUM_FIX_BIG:
1333     r = sexp_bignum_normalize(sexp_bignum_add_fixnum(ctx, b, a));
1334     break;
1335   case SEXP_NUM_FLO_FLO:
1336     r = sexp_fp_add(ctx, a, b);
1337     break;
1338   case SEXP_NUM_FLO_BIG:
1339     r = sexp_make_flonum(ctx, sexp_flonum_value(a) + sexp_bignum_to_double(b));
1340     break;
1341   case SEXP_NUM_BIG_BIG:
1342     r = sexp_bignum_normalize(sexp_bignum_add(ctx, NULL, b, a));
1343     break;
1344 #if SEXP_USE_RATIOS
1345   case SEXP_NUM_FLO_RAT:
1346     r = sexp_make_flonum(ctx, sexp_flonum_value(a) + sexp_ratio_to_double(ctx, b));
1347     break;
1348   case SEXP_NUM_FIX_RAT:
1349   case SEXP_NUM_BIG_RAT:
1350     a = tmp = sexp_make_ratio(ctx, a, SEXP_ONE);
1351     /* ... FALLTHROUGH ... */
1352   case SEXP_NUM_RAT_RAT:
1353     r = sexp_ratio_add(ctx, a, b);
1354     break;
1355 #endif
1356 #if SEXP_USE_COMPLEX
1357 #if SEXP_USE_RATIOS
1358   case SEXP_NUM_RAT_CPX:
1359 #endif
1360   case SEXP_NUM_FLO_CPX:
1361   case SEXP_NUM_FIX_CPX:
1362   case SEXP_NUM_BIG_CPX:
1363     a = tmp = sexp_make_complex(ctx, a, SEXP_ZERO);
1364     /* ... FALLTHROUGH ... */
1365   case SEXP_NUM_CPX_CPX:
1366     r = sexp_complex_add(ctx, a, b);
1367     break;
1368 #endif
1369   }
1370   sexp_gc_release1(ctx);
1371   return r;
1372 }
1373 
sexp_sub(sexp ctx,sexp a,sexp b)1374 sexp sexp_sub (sexp ctx, sexp a, sexp b) {
1375 #if SEXP_USE_FLONUMS
1376   int negatep=0;
1377 #endif
1378   int at=sexp_number_type(a), bt=sexp_number_type(b);
1379   sexp r=SEXP_VOID;
1380   sexp_gc_var2(tmp1, tmp2);
1381   sexp_gc_preserve2(ctx, tmp1, tmp2);
1382   switch ((at * SEXP_NUM_NUMBER_TYPES) + bt) {
1383   case SEXP_NUM_NOT_NOT: case SEXP_NUM_NOT_FIX:
1384   case SEXP_NUM_NOT_FLO: case SEXP_NUM_NOT_BIG:
1385 #if SEXP_USE_RATIOS
1386   case SEXP_NUM_NOT_RAT:
1387 #endif
1388 #if SEXP_USE_COMPLEX
1389   case SEXP_NUM_NOT_CPX:
1390 #endif
1391     r = sexp_type_exception(ctx, NULL, SEXP_NUMBER, a);
1392     break;
1393   case SEXP_NUM_FIX_NOT: case SEXP_NUM_FLO_NOT: case SEXP_NUM_BIG_NOT:
1394 #if SEXP_USE_RATIOS
1395   case SEXP_NUM_RAT_NOT:
1396 #endif
1397 #if SEXP_USE_COMPLEX
1398   case SEXP_NUM_CPX_NOT:
1399 #endif
1400     r = sexp_type_exception(ctx, NULL, SEXP_NUMBER, b);
1401     break;
1402   case SEXP_NUM_FIX_FIX:
1403     r = sexp_fx_sub(a, b);      /* VM catches this case */
1404     break;
1405   case SEXP_NUM_FIX_FLO:
1406     r = sexp_make_flonum(ctx, a==SEXP_ZERO ? -sexp_flonum_value(b) : sexp_fixnum_to_double(a)-sexp_flonum_value(b));
1407     break;
1408   case SEXP_NUM_FIX_BIG:
1409     tmp1 = sexp_fixnum_to_bignum(ctx, a);
1410     r = sexp_bignum_sub(ctx, NULL, b, tmp1);
1411     sexp_negate_exact(r);
1412     r = sexp_bignum_normalize(r);
1413     break;
1414   case SEXP_NUM_FLO_FIX:
1415     r = sexp_make_flonum(ctx, sexp_flonum_value(a) - sexp_fixnum_to_double(b));
1416     break;
1417   case SEXP_NUM_FLO_FLO:
1418     r = sexp_fp_sub(ctx, a, b);
1419     break;
1420   case SEXP_NUM_FLO_BIG:
1421     r = sexp_make_flonum(ctx, sexp_flonum_value(a) - sexp_bignum_to_double(b));
1422     break;
1423   case SEXP_NUM_BIG_FIX:
1424     tmp1 = sexp_fixnum_to_bignum(ctx, b);
1425     r = sexp_bignum_normalize(sexp_bignum_sub(ctx, NULL, a, tmp1));
1426     break;
1427   case SEXP_NUM_BIG_FLO:
1428     r = sexp_make_flonum(ctx, sexp_bignum_to_double(a) - sexp_flonum_value(b));
1429     break;
1430   case SEXP_NUM_BIG_BIG:
1431     r = sexp_bignum_normalize(sexp_bignum_sub(ctx, NULL, a, b));
1432     break;
1433 #if SEXP_USE_RATIOS
1434   case SEXP_NUM_FLO_RAT:
1435     r = sexp_make_flonum(ctx, sexp_flonum_value(a) - sexp_ratio_to_double(ctx, b));
1436     break;
1437   case SEXP_NUM_RAT_FLO:
1438     r = sexp_make_flonum(ctx, sexp_ratio_to_double(ctx, a) - sexp_flonum_value(b));
1439     break;
1440   case SEXP_NUM_RAT_FIX:
1441   case SEXP_NUM_RAT_BIG:
1442     tmp1 = a; a = b; b = tmp1;
1443     negatep = 1;
1444     /* ... FALLTHROUGH ... */
1445   case SEXP_NUM_FIX_RAT:
1446   case SEXP_NUM_BIG_RAT:
1447     a = tmp1 = sexp_make_ratio(ctx, a, SEXP_ONE);
1448     /* ... FALLTHROUGH ... */
1449   case SEXP_NUM_RAT_RAT:
1450     tmp2 = sexp_make_ratio(ctx, sexp_ratio_numerator(b), sexp_ratio_denominator(b));
1451     if (sexp_bignump(sexp_ratio_numerator(tmp2)))
1452       sexp_ratio_numerator(tmp2) = sexp_copy_bignum(ctx, NULL, sexp_ratio_numerator(tmp2), 0);
1453     sexp_negate_exact(sexp_ratio_numerator(tmp2));
1454     r = sexp_ratio_add(ctx, a, tmp2);
1455     if (negatep) {
1456       if (sexp_ratiop(r)) {
1457         sexp_negate_exact(sexp_ratio_numerator(r));
1458       } else {
1459         sexp_negate_exact(r);
1460       }
1461     }
1462     break;
1463 #endif
1464 #if SEXP_USE_COMPLEX
1465 #if SEXP_USE_RATIOS
1466   case SEXP_NUM_RAT_CPX:
1467     a = tmp1 = sexp_make_flonum(ctx, sexp_ratio_to_double(ctx, a));
1468     goto complex_sub;
1469   case SEXP_NUM_CPX_RAT:
1470     b = tmp1 = sexp_make_flonum(ctx, sexp_ratio_to_double(ctx, b));
1471     /* ... FALLTHROUGH ... */
1472 #endif
1473   case SEXP_NUM_CPX_FLO:
1474   case SEXP_NUM_CPX_FIX:
1475   case SEXP_NUM_CPX_BIG:
1476     tmp1 = a; a = b; b = tmp1;
1477     negatep = 1;
1478     /* ... FALLTHROUGH ... */
1479   case SEXP_NUM_FLO_CPX:
1480   case SEXP_NUM_FIX_CPX:
1481   case SEXP_NUM_BIG_CPX:
1482     a = tmp1 = sexp_make_complex(ctx, a, SEXP_ZERO);
1483     /* ... FALLTHROUGH ... */
1484   case SEXP_NUM_CPX_CPX:
1485 #if SEXP_USE_RATIOS
1486   complex_sub:
1487 #endif
1488     r = sexp_complex_sub(ctx, a, b);
1489     if (negatep) {
1490       if (sexp_complexp(r)) {
1491         r = sexp_complex_copy(ctx, r);
1492         sexp_negate(sexp_complex_real(r));
1493         sexp_negate(sexp_complex_imag(r));
1494       } else {
1495         sexp_negate(r);
1496       }
1497     }
1498     break;
1499 #endif
1500   }
1501   sexp_gc_release2(ctx);
1502   return r;
1503 }
1504 
sexp_mul(sexp ctx,sexp a,sexp b)1505 sexp sexp_mul (sexp ctx, sexp a, sexp b) {
1506   sexp_lsint_t prod;
1507   int at=sexp_number_type(a), bt=sexp_number_type(b), t;
1508   sexp r=SEXP_VOID;
1509   sexp_gc_var1(tmp);
1510   sexp_gc_preserve1(ctx, tmp);
1511   if (at > bt) {r=a; a=b; b=r; t=at; at=bt; bt=t;}
1512   switch ((at * SEXP_NUM_NUMBER_TYPES) + bt) {
1513   case SEXP_NUM_NOT_NOT: case SEXP_NUM_NOT_FIX:
1514   case SEXP_NUM_NOT_FLO: case SEXP_NUM_NOT_BIG:
1515 #if SEXP_USE_RATIOS
1516   case SEXP_NUM_NOT_RAT:
1517 #endif
1518     r = sexp_type_exception(ctx, NULL, SEXP_NUMBER, a);
1519     break;
1520   case SEXP_NUM_FIX_FIX:
1521     prod = lsint_mul_sint(lsint_from_sint(sexp_unbox_fixnum(a)), sexp_unbox_fixnum(b));
1522     if (!lsint_is_fixnum(prod))
1523       r = sexp_mul(ctx, tmp=sexp_fixnum_to_bignum(ctx, a), b);
1524     else
1525       r = sexp_make_fixnum(lsint_to_sint(prod));
1526     break;
1527   case SEXP_NUM_FIX_FLO:
1528     r = (a==SEXP_ZERO ? a : sexp_make_flonum(ctx, sexp_fixnum_to_double(a)*sexp_flonum_value(b)));
1529     break;
1530   case SEXP_NUM_FIX_BIG:
1531     r = sexp_bignum_fxmul(ctx, NULL, b, sexp_unbox_fx_abs(a), 0);
1532     sexp_bignum_sign(r) = sexp_fx_sign(a) * sexp_bignum_sign(b);
1533     r = sexp_bignum_normalize(r);
1534     break;
1535   case SEXP_NUM_FLO_FLO:
1536     r = sexp_fp_mul(ctx, a, b);
1537     break;
1538   case SEXP_NUM_FLO_BIG:
1539     r = sexp_make_flonum(ctx, sexp_flonum_value(a) * sexp_bignum_to_double(b));
1540     break;
1541   case SEXP_NUM_BIG_BIG:
1542     r = sexp_bignum_normalize(sexp_bignum_mul(ctx, NULL, a, b));
1543     break;
1544 #if SEXP_USE_RATIOS
1545   case SEXP_NUM_FLO_RAT:
1546     r = sexp_make_flonum(ctx, sexp_flonum_value(a) * sexp_ratio_to_double(ctx, b));
1547     break;
1548   case SEXP_NUM_FIX_RAT:
1549   case SEXP_NUM_BIG_RAT:
1550     a = tmp = sexp_make_ratio(ctx, a, SEXP_ONE);
1551     /* ... FALLTHROUGH ... */
1552   case SEXP_NUM_RAT_RAT:
1553     r = sexp_ratio_mul(ctx, a, b);
1554     break;
1555 #endif
1556 #if SEXP_USE_COMPLEX
1557 #if SEXP_USE_RATIOS
1558   case SEXP_NUM_RAT_CPX:
1559 #endif
1560   case SEXP_NUM_FLO_CPX:
1561   case SEXP_NUM_FIX_CPX:
1562   case SEXP_NUM_BIG_CPX:
1563     a = tmp = sexp_make_complex(ctx, a, SEXP_ZERO);
1564     /* ... FALLTHROUGH ... */
1565   case SEXP_NUM_CPX_CPX:
1566     r = sexp_complex_mul(ctx, a, b);
1567     break;
1568 #endif
1569   }
1570   sexp_gc_release1(ctx);
1571   return r;
1572 }
1573 
sexp_div(sexp ctx,sexp a,sexp b)1574 sexp sexp_div (sexp ctx, sexp a, sexp b) {
1575   int at=sexp_number_type(a), bt=sexp_number_type(b);
1576 #if ! SEXP_USE_RATIOS
1577   double f;
1578 #endif
1579   sexp r=SEXP_VOID;
1580   sexp_gc_var1(tmp);
1581   sexp_gc_preserve1(ctx, tmp);
1582   switch ((at * SEXP_NUM_NUMBER_TYPES) + bt) {
1583   case SEXP_NUM_NOT_NOT: case SEXP_NUM_NOT_FIX:
1584   case SEXP_NUM_NOT_FLO: case SEXP_NUM_NOT_BIG:
1585 #if SEXP_USE_RATIOS
1586   case SEXP_NUM_NOT_RAT:
1587 #endif
1588     r = sexp_type_exception(ctx, NULL, SEXP_NUMBER, a);
1589     break;
1590   case SEXP_NUM_FIX_NOT: case SEXP_NUM_FLO_NOT: case SEXP_NUM_BIG_NOT:
1591 #if SEXP_USE_RATIOS
1592   case SEXP_NUM_RAT_NOT:
1593 #endif
1594     r = sexp_type_exception(ctx, NULL, SEXP_NUMBER, b);
1595     break;
1596   case SEXP_NUM_FIX_FIX:
1597 #if SEXP_USE_RATIOS
1598     tmp = sexp_make_ratio(ctx, a, b);
1599     r = sexp_ratio_normalize(ctx, tmp, SEXP_FALSE);
1600 #else
1601     f = sexp_fixnum_to_double(a) / sexp_fixnum_to_double(b);
1602     r = ((f == trunc(f)) ? sexp_make_fixnum((sexp_sint_t)f)
1603          : sexp_make_flonum(ctx, f));
1604 #endif
1605     break;
1606   case SEXP_NUM_FIX_FLO:
1607     r = sexp_make_flonum(ctx, sexp_fixnum_to_double(a)/sexp_flonum_value(b));
1608     break;
1609   case SEXP_NUM_FIX_BIG:
1610 #if SEXP_USE_RATIOS
1611     tmp = sexp_make_ratio(ctx, a, b);
1612     r = sexp_ratio_normalize(ctx, tmp, SEXP_FALSE);
1613 #else
1614     r = sexp_make_flonum(ctx, sexp_fixnum_to_double(a)/sexp_bignum_to_double(b));
1615 #endif
1616     break;
1617   case SEXP_NUM_FLO_FIX:
1618     r = sexp_make_flonum(ctx, sexp_flonum_value(a)/sexp_fixnum_to_double(b));
1619     break;
1620   case SEXP_NUM_FLO_FLO:
1621     r = sexp_fp_div(ctx, a, b);
1622     break;
1623   case SEXP_NUM_FLO_BIG:
1624     r = sexp_make_flonum(ctx, sexp_flonum_value(a) / sexp_bignum_to_double(b));
1625     break;
1626   case SEXP_NUM_BIG_FIX:
1627 #if SEXP_USE_RATIOS
1628     tmp = sexp_make_ratio(ctx, a, b);
1629     r = sexp_ratio_normalize(ctx, tmp, SEXP_FALSE);
1630     break;
1631 #else
1632     b = tmp = sexp_fixnum_to_bignum(ctx, b);
1633 #endif
1634     /* ... FALLTHROUGH if ! SEXP_USE_RATIOS ... */
1635   case SEXP_NUM_BIG_BIG:
1636 #if SEXP_USE_RATIOS
1637     tmp = sexp_make_ratio(ctx, a, b);
1638     r = sexp_ratio_normalize(ctx, tmp, SEXP_FALSE);
1639 #else
1640     r = sexp_bignum_quot_rem(ctx, &tmp, a, b);
1641     if (sexp_bignum_normalize(tmp) != SEXP_ZERO)
1642       r = sexp_make_flonum(ctx, sexp_bignum_to_double(a)
1643                            / sexp_bignum_to_double(b));
1644     else
1645       r = sexp_bignum_normalize(r);
1646 #endif
1647     break;
1648   case SEXP_NUM_BIG_FLO:
1649     r = sexp_make_flonum(ctx, sexp_bignum_to_double(a) / sexp_flonum_value(b));
1650     break;
1651 #if SEXP_USE_RATIOS
1652   case SEXP_NUM_FLO_RAT:
1653     r = sexp_make_flonum(ctx, sexp_flonum_value(a) / sexp_ratio_to_double(ctx, b));
1654     break;
1655   case SEXP_NUM_RAT_FLO:
1656     r = sexp_make_flonum(ctx, sexp_ratio_to_double(ctx, a) / sexp_flonum_value(b));
1657     break;
1658   case SEXP_NUM_RAT_FIX:
1659   case SEXP_NUM_RAT_BIG:
1660     b = tmp = sexp_make_ratio(ctx, b, SEXP_ONE);
1661     /* ... FALLTHROUGH ... */
1662   case SEXP_NUM_FIX_RAT:
1663   case SEXP_NUM_BIG_RAT:
1664     if (!sexp_ratiop(a))
1665       a = tmp = sexp_make_ratio(ctx, a, SEXP_ONE);
1666     /* ... FALLTHROUGH ... */
1667   case SEXP_NUM_RAT_RAT:
1668     r = sexp_ratio_div(ctx, a, b);
1669     break;
1670 #endif
1671 #if SEXP_USE_COMPLEX
1672 #if SEXP_USE_RATIOS
1673   case SEXP_NUM_CPX_RAT:
1674     b = tmp = sexp_make_flonum(ctx, sexp_ratio_to_double(ctx, b));
1675     /* ... FALLTHROUGH ... */
1676 #endif
1677   case SEXP_NUM_CPX_FLO:
1678   case SEXP_NUM_CPX_FIX:
1679   case SEXP_NUM_CPX_BIG:
1680     b = tmp = sexp_make_complex(ctx, b, SEXP_ZERO);
1681     /* ... FALLTHROUGH ... */
1682 #if SEXP_USE_RATIOS
1683   case SEXP_NUM_RAT_CPX:
1684     if (sexp_ratiop(a))
1685       a = tmp = sexp_make_flonum(ctx, sexp_ratio_to_double(ctx, a));
1686     /* ... FALLTHROUGH ... */
1687 #endif
1688   case SEXP_NUM_FLO_CPX:
1689   case SEXP_NUM_FIX_CPX:
1690   case SEXP_NUM_BIG_CPX:
1691     if (!sexp_complexp(a))
1692       a = tmp = sexp_make_complex(ctx, a, SEXP_ZERO);
1693     /* ... FALLTHROUGH ... */
1694   case SEXP_NUM_CPX_CPX:
1695     r = sexp_complex_div(ctx, a, b);
1696     break;
1697 #endif
1698   }
1699   sexp_gc_release1(ctx);
1700   return r;
1701 }
1702 
sexp_to_inexact(sexp ctx,sexp a)1703 sexp sexp_to_inexact (sexp ctx, sexp a) {
1704   if (sexp_fixnump(a)) return sexp_fixnum_to_flonum(ctx, a);
1705   if (sexp_bignump(a)) return sexp_make_flonum(ctx, sexp_bignum_to_double(a));
1706   return a;
1707 }
1708 
sexp_quotient(sexp ctx,sexp a,sexp b)1709 sexp sexp_quotient (sexp ctx, sexp a, sexp b) {
1710   int at=sexp_number_type(a), bt=sexp_number_type(b);
1711   sexp r=SEXP_VOID;
1712   sexp_gc_var1(tmp);
1713   if (b == SEXP_ONE) return a;
1714   sexp_gc_preserve1(ctx, tmp);
1715   switch ((at * SEXP_NUM_NUMBER_TYPES) + bt) {
1716   case SEXP_NUM_NOT_NOT: case SEXP_NUM_NOT_FIX:
1717   case SEXP_NUM_NOT_FLO: case SEXP_NUM_NOT_BIG:
1718     r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, a);
1719     break;
1720   case SEXP_NUM_FIX_NOT: case SEXP_NUM_FLO_NOT: case SEXP_NUM_BIG_NOT:
1721     r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, b);
1722     break;
1723   case SEXP_NUM_FLO_FIX: case SEXP_NUM_FLO_FLO: case SEXP_NUM_FLO_BIG:
1724 #if SEXP_USE_RATIOS
1725   case SEXP_NUM_FLO_RAT:
1726 #endif
1727     if (sexp_flonum_value(a) != trunc(sexp_flonum_value(a))) {
1728       r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, a);
1729     } else {
1730       tmp = sexp_bignum_normalize(sexp_double_to_bignum(ctx, sexp_flonum_value(a)));
1731       tmp = sexp_quotient(ctx, tmp, b);
1732       r = sexp_to_inexact(ctx, tmp);
1733     }
1734     break;
1735 #if SEXP_USE_RATIOS
1736   case SEXP_NUM_RAT_FIX: case SEXP_NUM_RAT_BIG: case SEXP_NUM_RAT_RAT:
1737 #endif
1738 #if SEXP_USE_COMPLEX
1739   case SEXP_NUM_FLO_CPX: case SEXP_NUM_CPX_FIX: case SEXP_NUM_CPX_FLO:
1740   case SEXP_NUM_CPX_BIG: case SEXP_NUM_CPX_CPX:
1741 #if SEXP_USE_RATIOS
1742   case SEXP_NUM_CPX_RAT:
1743 #endif
1744 #endif
1745     r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, a);
1746     break;
1747   case SEXP_NUM_FIX_FLO: case SEXP_NUM_BIG_FLO:
1748 #if SEXP_USE_RATIOS
1749   case SEXP_NUM_RAT_FLO:
1750 #endif
1751     if (sexp_flonum_value(b) != trunc(sexp_flonum_value(b))) {
1752       r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, b);
1753     } else {
1754       tmp = sexp_bignum_normalize(sexp_double_to_bignum(ctx, sexp_flonum_value(b)));
1755       tmp = sexp_quotient(ctx, a, tmp);
1756       r = sexp_to_inexact(ctx, tmp);
1757     }
1758     break;
1759 #if SEXP_USE_RATIOS
1760   case SEXP_NUM_FIX_RAT: case SEXP_NUM_BIG_RAT:
1761 #endif
1762 #if SEXP_USE_COMPLEX
1763   case SEXP_NUM_FIX_CPX: case SEXP_NUM_BIG_CPX:
1764 #endif
1765     r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, b);
1766     break;
1767   case SEXP_NUM_FIX_FIX:
1768     r = sexp_fx_div(a, b);
1769     break;
1770   case SEXP_NUM_FIX_BIG:
1771     r = SEXP_ZERO;
1772     break;
1773   case SEXP_NUM_BIG_FIX:
1774     b = tmp = sexp_fixnum_to_bignum(ctx, b);
1775     /* ... FALLTHROUGH ... */
1776   case SEXP_NUM_BIG_BIG:
1777     r = sexp_bignum_normalize(sexp_bignum_quotient(ctx, a, b));
1778     break;
1779   }
1780   sexp_gc_release1(ctx);
1781   return r;
1782 }
1783 
sexp_remainder(sexp ctx,sexp a,sexp b)1784 sexp sexp_remainder (sexp ctx, sexp a, sexp b) {
1785   int at=sexp_number_type(a), bt=sexp_number_type(b);
1786   sexp r=SEXP_VOID;
1787   sexp_gc_var1(tmp);
1788   if (b == SEXP_ONE) return SEXP_ZERO;
1789   sexp_gc_preserve1(ctx, tmp);
1790   switch ((at * SEXP_NUM_NUMBER_TYPES) + bt) {
1791   case SEXP_NUM_NOT_NOT: case SEXP_NUM_NOT_FIX:
1792   case SEXP_NUM_NOT_FLO: case SEXP_NUM_NOT_BIG:
1793     r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, a);
1794     break;
1795   case SEXP_NUM_FIX_NOT: case SEXP_NUM_FLO_NOT: case SEXP_NUM_BIG_NOT:
1796     r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, b);
1797     break;
1798   case SEXP_NUM_FLO_FIX: case SEXP_NUM_FLO_FLO: case SEXP_NUM_FLO_BIG:
1799 #if SEXP_USE_RATIOS
1800   case SEXP_NUM_FLO_RAT:
1801 #endif
1802     if (isinf(sexp_flonum_value(a)) ||
1803         sexp_flonum_value(a) != trunc(sexp_flonum_value(a))) {
1804       r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, a);
1805     } else if (bt == SEXP_NUM_FLO && isinf(sexp_flonum_value(b))) {
1806       r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, b);
1807     } else {
1808       tmp = sexp_bignum_normalize(sexp_double_to_bignum(ctx, sexp_flonum_value(a)));
1809       tmp = sexp_remainder(ctx, tmp, b);
1810       r = sexp_to_inexact(ctx, tmp);
1811     }
1812     break;
1813 #if SEXP_USE_RATIOS
1814   case SEXP_NUM_RAT_FIX: case SEXP_NUM_RAT_BIG: case SEXP_NUM_RAT_RAT:
1815 #endif
1816 #if SEXP_USE_COMPLEX
1817   case SEXP_NUM_FLO_CPX: case SEXP_NUM_CPX_FIX: case SEXP_NUM_CPX_FLO:
1818   case SEXP_NUM_CPX_BIG: case SEXP_NUM_CPX_CPX:
1819 #if SEXP_USE_RATIOS
1820   case SEXP_NUM_CPX_RAT:
1821 #endif
1822 #endif
1823     r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, a);
1824     break;
1825   case SEXP_NUM_FIX_FLO: case SEXP_NUM_BIG_FLO:
1826 #if SEXP_USE_RATIOS
1827   case SEXP_NUM_RAT_FLO:
1828 #endif
1829     if (isinf(sexp_flonum_value(b)) ||
1830         sexp_flonum_value(b) != trunc(sexp_flonum_value(b))) {
1831       r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, b);
1832     } else {
1833       tmp = sexp_bignum_normalize(sexp_double_to_bignum(ctx, sexp_flonum_value(b)));
1834       tmp = sexp_remainder(ctx, a, tmp);
1835       r = sexp_to_inexact(ctx, tmp);
1836     }
1837     break;
1838 #if SEXP_USE_RATIOS
1839   case SEXP_NUM_FIX_RAT: case SEXP_NUM_BIG_RAT:
1840 #endif
1841 #if SEXP_USE_COMPLEX
1842   case SEXP_NUM_FIX_CPX: case SEXP_NUM_BIG_CPX:
1843 #endif
1844     r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, b);
1845     break;
1846   case SEXP_NUM_FIX_FIX:
1847     r = sexp_fx_rem(a, b);
1848     break;
1849   case SEXP_NUM_FIX_BIG:
1850     r = a;
1851     break;
1852   case SEXP_NUM_BIG_FIX:
1853     r = sexp_bignum_fxrem(ctx, a, sexp_unbox_fixnum(b));
1854     break;
1855   case SEXP_NUM_BIG_BIG:
1856     r = sexp_bignum_normalize(sexp_bignum_remainder(ctx, a, b));
1857     break;
1858   }
1859   sexp_gc_release1(ctx);
1860   return r;
1861 }
1862 
sexp_compare(sexp ctx,sexp a,sexp b)1863 sexp sexp_compare (sexp ctx, sexp a, sexp b) {
1864   int at=sexp_number_type(a), bt=sexp_number_type(b);
1865   sexp r=SEXP_VOID;
1866   double f, g;
1867   sexp_gc_var1(tmp);
1868   sexp_gc_preserve1(ctx, tmp);
1869   if (at > bt) {
1870     r = sexp_compare(ctx, b, a);
1871     sexp_negate(r);
1872   } else {
1873     switch ((at * SEXP_NUM_NUMBER_TYPES) + bt) {
1874     case SEXP_NUM_NOT_NOT: case SEXP_NUM_NOT_FIX:
1875     case SEXP_NUM_NOT_FLO: case SEXP_NUM_NOT_BIG:
1876 #if SEXP_USE_COMPLEX
1877     case SEXP_NUM_CPX_CPX: case SEXP_NUM_CPX_FIX:
1878     case SEXP_NUM_CPX_FLO: case SEXP_NUM_CPX_BIG:
1879 #if SEXP_USE_RATIOS
1880     case SEXP_NUM_CPX_RAT:
1881 #endif
1882 #endif
1883       r = sexp_type_exception(ctx, NULL, SEXP_NUMBER, a);
1884       break;
1885     case SEXP_NUM_FIX_FIX:
1886       r = sexp_make_fixnum(sexp_unbox_fixnum(a) - sexp_unbox_fixnum(b));
1887       break;
1888     case SEXP_NUM_FIX_FLO:
1889       f = sexp_fixnum_to_double(a);
1890       g = sexp_flonum_value(b);
1891       if (isnan(g))
1892         r = sexp_xtype_exception(ctx, NULL, "can't compare NaN", b);
1893       else
1894         r = sexp_make_fixnum(f < g ? -1 : f == g ? 0 : 1);
1895       break;
1896     case SEXP_NUM_FIX_BIG:
1897       if ((sexp_bignum_hi(b) > 1) ||
1898           (sexp_bignum_data(b)[0] > SEXP_MAX_FIXNUM))
1899         r = sexp_make_fixnum(sexp_bignum_sign(b) < 0 ? 1 : -1);
1900       else
1901         r = sexp_make_fixnum(sexp_unbox_fixnum(a) - (sexp_sint_t)sexp_bignum_data(b)[0]);
1902       break;
1903     case SEXP_NUM_FLO_FLO:
1904       f = sexp_flonum_value(a);
1905       g = sexp_flonum_value(b);
1906       if (isnan(f))
1907         r = sexp_xtype_exception(ctx, NULL, "can't compare NaN", a);
1908       else if (isnan(g))
1909         r = sexp_xtype_exception(ctx, NULL, "can't compare NaN", b);
1910       else
1911         r = sexp_make_fixnum(f < g ? -1 : f == g ? 0 : 1);
1912       break;
1913     case SEXP_NUM_FLO_BIG:
1914       f = sexp_flonum_value(a);
1915       if (isinf(f)) {
1916         r = f > 0 ? SEXP_ONE : SEXP_NEG_ONE;
1917         break;
1918       } else if (isnan(f)) {
1919         r = sexp_xtype_exception(ctx, NULL, "can't compare NaN", a);
1920         break;
1921       } else {
1922         a = tmp = sexp_double_to_bignum(ctx, f);
1923       }
1924       /* ... FALLTHROUGH ... */
1925     case SEXP_NUM_BIG_BIG:
1926       r = sexp_make_fixnum(sexp_bignum_compare(a, b));
1927       break;
1928 #if SEXP_USE_RATIOS
1929     case SEXP_NUM_FLO_RAT:
1930       f = sexp_flonum_value(a);
1931       if (isinf(f)) {
1932         r = f > 0 ? SEXP_ONE : SEXP_NEG_ONE;
1933       } else if (isnan(f)) {
1934         r = sexp_xtype_exception(ctx, NULL, "can't compare NaN", a);
1935       } else {
1936         g = sexp_ratio_to_double(ctx, b);
1937         r = sexp_make_fixnum(f < g ? -1 : f == g ? 0 : 1);
1938       }
1939       break;
1940     case SEXP_NUM_FIX_RAT:
1941     case SEXP_NUM_BIG_RAT:
1942       a = tmp = sexp_make_ratio(ctx, a, SEXP_ONE);
1943       /* ... FALLTHROUGH ... */
1944     case SEXP_NUM_RAT_RAT:
1945       r = sexp_ratio_compare(ctx, a, b);
1946       break;
1947 #endif
1948     }
1949   }
1950   sexp_gc_release1(ctx);
1951   return r;
1952 }
1953 
1954 #endif
1955