source: project/release/4/numbers/trunk/numbers-c.c @ 32555

Last change on this file since 32555 was 32555, checked in by sjamaan, 5 years ago

numbers: Improve stability of Burnikel-Ziegler performance.

This changes Burnikel-Ziegler division's threshold check to look at
the *difference* between the two arguments instead of the size of the
smallest one. This idea is taken from MCA, algorithm 1.9.

This seems to make a big difference on the "pidigits" benchmark from
the Great Language Shootout game, which somehow triggers such
pathological behaviour of the B/Z division routine that using plain
gradebook division would be up to three times faster on that
benchmark. Other benchmarks are mostly unaffected by this issue.

Many thanks to "balkenbrij" on Reddit for finding this issue!

File size: 112.5 KB
Line 
1/* numbers-c.c
2 *
3 * Copyright 2010-2015 The CHICKEN Team
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are
7 * met:
8 *
9 *    1. Redistributions of source code must retain the above copyright
10 *       notice, this list of conditions and the following disclaimer.
11 *
12 *    2. Redistributions in binary form must reproduce the above
13 *       copyright notice, this list of conditions and the following
14 *       disclaimer in the documentation and/or other materials provided
15 *       with the distribution.
16 *
17 *    3. The name of the author may not be used to endorse or promote
18 *       products derived from this software without specific prior
19 *       written permission.
20 *
21 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
22 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
23 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
24 * DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT,
25 * INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
26 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
27 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
28 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
29 * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
30 * IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31 * POSSIBILITY OF SUCH DAMAGE.
32 */
33
34#include <assert.h>
35#include <errno.h>
36#include <math.h> /* frexp() */
37
38#define nmax(x, y)     ((x) > (y) ? (x) : (y)) /* From runtime.c */
39#define nmin(x, y)     ((x) < (y) ? (x) : (y)) /* From runtime.c */
40#define free_tmp_bignum(b)    C_free((void *)(b))
41
42static void *tags;
43
44#include "numbers-c.h"
45
46/* The bignum digit representation is fullword- little endian, so on
47 * LE machines the halfdigits are numbered in the same order.  On BE
48 * machines, we must swap the odd and even positions.
49 */
50#ifdef C_BIG_ENDIAN
51#define C_uhword_ref(x, p)           ((C_uhword *)(x))[(p)^1]
52#else
53#define C_uhword_ref(x, p)           ((C_uhword *)(x))[(p)]
54#endif
55#define C_uhword_set(x, p, d)        (C_uhword_ref(x,p) = (d))
56
57static C_word init_tags(___scheme_value tagvec);
58static void bignum_negate_2(C_word c, C_word self, C_word new_big) C_noret;
59static C_word rat_cmp(C_word x, C_word y);
60static C_word basic_cmp(C_word x, C_word y, char *loc, int eqp);
61static void allocate_bignum_2(C_word c, C_word self, C_word bigvec) C_noret;
62static C_word allocate_tmp_bignum(C_word size, C_word negp, C_word initp);
63static C_uword bignum_digits_destructive_scale_up_with_carry(C_uword *start, C_uword *end, C_uword factor, C_uword carry);
64static C_uword bignum_digits_destructive_scale_down(C_uword *start, C_uword *end, C_uword denominator);
65static C_uword bignum_digits_destructive_shift_right(C_uword *start, C_uword *end, int shift_right, int negp);
66static C_uword bignum_digits_destructive_shift_left(C_uword *start, C_uword *end, int shift_left);
67static void bignum_plus_unsigned(C_word k, C_word x, C_word y, C_word negp) C_noret;
68static void bignum_plus_unsigned_2(C_word c, C_word self, C_word result) C_noret;
69static C_word int_flo_cmp(C_word intnum, C_word flonum);
70static C_word flo_int_cmp(C_word flonum, C_word intnum);
71static C_word rat_flo_cmp(C_word ratnum, C_word flonum);
72static C_word flo_rat_cmp(C_word flonum, C_word ratnum);
73static int bignum_cmp_unsigned(C_word x, C_word y);
74static void bignum_minus_unsigned(C_word k, C_word x, C_word y) C_noret;
75static void bignum_minus_unsigned_2(C_word c, C_word self, C_word result) C_noret;
76static void integer_times_2(C_word c, C_word self, C_word new_big) C_noret;
77static C_regparm void bignum_digits_multiply(C_word x, C_word y, C_word result);
78static void bignum_times_bignum_unsigned(C_word k, C_word x, C_word y, C_word negp) C_noret;
79static void bignum_times_bignum_unsigned_2(C_word c, C_word self, C_word result) C_noret;
80static void digits_to_integer_2(C_word c, C_word self, C_word result) C_noret;
81static C_regparm C_word str_to_bignum(C_word bignum, char *str, char *str_end, int radix);
82static void bignum_to_str_2(C_word c, C_word self, C_word string) C_noret;
83static void fabs_frexp_to_digits(C_uword exp, double sign, C_uword *start, C_uword *scan);
84static C_word flo_to_tmp_bignum(C_word x);
85static void flo_to_int_2(C_word c, C_word self, C_word result) C_noret;
86static void bignum_actual_shift(C_word c, C_word self, C_word result) C_noret;
87static void bignum_actual_extraction(C_word c, C_word self, C_word result) C_noret;
88static void bignum_random_2(C_word c, C_word self, C_word result) C_noret;
89static void bignum_bitwise_and_2(C_word c, C_word self, C_word result) C_noret;
90static void bignum_bitwise_ior_2(C_word c, C_word self, C_word result) C_noret;
91static void bignum_bitwise_xor_2(C_word c, C_word self, C_word result) C_noret;
92static void bignum_digits_destructive_negate(C_word result);
93static C_regparm void basic_divrem(C_word c, C_word self, C_word k, C_word x, C_word y, C_word return_r, C_word return_q) C_noret;
94static C_regparm void integer_divrem(C_word c, C_word self, C_word k, C_word x, C_word y, C_word return_q, C_word return_r) C_noret;
95static C_regparm void bignum_divrem(C_word c, C_word self, C_word k, C_word x, C_word y, C_word return_q, C_word return_r) C_noret;
96static void divrem_intflo_2(C_word c, C_word self, ...) C_noret;
97static void bignum_divrem_fixnum_2(C_word c, C_word self, C_word negated_big) C_noret;
98static void bignum_negneg_bitwise_op(C_word c, C_word self, C_word result) C_noret;
99static void bignum_posneg_bitwise_op(C_word c, C_word self, C_word result) C_noret;
100static void bignum_pospos_bitwise_op(C_word c, C_word self, C_word result) C_noret;
101static C_word bignum_remainder_unsigned_halfdigit(C_word num, C_word den);
102static void bignum_destructive_divide_unsigned_small(C_word c, C_word self, C_word quotient);
103static void bignum_divide_2_unsigned(C_word c, C_word self, C_word quotient) C_noret;
104static void bignum_divide_2_unsigned_2(C_word c, C_word self, C_word remainder) C_noret;
105static void bignum_destructive_divide_normalized(C_word big_u, C_word big_v, C_word big_q);
106
107static void barf(int code, char *loc, ...) C_noret;
108static void try_extended_number(char *ext_proc_name, C_word c, C_word k, ...) C_noret;
109
110/* Use high numbers for when CHICKEN 4 grows more error codes! */
111#define C_BAD_ARGUMENT_TYPE_NO_REAL_ERROR                99948
112#define C_BAD_ARGUMENT_TYPE_COMPLEX_NO_ORDERING_ERROR    99949
113
114/* XXX THIS IS DUPLICATED HERE FROM runtime.c, but should be ripped out */
115static void barf(int code, char *loc, ...)
116{
117  char *msg;
118  int c, i;
119  va_list v;
120  /* Just take a random size that will "always" fit... */
121  C_word err, ab[C_SIZEOF_STRING(64)], *a = ab;
122
123  err = C_lookup_symbol(C_intern2(&a, C_text("numbers#@error-hook")));
124
125  switch(code) {
126  case C_BAD_ARGUMENT_COUNT_ERROR:
127    msg = C_text("bad argument count");
128    c = 3;
129    break;
130
131  case C_BAD_MINIMUM_ARGUMENT_COUNT_ERROR:
132    msg = C_text("too few arguments");
133    c = 3;
134    break;
135 
136  case C_BAD_ARGUMENT_TYPE_ERROR:
137    msg = C_text("bad argument type");
138    c = 1;
139    break;
140
141  case C_DIVISION_BY_ZERO_ERROR:
142    msg = C_text("division by zero");
143    c = 0;
144    break;
145
146  case C_OUT_OF_RANGE_ERROR:
147    msg = C_text("out of range");
148    c = 2;
149    break;
150
151  case C_CANT_REPRESENT_INEXACT_ERROR:
152    msg = C_text("inexact number cannot be represented as an exact number");
153    c = 1;
154    break;
155
156  case C_BAD_ARGUMENT_TYPE_NO_FIXNUM_ERROR:
157    msg = C_text("bad argument type - not a fixnum");
158    c = 1;
159    break;
160
161  case C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR:
162    msg = C_text("bad argument type - not a number");
163    c = 1;
164    break;
165
166  case C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR:
167    msg = C_text("bad argument type - not an integer");
168    c = 1;
169    break;
170
171  case C_BAD_ARGUMENT_TYPE_NO_UINTEGER_ERROR:
172    msg = C_text("bad argument type - not an unsigned integer");
173    c = 1;
174    break;
175
176  case C_BAD_ARGUMENT_TYPE_NO_FLONUM_ERROR:
177    msg = C_text("bad argument type - not a flonum");
178    c = 1;
179    break;
180
181  case C_BAD_ARGUMENT_TYPE_BAD_BASE_ERROR:
182    msg = C_text("bad argument type - invalid base");
183    c = 1;
184    break;
185
186  case C_BAD_ARGUMENT_TYPE_NO_REAL_ERROR:
187    msg = C_text("bad argument type - not a real");
188    c = 1;
189    break;
190
191  case C_BAD_ARGUMENT_TYPE_COMPLEX_NO_ORDERING_ERROR:
192    msg = C_text("bad argument type - complex number has no ordering");
193    c = 1;
194    break;
195   
196  default:
197    fprintf(stderr, "Unknown error");
198    abort();
199  }
200 
201  if(!C_immediatep(err)) {
202    C_save(C_fix(code));
203   
204    C_save(C_intern2(&a, loc));
205   
206    va_start(v, loc);
207    i = c;
208
209    while(i--)
210      C_save(va_arg(v, C_word));
211
212    va_end(v);
213    /* No continuation is passed: '##sys#error-hook' may not return: */
214    C_do_apply(c + 2, C_block_item(err, 0), C_SCHEME_UNDEFINED); 
215  } else {
216    fprintf(stderr, "No error hook!");
217    abort();
218  }
219}
220
221void C_not_an_integer_error(char *loc, C_word x)
222{
223  barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, loc, x);
224}
225
226/* This exists because numbers_div_by_zero_error doesn't exist in CHICKENs < 4.9.0 */
227void numbers_div_by_zero_error(char *loc)
228{
229  barf(C_DIVISION_BY_ZERO_ERROR, loc);
230}
231
232
233/* Never use extended number hook procedure names longer than this! */
234/* Current longest name: numbers#@bignum-2-divrem-burnikel-ziegler */
235#define MAX_EXTNUM_HOOK_NAME 64
236
237/* This exists so that we don't have to create any extra closures */
238static void try_extended_number(char *ext_proc_name, C_word c, C_word k, ...)
239{
240  static C_word ab[C_SIZEOF_STRING(MAX_EXTNUM_HOOK_NAME)];
241  int i;
242  va_list v;
243  C_word ext_proc_sym, ext_proc = C_SCHEME_FALSE, *a = ab;
244
245  ext_proc_sym = C_lookup_symbol(C_intern2(&a, ext_proc_name));
246
247  if(!C_immediatep(ext_proc_sym))
248    ext_proc = C_block_item(ext_proc_sym, 0);
249
250  if (!C_immediatep(ext_proc) && C_closurep(ext_proc)) {
251    va_start(v, k);
252    i = c - 1;
253
254    while(i--)
255      C_save(va_arg(v, C_word));
256
257    va_end(v);
258    C_do_apply(c - 1, ext_proc, k);
259  } else {
260    /* TODO: Convert to barf(), add new error code */
261    fprintf(stderr, "No extended number hook for %s!\n", ext_proc_name);
262    abort();
263  }
264}
265
266static C_word
267init_tags(___scheme_value tagvec)
268{
269  tags = CHICKEN_new_gc_root();
270  CHICKEN_gc_root_set(tags, tagvec);
271  return C_SCHEME_UNDEFINED;
272}
273
274#define ratnum_type_tag C_block_item(CHICKEN_gc_root_ref(tags), RAT_TAG)
275#define compnum_type_tag C_block_item(CHICKEN_gc_root_ref(tags), COMP_TAG)
276
277C_regparm C_word C_fcall C_i_numbers_numberp(C_word x)
278{
279  return C_mk_bool((x & C_FIXNUM_BIT) ||
280                   (!C_immediatep(x) && 
281                    (C_block_header(x) == C_FLONUM_TAG ||
282                     C_IS_BIGNUM_TYPE(x) ||
283                     (C_header_bits(x) == C_STRUCTURE_TYPE &&
284                      (C_block_item(x, 0) == ratnum_type_tag ||
285                       C_block_item(x, 0) == compnum_type_tag)))));
286}
287
288/* TODO: Rename to  C_i_integerp */
289C_regparm C_word C_fcall C_i_numbers_integerp(C_word x)
290{
291  return C_mk_bool((x & C_FIXNUM_BIT) ||
292                   (!C_immediatep(x) && 
293                    ((C_block_header(x) == C_FLONUM_TAG &&
294                      C_truep(C_u_i_fpintegerp_fixed(x))) ||
295                     C_IS_BIGNUM_TYPE(x))));
296}
297
298C_inline C_word basic_eqvp(C_word x, C_word y)
299{
300  return (x == y ||
301
302          (!C_immediatep(x) && !C_immediatep(y) &&
303           C_block_header(x) == C_block_header(y) &&
304           
305           ((C_block_header(x) == C_FLONUM_TAG &&
306             C_flonum_magnitude(x) == C_flonum_magnitude(y)) ||
307           
308            (C_IS_BIGNUM_TYPE(x) && C_u_i_bignum_cmp(x, y) == C_fix(0)))));
309}
310
311C_regparm C_word C_fcall C_i_basic_positivep(C_word x)
312{
313  if (x & C_FIXNUM_BIT)
314    return C_u_i_fixnum_positivep(x);
315  else if (C_immediatep(x))
316    barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "positive?", x);
317  else if (C_block_header(x) == C_FLONUM_TAG)
318    return C_mk_bool(C_flonum_magnitude(x) > 0.0);
319  else if (C_IS_BIGNUM_TYPE(x))
320    return C_mk_nbool(C_bignum_negativep(x));
321  else if (C_header_bits(x) == C_STRUCTURE_TYPE &&
322           (C_block_item(x, 0) == ratnum_type_tag))
323    return C_u_i_integer_positivep(C_block_item(x, 1));
324  else if (C_header_bits(x) == C_STRUCTURE_TYPE &&
325           (C_block_item(x, 0) == compnum_type_tag))
326    barf(C_BAD_ARGUMENT_TYPE_NO_REAL_ERROR, "positive?", x);
327  else
328    barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "positive?", x);
329}
330
331C_regparm C_word C_fcall C_u_i_integer_positivep(C_word x)
332{
333  if (x & C_FIXNUM_BIT) return C_u_i_fixnum_positivep(x);
334  else return C_mk_nbool(C_bignum_negativep(x));
335}
336
337C_regparm C_word C_fcall C_i_basic_negativep(C_word x)
338{
339  if (x & C_FIXNUM_BIT)
340    return C_u_i_fixnum_negativep(x);
341  else if (C_immediatep(x))
342    barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "negative?", x);
343  else if (C_block_header(x) == C_FLONUM_TAG)
344    return C_mk_bool(C_flonum_magnitude(x) < 0.0);
345  else if (C_IS_BIGNUM_TYPE(x))
346    return C_mk_bool(C_bignum_negativep(x));
347  else if (C_header_bits(x) == C_STRUCTURE_TYPE &&
348           (C_block_item(x, 0) == ratnum_type_tag))
349    return C_u_i_integer_negativep(C_block_item(x, 1));
350  else if (C_header_bits(x) == C_STRUCTURE_TYPE &&
351           (C_block_item(x, 0) == compnum_type_tag))
352    barf(C_BAD_ARGUMENT_TYPE_NO_REAL_ERROR, "negative?", x);
353  else
354    barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "negative?", x);
355}
356
357C_regparm C_word C_fcall C_u_i_integer_negativep(C_word x)
358{
359  if (x & C_FIXNUM_BIT) return C_u_i_fixnum_negativep(x);
360  else return C_mk_bool(C_bignum_negativep(x));
361}
362
363/* TODO: Rename to C_i_eqvp */
364C_regparm C_word C_fcall
365C_i_numbers_eqvp(C_word x, C_word y)
366{
367  return
368    C_mk_bool(basic_eqvp(x, y) ||
369              (!C_immediatep(x) && !C_immediatep(y) &&
370               (C_block_header(x) == C_block_header(y) &&
371                C_header_bits(x) == C_STRUCTURE_TYPE &&
372                C_block_item(x, 0) == C_block_item(y, 0) &&
373                (C_block_item(x, 0) == ratnum_type_tag ||
374                 C_block_item(x, 0) == compnum_type_tag) &&
375                basic_eqvp(C_block_item(x, 1), C_block_item(y, 1)) &&
376                basic_eqvp(C_block_item(x, 2), C_block_item(y, 2)))));
377}
378
379C_regparm C_word C_fcall C_i_nanp(C_word x)
380{
381  if (x & C_FIXNUM_BIT) {
382    return C_SCHEME_FALSE;
383  } else if (C_immediatep(x)) {
384    barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "nan?", x);
385  } else if (C_block_header(x) == C_FLONUM_TAG) {
386    return C_u_i_flonum_nanp(x);
387  } else if (C_IS_BIGNUM_TYPE(x)) {
388    return C_SCHEME_FALSE;
389  } else if (C_header_bits(x) == C_STRUCTURE_TYPE) {
390    /* To make this inlineable we don't call try_extended_number */
391    if (C_block_item(x, 0) == ratnum_type_tag)
392      return C_SCHEME_FALSE;
393    else if (C_block_item(x, 0) == compnum_type_tag)
394      return C_mk_bool(C_truep(C_i_nanp(C_block_item(x, 1))) ||
395                       C_truep(C_i_nanp(C_block_item(x, 2))));
396    else
397      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "nan?", x);
398  } else {
399    barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "nan?", x);
400  }
401}
402
403/* TODO: Rename to C_i_finitep */
404C_regparm C_word C_fcall C_i_numbers_finitep(C_word x)
405{
406  if (x & C_FIXNUM_BIT) {
407    return C_SCHEME_TRUE;
408  } else if (C_immediatep(x)) {
409    barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "finite?", x);
410  } else if (C_block_header(x) == C_FLONUM_TAG) {
411    return C_u_i_flonum_finitep(x);
412  } else if (C_IS_BIGNUM_TYPE(x)) {
413    return C_SCHEME_TRUE;
414  } else if (C_header_bits(x) == C_STRUCTURE_TYPE) {
415    /* To make this inlineable we don't use try_extended_number */
416    if (C_block_item(x, 0) == ratnum_type_tag)
417      return C_SCHEME_TRUE;
418    else if (C_block_item(x, 0) == compnum_type_tag)
419      return C_and(C_i_numbers_finitep(C_block_item(x, 1)),
420                   C_i_numbers_finitep(C_block_item(x, 2)));
421    else
422      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "finite?", x);
423  } else {
424    barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "finite?", x);
425  }
426}
427
428/* TODO: Rename to C_i_infinitep */
429C_regparm C_word C_fcall C_i_numbers_infinitep(C_word x)
430{
431  if (x & C_FIXNUM_BIT) {
432    return C_SCHEME_FALSE;
433  } else if (C_immediatep(x)) {
434    barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "infinite?", x);
435  } else if (C_block_header(x) == C_FLONUM_TAG) {
436    return C_u_i_flonum_infinitep(x);
437  } else if (C_IS_BIGNUM_TYPE(x)) {
438    return C_SCHEME_FALSE;
439  } else if (C_header_bits(x) == C_STRUCTURE_TYPE) {
440    /* To make this inlineable we don't use try_extended_number */
441    if (C_block_item(x, 0) == ratnum_type_tag)
442      return C_SCHEME_FALSE;
443    else if (C_block_item(x, 0) == compnum_type_tag)
444      return C_mk_bool(C_truep(C_i_numbers_infinitep(C_block_item(x, 1))) ||
445                       C_truep(C_i_numbers_infinitep(C_block_item(x, 2))));
446    else
447      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "infinite?", x);
448  } else {
449    barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "infinite?", x);
450  }
451}
452
453/* TODO: Rename to C_i_zerop */
454C_regparm C_word C_fcall C_i_numbers_zerop(C_word x)
455{
456  if (x & C_FIXNUM_BIT) {
457    return C_mk_bool(x == C_fix(0));
458  } else if (C_immediatep(x)) {
459    barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "zero?", x);
460  } else if (C_block_header(x) == C_FLONUM_TAG) {
461    return C_mk_bool(C_flonum_magnitude(x) == 0.0);
462  } else if (C_IS_BIGNUM_TYPE(x) ||
463             (C_header_bits(x) == C_STRUCTURE_TYPE &&
464              (C_block_item(x, 0) == ratnum_type_tag ||
465               C_block_item(x, 0) == compnum_type_tag))) {
466    return C_SCHEME_FALSE;
467  } else {
468    barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "zero?", x);
469  }
470}
471
472/* Copy all the digits from source to target, obliterating what was
473 * there.  If target is larger than source, the most significant
474 * digits will remain untouched.
475 */
476C_inline void bignum_digits_destructive_copy(C_word target, C_word source)
477{
478  C_memcpy(C_bignum_digits(target), C_bignum_digits(source),
479           C_wordstobytes(C_bignum_size(source)));
480}
481
482void C_ccall
483C_2_basic_plus(C_word c, C_word self, C_word k, C_word x, C_word y)
484{
485  if (x & C_FIXNUM_BIT) {
486    if (y & C_FIXNUM_BIT) {
487      C_word *a = C_alloc(C_SIZEOF_FIX_BIGNUM);
488      C_kontinue(k, C_a_u_i_2_fixnum_plus(&a, 2, x, y));
489    } else if (C_immediatep(y)) {
490      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "+", y);
491    } else if (C_block_header(y) == C_FLONUM_TAG) {
492      C_word *a = C_alloc(C_SIZEOF_FLONUM);
493      C_kontinue(k, C_flonum(&a, (double)C_unfix(x) + C_flonum_magnitude(y)));
494    } else if (C_IS_BIGNUM_TYPE(y)) {
495      C_u_2_integer_plus(4, (C_word)NULL, k, x, y);
496    } else {
497      try_extended_number("numbers#@extended-2-plus", 3, k, x, y);
498    }
499  } else if (C_immediatep(x)) {
500    barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "+", x);
501  } else if (C_block_header(x) == C_FLONUM_TAG) {
502    C_word *a = C_alloc(C_SIZEOF_FLONUM);
503    if (y & C_FIXNUM_BIT) {
504      C_kontinue(k, C_flonum(&a, C_flonum_magnitude(x) + (double)C_unfix(y)));
505    } else if (C_immediatep(y)) {
506      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "+", y);
507    } else if (C_block_header(y) == C_FLONUM_TAG) {
508      C_kontinue(k, C_a_i_flonum_plus(&a, 2, x, y));
509    } else if (C_IS_BIGNUM_TYPE(y)) {
510      C_kontinue(k, C_flonum(&a, C_flonum_magnitude(x)+C_bignum_to_double(y)));
511    } else {
512      try_extended_number("numbers#@extended-2-plus", 3, k, x, y);
513    }
514  } else if (C_IS_BIGNUM_TYPE(x)) {
515    if (y & C_FIXNUM_BIT) {
516      C_u_2_integer_plus(4, (C_word)NULL, k, x, y);
517    } else if (C_immediatep(y)) {
518      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "+", y);
519    } else if (C_block_header(y) == C_FLONUM_TAG) {
520      C_word *a = C_alloc(C_SIZEOF_FLONUM);
521      C_kontinue(k, C_flonum(&a, C_bignum_to_double(x)+C_flonum_magnitude(y)));
522    } else if (C_IS_BIGNUM_TYPE(y)) {
523      C_u_2_integer_plus(4, (C_word)NULL, k, x, y);
524    } else {
525      try_extended_number("numbers#@extended-2-plus", 3, k, x, y);
526    }
527  } else {
528    try_extended_number("numbers#@extended-2-plus", 3, k, x, y);
529  }
530}
531
532void C_ccall
533C_u_2_integer_plus(C_word c, C_word self, C_word k, C_word x, C_word y)
534{
535  if ((x & y) & C_FIXNUM_BIT) {
536    C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab;
537    C_kontinue(k, C_a_u_i_2_fixnum_plus(&a, 2, x, y));
538  } else {
539    C_word ab[C_SIZEOF_FIX_BIGNUM * 2], *a = ab;
540    if (x & C_FIXNUM_BIT) x = C_a_u_i_fix_to_big(&a, x);
541    if (y & C_FIXNUM_BIT) y = C_a_u_i_fix_to_big(&a, y);
542
543    if (C_bignum_negativep(x)) {
544      if (C_bignum_negativep(y)) bignum_plus_unsigned(k, x, y, C_SCHEME_TRUE);
545      else bignum_minus_unsigned(k, y, x);
546    } else {
547      if (C_bignum_negativep(y)) bignum_minus_unsigned(k, x, y);
548      else bignum_plus_unsigned(k, x, y, C_SCHEME_FALSE);
549    }
550  }
551}
552
553/* Needs C_SIZEOF_FIX_BIGNUM */
554C_regparm C_word C_fcall
555C_a_u_i_2_fixnum_plus(C_word **ptr, C_word n, C_word x, C_word y)
556{
557  /* Exceptional situation: this will cause a real underflow */
558  if (x == C_fix(C_MOST_NEGATIVE_FIXNUM) && y == C_fix(C_MOST_NEGATIVE_FIXNUM)) {
559    return C_bignum1(ptr, 1, ((C_uword)-C_MOST_NEGATIVE_FIXNUM) << 1);
560  } else {
561    C_word z = C_unfix(x) + C_unfix(y);
562
563    if(!C_fitsinfixnump(z)) {
564      /* TODO: function/macro returning either fixnum or bignum from a C int */
565      /* This should help with the C API/FFI too. */
566      return C_bignum1(ptr, z < 0, labs(z));
567    } else {
568      return C_fix(z);
569    }
570  }
571}
572
573void C_ccall
574C_2_basic_minus(C_word c, C_word self, C_word k, C_word x, C_word y)
575{
576  if (x & C_FIXNUM_BIT) {
577    if (y & C_FIXNUM_BIT) {
578      C_word *a = C_alloc(C_SIZEOF_FIX_BIGNUM);
579      C_kontinue(k, C_a_u_i_2_fixnum_minus(&a, 2, x, y));
580    } else if (C_immediatep(y)) {
581      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "-", y);
582    } else if (C_block_header(y) == C_FLONUM_TAG) {
583      C_word *a = C_alloc(C_SIZEOF_FLONUM);
584      C_kontinue(k, C_flonum(&a, (double)C_unfix(x) - C_flonum_magnitude(y)));
585    } else if (C_IS_BIGNUM_TYPE(y)) {
586      C_u_2_integer_minus(4, (C_word)NULL, k, x, y);
587    } else {
588      try_extended_number("numbers#@extended-2-minus", 3, k, x, y);
589    }
590  } else if (C_immediatep(x)) {
591    barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "-", x);
592  } else if (C_block_header(x) == C_FLONUM_TAG) {
593    C_word *a = C_alloc(C_SIZEOF_FLONUM);
594    if (y & C_FIXNUM_BIT) {
595      C_kontinue(k, C_flonum(&a, C_flonum_magnitude(x) - (double)C_unfix(y)));
596    } else if (C_immediatep(y)) {
597      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "-", y);
598    } else if (C_block_header(y) == C_FLONUM_TAG) {
599      C_kontinue(k, C_a_i_flonum_difference(&a, 2, x, y)); /* XXX NAMING! */
600    } else if (C_IS_BIGNUM_TYPE(y)) {
601      C_kontinue(k, C_flonum(&a, C_flonum_magnitude(x)-C_bignum_to_double(y)));
602    } else {
603      try_extended_number("numbers#@extended-2-minus", 3, k, x, y);
604    }
605  } else if (C_IS_BIGNUM_TYPE(x)) {
606    if (y & C_FIXNUM_BIT) {
607      C_u_2_integer_minus(4, (C_word)NULL, k, x, y);
608    } else if (C_immediatep(y)) {
609      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "-", y);
610    } else if (C_block_header(y) == C_FLONUM_TAG) {
611      C_word *a = C_alloc(C_SIZEOF_FLONUM);
612      C_kontinue(k, C_flonum(&a, C_bignum_to_double(x)-C_flonum_magnitude(y)));
613    } else if (C_IS_BIGNUM_TYPE(y)) {
614      C_u_2_integer_minus(4, (C_word)NULL, k, x, y);
615    } else {
616      try_extended_number("numbers#@extended-2-minus", 3, k, x, y);
617    }
618  } else {
619    try_extended_number("numbers#@extended-2-minus", 3, k, x, y);
620  }
621}
622
623void C_ccall
624C_u_2_integer_minus(C_word c, C_word self, C_word k, C_word x, C_word y)
625{
626  if ((x & y) & C_FIXNUM_BIT) {
627    C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab;
628    C_kontinue(k, C_a_u_i_2_fixnum_minus(&a, 2, x, y));
629  } else {
630    C_word ab[C_SIZEOF_FIX_BIGNUM * 2], *a = ab;
631    if (x & C_FIXNUM_BIT) x = C_a_u_i_fix_to_big(&a, x);
632    if (y & C_FIXNUM_BIT) y = C_a_u_i_fix_to_big(&a, y);
633
634    if (C_bignum_negativep(x)) {
635      if (C_bignum_negativep(y)) bignum_minus_unsigned(k, y, x);
636      else bignum_plus_unsigned(k, x, y, C_SCHEME_TRUE);
637    } else {
638      if (C_bignum_negativep(y)) bignum_plus_unsigned(k, x, y, C_SCHEME_FALSE);
639      else bignum_minus_unsigned(k, x, y);
640    }
641  }
642}
643
644/* Needs C_SIZEOF_FIX_BIGNUM */
645C_regparm C_word C_fcall
646C_a_u_i_2_fixnum_minus(C_word **ptr, C_word n, C_word x, C_word y)
647{
648  C_word z = C_unfix(x) - C_unfix(y);
649
650  if(!C_fitsinfixnump(z)) {
651    /* TODO: function/macro returning either fixnum or bignum from a C int */
652    /* This should help with the C API/FFI too. */
653    return C_bignum1(ptr, z < 0, labs(z));
654  } else {
655    return C_fix(z);
656  }
657}
658
659static void
660bignum_plus_unsigned(C_word k, C_word x, C_word y, C_word negp)
661{
662  C_word kab[C_SIZEOF_CLOSURE(4)], *ka = kab, k2, size;
663
664  if (C_bignum_size(y) > C_bignum_size(x)) {  /* Ensure size(y) <= size(x) */
665    C_word z = x;
666    x = y;
667    y = z;
668  }
669
670  k2 = C_closure(&ka, 4, (C_word)bignum_plus_unsigned_2, k, x, y);
671 
672  size = C_fix(C_bignum_size(x) + 1); /* One more digit, for possible carry. */
673  C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
674}
675
676static void
677bignum_plus_unsigned_2(C_word c, C_word self, C_word result)
678{
679  C_word k = C_block_item(self, 1),
680         x = C_block_item(self, 2),
681         y = C_block_item(self, 3);
682  C_uword *scan_y = C_bignum_digits(y),
683          *end_y = scan_y + C_bignum_size(y),
684          *scan_r = C_bignum_digits(result),
685          *end_r = scan_r + C_bignum_size(result),
686          sum, digit;
687  int carry = 0;
688
689  /* Copy x into r so we can operate on two pointers, which is faster
690   * than three, and we can stop earlier after adding y.  It's slower
691   * if x and y have equal length.  On average it's slightly faster.
692   */
693  bignum_digits_destructive_copy(result, x);
694  *(end_r-1) = 0; /* Ensure most significant digit is initialised */
695
696  /* Move over x and y simultaneously, destructively adding digits w/ carry. */
697  while (scan_y < end_y) {
698    digit = *scan_r;
699    if (carry) {
700      sum = digit + *scan_y++ + 1;
701      carry = sum <= digit;
702    } else {
703      sum = digit + *scan_y++;
704      carry = sum < digit;
705    }
706    (*scan_r++) = sum;
707  }
708 
709  /* The end of y, the smaller number.  Propagate carry into the rest of x. */
710  while (carry) {
711    sum = (*scan_r) + 1;
712    carry = (sum == 0);
713    (*scan_r++) = sum;
714  }
715  assert(scan_r <= end_r);
716
717  C_kontinue(k, C_bignum_simplify(result));
718}
719
720static int
721bignum_cmp_unsigned(C_word x, C_word y)
722{
723  C_word xlen = C_bignum_size(x), ylen = C_bignum_size(y);
724
725  if (xlen < ylen) {
726    return -1;
727  } else if (xlen > ylen) {
728    return 1;
729  } else if (x == y) {
730    return 0;
731  } else {
732    C_uword *startx = C_bignum_digits(x),
733            *scanx = startx + xlen,
734            *scany = C_bignum_digits(y) + ylen;
735
736    while (startx < scanx) {
737      C_uword xdigit = (*--scanx), ydigit = (*--scany);
738      if (xdigit < ydigit)
739        return -1;
740      if (xdigit > ydigit)
741        return 1;
742    }
743    return 0;
744  }
745}
746
747static void
748bignum_minus_unsigned(C_word k, C_word x, C_word y)
749{
750  C_word kab[C_SIZEOF_CLOSURE(4)], *ka = kab, k2, size;
751
752  switch(bignum_cmp_unsigned(x, y)) {
753  case 0:             /* x = y, return 0 */
754    C_kontinue(k, C_fix(0));
755  case -1:            /* abs(x) < abs(y), return -(abs(y) - abs(x)) */
756    k2 = C_closure(&ka, 4, (C_word)bignum_minus_unsigned_2, k, y, x);
757   
758    size = C_fix(C_bignum_size(y)); /* Maximum size of result is length of y. */
759    C_allocate_bignum(5, (C_word)NULL, k2, size, C_SCHEME_TRUE, C_SCHEME_FALSE);
760  case 1:             /* abs(x) > abs(y), return abs(x) - abs(y) */
761  default:
762    k2 = C_closure(&ka, 4, (C_word)bignum_minus_unsigned_2, k, x, y);
763   
764    size = C_fix(C_bignum_size(x)); /* Maximum size of result is length of x. */
765    C_allocate_bignum(5, (C_word)NULL, k2, size, C_SCHEME_FALSE, C_SCHEME_FALSE);
766    break;
767  }
768}
769
770static void
771bignum_minus_unsigned_2(C_word c, C_word self, C_word result)
772{
773  C_word k = C_block_item(self, 1),
774         x = C_block_item(self, 2),
775         y = C_block_item(self, 3);
776  C_uword *scan_r = C_bignum_digits(result),
777          *end_r = scan_r + C_bignum_size(result),
778          *scan_y = C_bignum_digits(y),
779          *end_y = scan_y + C_bignum_size(y),
780          difference, digit;
781  int borrow = 0;
782
783  bignum_digits_destructive_copy(result, x); /* See bignum_plus_unsigned_2 */
784
785  /* Destructively subtract y's digits w/ borrow from and back into r. */
786  while (scan_y < end_y) {
787    digit = *scan_r;
788    if (borrow) {
789      difference = digit - *scan_y++ - 1;
790      borrow = difference >= digit;
791    } else {
792      difference = digit - *scan_y++;
793      borrow = difference > digit;
794    }
795    (*scan_r++) = difference;
796  }
797
798  /* The end of y, the smaller number.  Propagate borrow into the rest of x. */
799  while (borrow) {
800    digit = *scan_r;
801    difference = digit - borrow;
802    borrow = difference >= digit;
803    (*scan_r++) = difference;
804  }
805
806  assert(scan_r <= end_r);
807
808  C_kontinue(k, C_bignum_simplify(result));
809}
810
811C_word C_ccall
812C_u_i_2_fixnum_gcd(C_word x, C_word y)
813{
814   x = (x & C_INT_SIGN_BIT) ? -C_unfix(x) : C_unfix(x);
815   y = (y & C_INT_SIGN_BIT) ? -C_unfix(y) : C_unfix(y);
816   
817   while(y != 0) {
818     C_word r = x % y;
819     x = y;
820     y = r;
821   }
822   return C_fix(x);
823}
824
825C_word C_ccall
826C_a_u_i_2_flonum_gcd(C_word **p, C_word n, C_word x, C_word y)
827{
828   double xub, yub, r;
829
830   if (!C_truep(C_u_i_fpintegerp_fixed(x)))
831     barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, "gcd", x);
832   if (!C_truep(C_u_i_fpintegerp_fixed(y)))
833     barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, "gcd", y);
834
835   xub = C_flonum_magnitude(x);
836   yub = C_flonum_magnitude(y);
837
838   if (xub < 0.0) xub = -xub;
839   if (yub < 0.0) yub = -yub;
840   
841   while(yub != 0.0) {
842     r = fmod(xub, yub);
843     xub = yub;
844     yub = r;
845   }
846   return C_flonum(p, xub);
847}
848
849void C_ccall
850C_basic_abs(C_word c, C_word self, C_word k, C_word x)
851{
852  if (c != 3) {
853    C_bad_argc_2(c, 3, self);
854  } else if (x & C_FIXNUM_BIT) {
855    C_word *a = C_alloc(C_SIZEOF_FIX_BIGNUM);
856    C_kontinue(k, C_a_u_i_fixnum_abs(&a, 1, x));
857  } else if (C_immediatep(x)) {
858    barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "abs", x);
859  } else if (C_block_header(x) == C_FLONUM_TAG) {
860    C_word *a = C_alloc(C_SIZEOF_FLONUM);
861    C_kontinue(k, C_a_i_flonum_abs(&a, 1, x));
862  } else if (C_IS_BIGNUM_TYPE(x)) {
863    C_u_integer_abs(3, (C_word)NULL, k, x);
864  } else {
865    try_extended_number("numbers#@extended-abs", 2, k, x);
866  }
867}
868
869void C_ccall
870C_u_integer_abs(C_word c, C_word self, C_word k, C_word x)
871{
872  if (x & C_FIXNUM_BIT) {
873    C_word *a = C_alloc(C_SIZEOF_FIX_BIGNUM);
874    C_kontinue(k, C_a_u_i_fixnum_abs(&a, 1, x));
875  } else if (C_bignum_negativep(x)) {
876    C_u_integer_negate(3, (C_word)NULL, k, x);
877  } else {
878    C_kontinue(k, x);
879  }
880}
881
882void C_ccall
883C_basic_signum(C_word c, C_word self, C_word k, C_word x)
884{
885  if (c != 3) {
886    C_bad_argc_2(c, 3, self);
887  } else if (x & C_FIXNUM_BIT) {
888    C_kontinue(k, C_u_i_fixnum_signum(x));
889  } else if (C_immediatep(x)) {
890    barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "signum", x);
891  } else if (C_block_header(x) == C_FLONUM_TAG) {
892    C_word *a = C_alloc(C_SIZEOF_FLONUM);
893    C_kontinue(k, C_a_u_i_flonum_signum(&a, 1, x));
894  } else if (C_IS_BIGNUM_TYPE(x)) {
895    C_kontinue(k, C_bignum_negativep(x) ? C_fix(-1) : C_fix(1));
896  } else {
897    try_extended_number("numbers#@extended-signum", 2, k, x);
898  }
899}
900
901C_regparm C_word C_fcall C_u_i_integer_signum(C_word x)
902{
903  if (x & C_FIXNUM_BIT) return C_u_i_fixnum_signum(x);
904  else return (C_bignum_negativep(x) ? C_fix(-1) : C_fix(1));
905}
906
907C_regparm C_word C_fcall C_i_basic_evenp(C_word x)
908{
909  if(x & C_FIXNUM_BIT) {
910    return C_i_fixnumevenp(x);
911  } else if(C_immediatep(x)) {
912    barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, "even?", x);
913  } else if (C_block_header(x) == C_FLONUM_TAG) {
914    double val, dummy;
915    val = C_flonum_magnitude(x);
916    if(C_isnan(val) || C_isinf(val) || C_modf(val, &dummy) != 0.0)
917      barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, "even?", x);
918    else
919      return C_mk_bool(fmod(val, 2.0) == 0.0);
920  } else if (C_IS_BIGNUM_TYPE(x)) {
921    return C_mk_nbool(C_bignum_digits(x)[0] & 1);
922  } else { /* No need to try extended number */
923    barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, "even?", x);
924  }
925}
926
927C_regparm C_word C_fcall C_u_i_integer_evenp(C_word x)
928{
929  if (x & C_FIXNUM_BIT) return C_i_fixnumevenp(x);
930  return C_mk_nbool(C_bignum_digits(x)[0] & 1);
931}
932
933C_regparm C_word C_fcall C_i_basic_oddp(C_word x)
934{
935  if(x & C_FIXNUM_BIT) {
936    return C_i_fixnumoddp(x);
937  } else if(C_immediatep(x)) {
938    barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, "odd?", x);
939  } else if(C_block_header(x) == C_FLONUM_TAG) {
940    double val, dummy;
941    val = C_flonum_magnitude(x);
942    if(C_isnan(val) || C_isinf(val) || C_modf(val, &dummy) != 0.0)
943      barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, "odd?", x);
944    else
945      return C_mk_bool(fmod(val, 2.0) != 0.0);
946  } else if (C_IS_BIGNUM_TYPE(x)) {
947    return C_mk_bool(C_bignum_digits(x)[0] & 1);
948  } else {
949    barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, "odd?", x);
950  }
951}
952
953C_regparm C_word C_fcall C_u_i_integer_oddp(C_word x)
954{
955  if (x & C_FIXNUM_BIT) return C_i_fixnumoddp(x);
956  return C_mk_bool(C_bignum_digits(x)[0] & 1);
957}
958
959void C_ccall
960C_basic_negate(C_word c, C_word self, C_word k, C_word x)
961{
962  if (x & C_FIXNUM_BIT) {
963    C_word *a = C_alloc(C_SIZEOF_FIX_BIGNUM);
964    C_kontinue(k, C_a_u_i_fixnum_negate(&a, 1, x));
965  } else if (C_immediatep(x)) {
966    barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "-", x);
967  } else if (C_block_header(x) == C_FLONUM_TAG) {
968    C_word *a = C_alloc(C_SIZEOF_FLONUM);
969    C_kontinue(k, C_a_i_flonum_negate(&a, 1, x));
970  } else if (C_IS_BIGNUM_TYPE(x)) {
971    C_u_integer_negate(3, (C_word)NULL, k, x);
972  } else {
973    try_extended_number("numbers#@extended-negate", 2, k, x);
974  }
975}
976
977void C_ccall
978C_u_integer_negate(C_word c, C_word self, C_word k, C_word x)
979{
980  if (x & C_FIXNUM_BIT) {
981    C_word *a = C_alloc(C_SIZEOF_FIX_BIGNUM);
982    C_kontinue(k, C_a_u_i_fixnum_negate(&a, 1, x));
983  } else {
984    if (C_bignum_negated_fitsinfixnump(x)) {
985      C_kontinue(k, C_fix(C_MOST_NEGATIVE_FIXNUM));
986    } else {
987      C_word *ka, k2, negp = C_mk_nbool(C_bignum_negativep(x)),
988             size = C_fix(C_bignum_size(x));
989      ka = C_alloc(C_SIZEOF_CLOSURE(3));
990      k2 = C_closure(&ka, 3, (C_word)bignum_negate_2, k, x);
991      C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
992    }
993  }
994}
995
996static void
997bignum_negate_2(C_word c, C_word self, C_word new_big)
998{
999  C_word k = C_block_item(self, 1),
1000         old_big = C_block_item(self, 2);
1001  bignum_digits_destructive_copy(new_big, old_big);
1002  C_kontinue(k, C_bignum_simplify(new_big));
1003}
1004
1005C_regparm C_word C_fcall
1006C_a_u_i_fixnum_negate(C_word **ptr, C_word n, C_word x)
1007{
1008  /* Exceptional situation: this will cause an overflow to itself */
1009  if (x == C_fix(C_MOST_NEGATIVE_FIXNUM)) /* C_fitsinfixnump(x) */
1010    return C_bignum1(ptr, 0, -C_MOST_NEGATIVE_FIXNUM);
1011  else
1012    return C_fix(-C_unfix(x));
1013}
1014
1015/* TODO: Rename to C_nequalp */
1016void C_ccall C_numbers_nequalp(C_word c, C_word self, C_word k, ...)
1017{
1018  C_word x, y, result;
1019  va_list v;
1020
1021  if (c < 4) C_bad_argc_2(c, 4, self);
1022
1023  c -= 2; 
1024  va_start(v, k);
1025
1026  x = va_arg(v, C_word);
1027  while(--c) {
1028    y = va_arg(v, C_word);
1029    result = C_i_2_basic_equalp(x, y);
1030    if (result == C_SCHEME_FALSE) break;
1031  }
1032
1033  va_end(v);
1034  C_kontinue(k, result);
1035}
1036
1037/* Compare two numbers as ratnums.  Either may be rat-, fix- or bignums */
1038static C_word rat_cmp(C_word x, C_word y)
1039{
1040  C_word ab[C_SIZEOF_FIX_BIGNUM*4], *a = ab, x1, x2, y1, y2,
1041         s, t, ssize, tsize, result, negp;
1042  C_uword *scan;
1043
1044  /* Check for 1 or 0; if x or y is this, the other must be the ratnum */
1045  if (x == C_fix(0)) {        /* Only the sign of y1 matters */
1046    return basic_cmp(x, C_block_item(y, 1), "ratcmp", 0);
1047  } else if (x == C_fix(1)) { /* x1*y1 <> x2*y2 --> y2 <> y1 | x1/x2 = 1/1 */
1048    return basic_cmp(C_block_item(y, 2), C_block_item(y, 1), "ratcmp", 0);
1049  } else if (y == C_fix(0)) { /* Only the sign of x1 matters */
1050    return basic_cmp(C_block_item(x, 1), y, "ratcmp", 0);
1051  } else if (y == C_fix(1)) { /* x1*y1 <> x2*y2 --> x1 <> x2 | y1/y2 = 1/1 */
1052    return basic_cmp(C_block_item(x, 1), C_block_item(x, 2), "ratcmp", 0);
1053  }
1054
1055  /* Extract components x=x1/x2 and y=y1/y2 */
1056  if (x & C_FIXNUM_BIT || C_IS_BIGNUM_TYPE(x)) x1 = x, x2 = C_fix(1);
1057  else x1 = C_block_item(x, 1), x2 = C_block_item(x, 2);
1058
1059  if (y & C_FIXNUM_BIT || C_IS_BIGNUM_TYPE(y)) y1 = y, y2 = C_fix(1);
1060  else y1 = C_block_item(y, 1), y2 = C_block_item(y, 2);
1061
1062  /* We only want to deal with bignums (this is tricky enough) */
1063  if (x1 & C_FIXNUM_BIT) x1 = C_a_u_i_fix_to_big(&a, x1);
1064  if (x2 & C_FIXNUM_BIT) x2 = C_a_u_i_fix_to_big(&a, x2);
1065  if (y1 & C_FIXNUM_BIT) y1 = C_a_u_i_fix_to_big(&a, y1);
1066  if (y2 & C_FIXNUM_BIT) y2 = C_a_u_i_fix_to_big(&a, y2);
1067
1068  /* We multiply using schoolbook method, so this will be very slow in
1069   * extreme cases.  This is a tradeoff we make so that comparisons
1070   * are inlineable, which makes a big difference for the common case.
1071   */
1072  ssize = C_bignum_size(x1) + C_bignum_size(y2);
1073  negp = C_mk_bool(C_bignum_negativep(x1));
1074  s = allocate_tmp_bignum(C_fix(ssize), negp, C_SCHEME_TRUE);
1075  bignum_digits_multiply(x1, y2, s); /* Swap args if x1 < y2? */
1076
1077  tsize = C_bignum_size(y1) + C_bignum_size(x2);
1078  negp = C_mk_bool(C_bignum_negativep(y1));
1079  t = allocate_tmp_bignum(C_fix(tsize), negp, C_SCHEME_TRUE);
1080  bignum_digits_multiply(y1, x2, t); /* Swap args if y1 < x2? */
1081
1082  /* Shorten the numbers if needed */
1083  for (scan = C_bignum_digits(s)+ssize-1; *scan == 0; scan--) ssize--;
1084  C_bignum_mutate_size(s, ssize);
1085  for (scan = C_bignum_digits(t)+tsize-1; *scan == 0; scan--) tsize--;
1086  C_bignum_mutate_size(t, tsize);
1087
1088  result = C_u_i_bignum_cmp(s, t);
1089
1090  free_tmp_bignum(t);
1091  free_tmp_bignum(s);
1092  return result;
1093}
1094
1095/* The primitive comparison operator.  eqp should be 1 if we're only
1096 * interested in equality testing (can speed things up and in case of
1097 * compnums, equality checking is the only available operation).  This
1098 * may return #f, in case there is no answer (for NaNs) or as a quick
1099 * and dirty non-zero answer when eqp is true.  Ugly but effective :)
1100 */
1101static C_word basic_cmp(C_word x, C_word y, char *loc, int eqp)
1102{
1103  if (x & C_FIXNUM_BIT) {
1104    if (y & C_FIXNUM_BIT) {
1105      return C_fix((x < y) ? -1 : ((x > y) ? 1 : 0));
1106    } else if (C_immediatep(y)) {
1107      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y);
1108    } else if (C_block_header(y) == C_FLONUM_TAG) {
1109      return int_flo_cmp(x, y);
1110    } else if (C_IS_BIGNUM_TYPE(y)) {
1111      C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab;
1112      return C_u_i_bignum_cmp(C_a_u_i_fix_to_big(&a, x), y);
1113    } else if (C_header_bits(y) == C_STRUCTURE_TYPE &&
1114               C_block_item(y, 0) == ratnum_type_tag) {
1115      if (eqp) return C_SCHEME_FALSE;
1116      else return rat_cmp(x, y);
1117    } else if (C_block_item(y, 0) == compnum_type_tag) {
1118      if (eqp) return C_SCHEME_FALSE;
1119      else barf(C_BAD_ARGUMENT_TYPE_COMPLEX_NO_ORDERING_ERROR, loc, y);
1120    } else {
1121      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y);
1122    }
1123  } else if (C_immediatep(x)) {
1124    barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, x);
1125  } else if (C_block_header(x) == C_FLONUM_TAG) {
1126    if (y & C_FIXNUM_BIT) {
1127      return flo_int_cmp(x, y);
1128    } else if (C_immediatep(y)) {
1129      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y);
1130    } else if (C_block_header(y) == C_FLONUM_TAG) {
1131      double a = C_flonum_magnitude(x), b = C_flonum_magnitude(y);
1132      if (C_isnan(a) || C_isnan(b)) return C_SCHEME_FALSE; /* "mu" */
1133      else return C_fix((a < b) ? -1 : ((a > b) ? 1 : 0));
1134    } else if (C_IS_BIGNUM_TYPE(y)) {
1135      return flo_int_cmp(x, y);
1136    } else if (C_header_bits(y) == C_STRUCTURE_TYPE) {
1137      if (C_block_item(y, 0) == ratnum_type_tag) {
1138        return flo_rat_cmp(x, y);
1139      } else if (C_block_item(y, 0) == compnum_type_tag) {
1140        if (eqp) return C_SCHEME_FALSE;
1141        else barf(C_BAD_ARGUMENT_TYPE_COMPLEX_NO_ORDERING_ERROR, loc, y);
1142      } else {
1143        barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y);
1144      }
1145    } else {
1146      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y);
1147    }
1148  } else if (C_IS_BIGNUM_TYPE(x)) {
1149    if (y & C_FIXNUM_BIT) {
1150      C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab;
1151      return C_u_i_bignum_cmp(x, C_a_u_i_fix_to_big(&a, y));
1152    } else if (C_immediatep(y)) {
1153      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y);
1154    } else if (C_block_header(y) == C_FLONUM_TAG) {
1155      return int_flo_cmp(x, y);
1156    } else if (C_IS_BIGNUM_TYPE(y)) {
1157      return C_u_i_bignum_cmp(x, y);
1158    } else if (C_header_bits(y) == C_STRUCTURE_TYPE &&
1159               C_block_item(y, 0) == ratnum_type_tag) {
1160      if (eqp) return C_SCHEME_FALSE;
1161      else return rat_cmp(x, y);
1162    } else if (C_block_item(y, 0) == compnum_type_tag) {
1163      if (eqp) return C_SCHEME_FALSE;
1164      else barf(C_BAD_ARGUMENT_TYPE_COMPLEX_NO_ORDERING_ERROR, loc, y);
1165    } else {
1166      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y);
1167    }
1168  } else if (C_header_bits(x) == C_STRUCTURE_TYPE &&
1169             (C_block_item(x, 0) == ratnum_type_tag)) {
1170    if (y & C_FIXNUM_BIT) {
1171      if (eqp) return C_SCHEME_FALSE;
1172      else return rat_cmp(x, y);
1173    } else if (C_immediatep(y)) {
1174      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y);
1175    } else if (C_block_header(y) == C_FLONUM_TAG) {
1176      return rat_flo_cmp(x, y);
1177    } else if (C_IS_BIGNUM_TYPE(y)) {
1178      if (eqp) return C_SCHEME_FALSE;
1179      else return rat_cmp(x, y);
1180    } else if (C_header_bits(y) == C_STRUCTURE_TYPE &&
1181               (C_block_item(y, 0) == ratnum_type_tag)) {
1182      if (eqp) {
1183        return C_and(C_and(C_u_i_2_integer_equalp(C_block_item(x, 1),
1184                                                  C_block_item(y, 1)),
1185                           C_u_i_2_integer_equalp(C_block_item(x, 2),
1186                                                  C_block_item(y, 2))),
1187                     C_fix(0));
1188      } else {
1189        return rat_cmp(x, y);
1190      }
1191    } else if (C_header_bits(y) == C_STRUCTURE_TYPE &&
1192               (C_block_item(y, 0) == compnum_type_tag)) {
1193      if (eqp) return C_SCHEME_FALSE;
1194      else barf(C_BAD_ARGUMENT_TYPE_COMPLEX_NO_ORDERING_ERROR, loc, y);
1195    } else {
1196      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y);
1197    }
1198  } else if (C_header_bits(x) == C_STRUCTURE_TYPE &&
1199             (C_block_item(x, 0) == compnum_type_tag)) {
1200    if (!eqp) {
1201      barf(C_BAD_ARGUMENT_TYPE_COMPLEX_NO_ORDERING_ERROR, loc, x);
1202    } else if (y & C_FIXNUM_BIT) {
1203      return C_SCHEME_FALSE;
1204    } else if (C_immediatep(y)) {
1205      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y);
1206    } else if (C_block_header(y) == C_FLONUM_TAG ||
1207               C_IS_BIGNUM_TYPE(y) ||
1208               (C_header_bits(y) == C_STRUCTURE_TYPE &&
1209                (C_block_item(y, 0) == ratnum_type_tag))) {
1210      return C_SCHEME_FALSE;
1211    } else if (C_header_bits(y) == C_STRUCTURE_TYPE &&
1212               (C_block_item(y, 0) == compnum_type_tag)) {
1213      return C_and(C_and(C_i_2_basic_equalp(C_block_item(x, 1),
1214                                            C_block_item(y, 1)),
1215                         C_i_2_basic_equalp(C_block_item(x, 2),
1216                                            C_block_item(y, 2))),
1217                   C_fix(0));
1218    } else {
1219      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y);
1220    }
1221  } else {
1222    barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, x);
1223  }
1224}
1225
1226C_regparm C_word C_fcall
1227C_i_2_basic_equalp(C_word x, C_word y)
1228{
1229   return C_mk_bool(basic_cmp(x, y, "=", 1) == C_fix(0));
1230}
1231
1232C_word C_ccall C_u_i_2_integer_equalp(C_word x, C_word y)
1233{
1234  if (x & C_FIXNUM_BIT)
1235    return C_mk_bool(x == y);
1236  else if (y & C_FIXNUM_BIT)
1237    return C_SCHEME_FALSE;
1238  else
1239    return C_mk_bool(C_u_i_bignum_cmp(x, y) == C_fix(0));
1240}
1241
1242/* TODO: Rename to C_lessp */
1243void C_ccall C_numbers_lessp(C_word c, C_word self, C_word k, ...)
1244{
1245  C_word x, y, result;
1246  va_list v;
1247
1248  if (c < 4) C_bad_argc_2(c, 4, self);
1249
1250  c -= 2; 
1251  va_start(v, k);
1252
1253  x = va_arg(v, C_word);
1254  while(--c) {
1255    y = va_arg(v, C_word);
1256    result = C_i_2_basic_lessp(x, y);
1257    if (result == C_SCHEME_FALSE) break;
1258    x = y;
1259  }
1260
1261  va_end(v);
1262  C_kontinue(k, result);
1263}
1264
1265/* TODO: Rename to C_less_or_equal_p */
1266void C_ccall C_numbers_less_or_equal_p(C_word c, C_word self, C_word k, ...)
1267{
1268  C_word x, y, result;
1269  va_list v;
1270
1271  if (c < 4) C_bad_argc_2(c, 4, self);
1272
1273  c -= 2; 
1274  va_start(v, k);
1275
1276  x = va_arg(v, C_word);
1277  while(--c) {
1278    y = va_arg(v, C_word);
1279    result = C_i_2_basic_less_or_equalp(x, y);
1280    if (result == C_SCHEME_FALSE) break;
1281    x = y;
1282  }
1283
1284  va_end(v);
1285  C_kontinue(k, result);
1286}
1287
1288C_regparm C_word C_fcall C_i_2_basic_lessp(C_word x, C_word y)
1289{
1290   return C_mk_bool(basic_cmp(x, y, "<", 0) == C_fix(-1));
1291}
1292
1293C_regparm C_word C_fcall C_i_2_basic_less_or_equalp(C_word x, C_word y)
1294{
1295   C_word res = basic_cmp(x, y, "<=", 0);
1296   return C_mk_bool(res == C_fix(0) || res == C_fix(-1));
1297}
1298
1299C_word C_ccall
1300C_u_i_2_integer_lessp(C_word x, C_word y)
1301{
1302  if (x & C_FIXNUM_BIT) {
1303    if (y & C_FIXNUM_BIT) {
1304      return C_mk_bool(C_unfix(x) < C_unfix(y));
1305    } else {
1306      return C_mk_nbool(C_bignum_negativep(y));
1307    }
1308  } else if (y & C_FIXNUM_BIT) {
1309    return C_mk_bool(C_bignum_negativep(x));
1310  } else {
1311    return C_mk_bool(C_u_i_bignum_cmp(x, y) == C_fix(-1));
1312  }
1313}
1314
1315C_word C_ccall
1316C_u_i_2_integer_less_or_equalp(C_word x, C_word y)
1317{
1318  if (x & C_FIXNUM_BIT) {
1319    if (y & C_FIXNUM_BIT) {
1320      return C_mk_bool(C_unfix(x) <= C_unfix(y));
1321    } else {
1322      return C_mk_nbool(C_bignum_negativep(y));
1323    }
1324  } else if (y & C_FIXNUM_BIT) {
1325    return C_mk_bool(C_bignum_negativep(x));
1326  } else {
1327    C_word res = C_u_i_bignum_cmp(x, y);
1328    return C_mk_bool(res == C_fix(0) || res == C_fix(-1));
1329  }
1330}
1331
1332/* TODO: Rename to C_greater_p */
1333void C_ccall C_numbers_greaterp(C_word c, C_word self, C_word k, ...)
1334{
1335  C_word x, y, result;
1336  va_list v;
1337
1338  if (c < 4) C_bad_argc_2(c, 4, self);
1339
1340  c -= 2; 
1341  va_start(v, k);
1342
1343  x = va_arg(v, C_word);
1344  while(--c) {
1345    y = va_arg(v, C_word);
1346    result = C_i_2_basic_greaterp(x, y);
1347    if (result == C_SCHEME_FALSE) break;
1348    x = y;
1349  }
1350
1351  va_end(v);
1352  C_kontinue(k, result);
1353}
1354
1355/* TODO: Rename to C_greater_or_equal_p */
1356void C_ccall C_numbers_greater_or_equal_p(C_word c, C_word self, C_word k, ...)
1357{
1358  C_word x, y, result;
1359  va_list v;
1360
1361  if (c < 4) C_bad_argc_2(c, 4, self);
1362
1363  c -= 2; 
1364  va_start(v, k);
1365
1366  x = va_arg(v, C_word);
1367  while(--c) {
1368    y = va_arg(v, C_word);
1369    result = C_i_2_basic_greater_or_equalp(x, y);
1370    if (result == C_SCHEME_FALSE) break;
1371    x = y;
1372  }
1373
1374  va_end(v);
1375  C_kontinue(k, result);
1376}
1377
1378C_regparm C_word C_fcall C_i_2_basic_greaterp(C_word x, C_word y)
1379{
1380   return C_mk_bool(basic_cmp(x, y, ">", 0) == C_fix(1));
1381}
1382
1383C_regparm C_word C_fcall C_i_2_basic_greater_or_equalp(C_word x, C_word y)
1384{
1385   C_word res = basic_cmp(x, y, ">=", 0);
1386   return C_mk_bool(res == C_fix(0) || res == C_fix(1));
1387}
1388
1389C_word C_ccall
1390C_u_i_2_integer_greaterp(C_word x, C_word y)
1391{
1392  if (x & C_FIXNUM_BIT) {
1393    if (y & C_FIXNUM_BIT) {
1394      return C_mk_bool(C_unfix(x) > C_unfix(y));
1395    } else {
1396      return C_mk_bool(C_bignum_negativep(y));
1397    }
1398  } else if (y & C_FIXNUM_BIT) {
1399    return C_mk_nbool(C_bignum_negativep(x));
1400  } else {
1401    return C_mk_bool(C_u_i_bignum_cmp(x, y) == C_fix(1));
1402  }
1403}
1404
1405C_word C_ccall
1406C_u_i_2_integer_greater_or_equalp(C_word x, C_word y)
1407{
1408  if (x & C_FIXNUM_BIT) {
1409    if (y & C_FIXNUM_BIT) {
1410      return C_mk_bool(C_unfix(x) >= C_unfix(y));
1411    } else {
1412      return C_mk_bool(C_bignum_negativep(y));
1413    }
1414  } else if (y & C_FIXNUM_BIT) {
1415    return C_mk_nbool(C_bignum_negativep(x));
1416  } else {
1417    C_word res = C_u_i_bignum_cmp(x, y);
1418    return C_mk_bool(res == C_fix(0) || res == C_fix(1));
1419  }
1420}
1421
1422/* This is a bit weird: We have to compare flonums as bignums due to
1423 * precision loss on 64-bit platforms.  For simplicity, we convert
1424 * fixnums to bignums here.
1425 */
1426static C_word int_flo_cmp(C_word intnum, C_word flonum)
1427{
1428  C_word ab[C_SIZEOF_FIX_BIGNUM + C_SIZEOF_FLONUM], *a = ab, x, y, res;
1429  double i, f;
1430
1431  f = C_flonum_magnitude(flonum);
1432
1433  if (C_isnan(f)) {
1434    return C_SCHEME_FALSE; /* "mu" */
1435  } else if (C_isinf(f)) {
1436    return C_fix((f > 0.0) ? -1 : 1); /* x is smaller if f is +inf.0 */
1437  } else {
1438    f = modf(f, &i);
1439
1440    x = (intnum & C_FIXNUM_BIT) ? C_a_u_i_fix_to_big(&a, intnum) : intnum;
1441    y = flo_to_tmp_bignum(C_flonum(&a, i));
1442
1443    res = C_u_i_bignum_cmp(x, y);
1444    free_tmp_bignum(y);
1445
1446    if (res == C_fix(0)) /* Use fraction to break tie. If f > 0, x is smaller */
1447      return C_fix((f > 0.0) ? -1 : ((f < 0.0) ? 1 : 0));
1448    else
1449      return res;
1450  }
1451}
1452
1453/* For convenience (ie, to reduce the degree of mindfuck) */
1454static C_word flo_int_cmp(C_word flonum, C_word intnum)
1455{
1456  C_word res = int_flo_cmp(intnum, flonum);
1457  switch(res) {
1458  case C_fix(1): return C_fix(-1);
1459  case C_fix(-1): return C_fix(1);
1460  default: return res; /* Can be either C_fix(0) or C_SCHEME_FALSE(!) */
1461  }
1462}
1463
1464/* This code is completely braindead, but at least it allows us to do
1465 * inline comparisons!
1466 */
1467static C_word rat_flo_cmp(C_word ratnum, C_word flonum)
1468{
1469  C_word ab[C_SIZEOF_FIX_BIGNUM * 2 + C_SIZEOF_FLONUM], *a = ab,
1470         num, denom, ibig, res, nscaled, iscaled, negp;
1471  C_uword *scan;
1472  int shift_amount, ilen, nlen;
1473  double i, f;
1474
1475  f = C_flonum_magnitude(flonum);
1476
1477  if (C_isnan(f)) {
1478    return C_SCHEME_FALSE; /* "mu" */
1479  } else if (C_isinf(f)) {
1480    return C_fix((f > 0.0) ? -1 : 1); /* x is smaller if f is +inf.0 */
1481  } else {
1482    /* Scale up the floating-point number to become a whole integer,
1483     * and remember power of two (# of bits) to shift the numerator.
1484     */
1485    shift_amount = 0;
1486
1487    /* TODO: This doesn't work for denormalized flonums! */
1488    while (modf(f, &i) != 0.0) {
1489      f = ldexp(f, 1);
1490      shift_amount++;
1491    }
1492
1493    i = f; /* TODO: split i and f so it'll work for denormalized flonums */
1494
1495    num = C_block_item(ratnum, 1);
1496    num = (num & C_FIXNUM_BIT) ? C_a_u_i_fix_to_big(&a, num) : num;
1497
1498    if (C_bignum_negativep(num) && i >= 0.0) { /* Save time if signs differ */
1499      return C_fix(-1);
1500    } else if (!C_bignum_negativep(num) && i <= 0.0) { /* num is never 0 */
1501      return C_fix(1);
1502    } else {
1503      negp = C_mk_bool(C_bignum_negativep(num));
1504
1505      denom = C_block_item(ratnum, 2);
1506      denom = (denom & C_FIXNUM_BIT) ? C_a_u_i_fix_to_big(&a, denom) : denom;
1507
1508      ibig = flo_to_tmp_bignum(C_flonum(&a, i));
1509
1510      nlen = C_bignum_size(num) + C_bignum_size(denom);
1511      ilen = C_bignum_size(ibig) + C_bignum_size(denom);
1512
1513      /* Now, multiply the scaled flonum by the denominator, so we can
1514       * compare it directly to the scaled numerator.  Unfortunately,
1515       * this won't use Karatsuba multiplication, so for large numbers
1516       * it will be slower than it could be if comparisons were done
1517       * in CPS context.
1518       */
1519      iscaled = allocate_tmp_bignum(C_fix(ilen), negp, C_SCHEME_TRUE);
1520      bignum_digits_multiply(denom, ibig, iscaled); /* Swap args if i < d? */
1521      free_tmp_bignum(ibig);
1522
1523      nlen += C_BIGNUM_BITS_TO_DIGITS(shift_amount);
1524      nscaled = allocate_tmp_bignum(C_fix(nlen), negp, C_SCHEME_TRUE);
1525
1526      scan = C_bignum_digits(nscaled) + shift_amount / C_BIGNUM_DIGIT_LENGTH;
1527      C_memcpy(scan, C_bignum_digits(num), C_wordstobytes(C_bignum_size(num)));
1528      shift_amount = shift_amount % C_BIGNUM_DIGIT_LENGTH;
1529      if(shift_amount > 0) {
1530        bignum_digits_destructive_shift_left(
1531         scan, C_bignum_digits(nscaled) + nlen, shift_amount);
1532      }
1533
1534      /* Shorten the numbers if needed */
1535      for (scan = C_bignum_digits(iscaled)+ilen-1; *scan == 0; scan--) ilen--;
1536      C_bignum_mutate_size(iscaled, ilen);
1537      for (scan = C_bignum_digits(nscaled)+nlen-1; *scan == 0; scan--) nlen--;
1538      C_bignum_mutate_size(nscaled, nlen);
1539
1540      /* Finally, we're ready to compare them! */
1541      res = C_u_i_bignum_cmp(nscaled, iscaled);
1542      free_tmp_bignum(nscaled);
1543      free_tmp_bignum(iscaled);
1544
1545      return res;
1546    }
1547  }
1548}
1549
1550static C_word flo_rat_cmp(C_word flonum, C_word ratnum)
1551{
1552  C_word res = rat_flo_cmp(ratnum, flonum);
1553  switch(res) {
1554  case C_fix(1): return C_fix(-1);
1555  case C_fix(-1): return C_fix(1);
1556  default: return res; /* Can be either C_fix(0) or C_SCHEME_FALSE(!) */
1557  }
1558}
1559
1560C_word
1561C_u_i_bignum_cmp(C_word x, C_word y)
1562{
1563  if (C_bignum_negativep(x)) {
1564    if (C_bignum_negativep(y)) { /* Largest negative number is smallest */
1565      return C_fix(bignum_cmp_unsigned(y, x));
1566    } else {
1567      return C_fix(-1);
1568    }
1569  } else {
1570    if (C_bignum_negativep(y)) {
1571      return C_fix(1);
1572    } else {
1573      return C_fix(bignum_cmp_unsigned(x, y));
1574    }
1575  }
1576}
1577
1578void C_ccall
1579C_allocate_bignum(C_word c, C_word self, C_word k, C_word size, C_word negp, C_word initp)
1580{
1581  C_word kab[C_SIZEOF_CLOSURE(3)], *ka = kab, k2, init;
1582  k2 = C_closure(&ka, 3, (C_word)allocate_bignum_2, k, negp);
1583
1584  init = C_and(initp, C_make_character('\0'));
1585  C_allocate_vector(6, (C_word)NULL, k2,
1586                    C_bytes(C_fixnum_plus(size, C_fix(1))), /* Add header */
1587                    /* Byte vec, initialization, align at 8 bytes */
1588                    C_SCHEME_TRUE, init, C_SCHEME_FALSE);
1589}
1590
1591static void
1592allocate_bignum_2(C_word c, C_word self, C_word bigvec)
1593{
1594  C_word ab[C_SIZEOF_STRUCTURE(2)], *a = ab, bignum,
1595         k = C_block_item(self, 1),
1596         negp = C_truep(C_block_item(self, 2)) ? 1 : 0,
1597         tagvec = CHICKEN_gc_root_ref(tags);
1598
1599  C_set_block_item(bigvec, 0, negp);
1600
1601  bignum = C_a_i_record2(&a, 2, C_block_item(tagvec, BIG_TAG), bigvec);
1602  C_kontinue(k, bignum);
1603}
1604
1605static C_word
1606allocate_tmp_bignum(C_word size, C_word negp, C_word initp)
1607{
1608  C_word *mem = malloc(C_wordstobytes(C_SIZEOF_BIGNUM(C_unfix(size)))),
1609          bigvec = (C_word)(mem + C_SIZEOF_STRUCTURE(2)),
1610          tagvec = CHICKEN_gc_root_ref(tags);
1611  if (mem == NULL) abort();     /* TODO: panic */
1612 
1613  C_block_header_init(bigvec, (C_STRING_TYPE | C_wordstobytes(C_unfix(size)+1)));
1614  C_set_block_item(bigvec, 0, C_truep(negp));
1615
1616  if (C_truep(initp)) {
1617    C_memset(((C_uword *)C_data_pointer(bigvec))+1,
1618             0, C_wordstobytes(C_unfix(size)));
1619  }
1620
1621  return C_a_i_record2(&mem, 2, C_block_item(tagvec, BIG_TAG), bigvec);
1622}
1623
1624/* Simplification: scan trailing zeroes, then return a fixnum if the
1625 * value fits, or trim the bignum's length. */
1626C_word C_ccall
1627C_bignum_simplify(C_word big)
1628{
1629  C_uword *start = C_bignum_digits(big),
1630          *last_digit = start + C_bignum_size(big) - 1,
1631          *scan = last_digit, tmp;
1632  int length;
1633
1634  while (scan >= start && *scan == 0)
1635    scan--;
1636  length = scan - start + 1;
1637 
1638  switch(length) {
1639  case 0:
1640    return C_fix(0);
1641  case 1:
1642    tmp = *start;
1643    if (C_bignum_negativep(big) ?
1644        !(tmp & C_INT_SIGN_BIT) && C_fitsinfixnump(-(C_word)tmp) :
1645        C_ufitsinfixnump(tmp))
1646      return C_bignum_negativep(big) ? C_fix(-(C_word)tmp) : C_fix(tmp);
1647    /* FALLTHROUGH */
1648  default:
1649    if (scan < last_digit) C_bignum_mutate_size(big, length);
1650    return big;
1651  }
1652}
1653
1654static C_uword
1655bignum_digits_destructive_scale_up_with_carry(C_uword *start, C_uword *end, C_uword factor, C_uword carry)
1656{
1657  C_uword digit, p;
1658
1659  assert(C_fitsinbignumhalfdigitp(carry));
1660  assert(C_fitsinbignumhalfdigitp(factor));
1661
1662  /* See fixnum_times.  Substitute xlo = factor, xhi = 0, y = digit
1663   * and simplify the result to reduce variable usage.
1664   */
1665  while (start < end) {
1666    digit = (*start);
1667
1668    p = factor * C_BIGNUM_DIGIT_LO_HALF(digit) + carry;
1669    carry = C_BIGNUM_DIGIT_LO_HALF(p);
1670
1671    p = factor * C_BIGNUM_DIGIT_HI_HALF(digit) + C_BIGNUM_DIGIT_HI_HALF(p);
1672    (*start++) = C_BIGNUM_DIGIT_COMBINE(C_BIGNUM_DIGIT_LO_HALF(p), carry);
1673    carry = C_BIGNUM_DIGIT_HI_HALF(p);
1674  }
1675  return carry;
1676}
1677
1678static C_uword
1679bignum_digits_destructive_scale_down(C_uword *start, C_uword *end, C_uword denominator)
1680{
1681  C_uword digit, k = 0;
1682  C_uhword q_j_hi, q_j_lo;
1683
1684  /* Single digit divisor case from Hacker's Delight, Figure 9-1,
1685   * adapted to modify u[] in-place instead of writing to q[].
1686   */
1687  while (start < end) {
1688    digit = (*--end);
1689
1690    k = C_BIGNUM_DIGIT_COMBINE(k, C_BIGNUM_DIGIT_HI_HALF(digit)); /* j */
1691    q_j_hi = k / denominator;
1692    k -= q_j_hi * denominator;
1693
1694    k = C_BIGNUM_DIGIT_COMBINE(k, C_BIGNUM_DIGIT_LO_HALF(digit)); /* j-1 */
1695    q_j_lo = k / denominator;
1696    k -= q_j_lo * denominator;
1697   
1698    *end = C_BIGNUM_DIGIT_COMBINE(q_j_hi, q_j_lo);
1699  }
1700  return k;
1701}
1702
1703static C_uword
1704bignum_digits_destructive_shift_right(C_uword *start, C_uword *end, int shift_right, int negp)
1705{
1706  int shift_left = C_BIGNUM_DIGIT_LENGTH - shift_right;
1707  C_uword digit, carry = negp ? ((~(C_uword)0) << shift_left) : 0;
1708
1709  assert(shift_right < C_BIGNUM_DIGIT_LENGTH);
1710
1711  while (start < end) {
1712    digit = *(--end);
1713    *end = (digit >> shift_right) | carry;
1714    carry = digit << shift_left;
1715  }
1716  return carry >> shift_left; /* The bits that were shifted out to the right */
1717}
1718
1719static C_uword
1720bignum_digits_destructive_shift_left(C_uword *start, C_uword *end, int shift_left)
1721{
1722  C_uword carry = 0, digit;
1723  int shift_right = C_BIGNUM_DIGIT_LENGTH - shift_left;
1724
1725  assert(shift_left < C_BIGNUM_DIGIT_LENGTH);
1726
1727  while (start < end) {
1728    digit = *start;
1729    (*start++) = (digit << shift_left) | carry;
1730    carry = digit >> shift_right;
1731  }
1732  return carry;  /* This would end up as most significant digit if it fit */
1733}
1734
1735void C_ccall
1736C_2_basic_times(C_word c, C_word self, C_word k, C_word x, C_word y)
1737{
1738  if (x & C_FIXNUM_BIT) {
1739    if (y & C_FIXNUM_BIT) {
1740      C_word *a = C_alloc(C_SIZEOF_BIGNUM(2));
1741      C_kontinue(k, C_a_u_i_2_fixnum_times(&a, 2, x, y));
1742    } else if (C_immediatep(y)) {
1743      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", y);
1744    } else if (C_block_header(y) == C_FLONUM_TAG) {
1745      C_word *a = C_alloc(C_SIZEOF_FLONUM);
1746      C_kontinue(k, C_flonum(&a, (double)C_unfix(x) * C_flonum_magnitude(y)));
1747    } else if (C_IS_BIGNUM_TYPE(y)) {
1748      C_u_2_integer_times(4, (C_word)NULL, k, x, y);
1749    } else {
1750      try_extended_number("numbers#@extended-2-times", 3, k, x, y);
1751    }
1752  } else if (C_immediatep(x)) {
1753    barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", x);
1754  } else if (C_block_header(x) == C_FLONUM_TAG) {
1755    C_word *a = C_alloc(C_SIZEOF_FLONUM);
1756    if (y & C_FIXNUM_BIT) {
1757      C_kontinue(k, C_flonum(&a, C_flonum_magnitude(x) * (double)C_unfix(y)));
1758    } else if (C_immediatep(y)) {
1759      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", y);
1760    } else if (C_block_header(y) == C_FLONUM_TAG) {
1761      C_kontinue(k, C_a_i_flonum_times(&a, 2, x, y));
1762    } else if (C_IS_BIGNUM_TYPE(y)) {
1763      C_kontinue(k, C_flonum(&a, C_flonum_magnitude(x)*C_bignum_to_double(y)));
1764    } else {
1765      try_extended_number("numbers#@extended-2-times", 3, k, x, y);
1766    }
1767  } else if (C_IS_BIGNUM_TYPE(x)) {
1768    if (y & C_FIXNUM_BIT) {
1769      C_u_2_integer_times(4, (C_word)NULL, k, x, y);
1770    } else if (C_immediatep(y)) {
1771      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", x);
1772    } else if (C_block_header(y) == C_FLONUM_TAG) {
1773      C_word *a = C_alloc(C_SIZEOF_FLONUM);
1774      C_kontinue(k, C_flonum(&a, C_bignum_to_double(x)*C_flonum_magnitude(y)));
1775    } else if (C_IS_BIGNUM_TYPE(y)) {
1776      C_u_2_integer_times(4, (C_word)NULL, k, x, y);
1777    } else {
1778      try_extended_number("numbers#@extended-2-times", 3, k, x, y);
1779    }
1780  } else {
1781    try_extended_number("numbers#@extended-2-times", 3, k, x, y);
1782  }
1783}
1784
1785/* Needs SIZEOF_BIGNUM(2) */
1786C_regparm C_word C_fcall
1787C_a_u_i_2_fixnum_times(C_word **ptr, C_word n, C_word x, C_word y)
1788{
1789  C_uword negp, xhi, xlo, yhi, ylo, p, rhi, rlo;
1790
1791  negp = ((x & C_INT_SIGN_BIT) ? !(y & C_INT_SIGN_BIT) : (y & C_INT_SIGN_BIT));
1792  x = (x & C_INT_SIGN_BIT) ? -C_unfix(x) : C_unfix(x);
1793  y = (y & C_INT_SIGN_BIT) ? -C_unfix(y) : C_unfix(y);
1794
1795  xhi = C_BIGNUM_DIGIT_HI_HALF(x); xlo = C_BIGNUM_DIGIT_LO_HALF(x);
1796  yhi = C_BIGNUM_DIGIT_HI_HALF(y); ylo = C_BIGNUM_DIGIT_LO_HALF(y);
1797 
1798  /* This is simply bignum_digits_multiply unrolled for 2x2 halfdigits */
1799  p = xlo * ylo;
1800  rlo = C_BIGNUM_DIGIT_LO_HALF(p);
1801
1802  p = xhi * ylo + C_BIGNUM_DIGIT_HI_HALF(p);
1803  rhi = C_BIGNUM_DIGIT_HI_HALF(p);
1804
1805  p = xlo * yhi + C_BIGNUM_DIGIT_LO_HALF(p);
1806  rlo = C_BIGNUM_DIGIT_COMBINE(C_BIGNUM_DIGIT_LO_HALF(p), rlo);
1807
1808  rhi = xhi * yhi + C_BIGNUM_DIGIT_HI_HALF(p) + rhi;
1809
1810  if (rhi) {
1811    return C_bignum2(ptr, negp != 0, rlo, rhi);
1812  } else if (negp ?
1813             ((rlo & C_INT_SIGN_BIT) || !C_fitsinfixnump(-(C_word)rlo)) :
1814             !C_ufitsinfixnump(rlo)) {
1815    return C_bignum1(ptr, negp != 0, rlo);
1816  } else {
1817    return C_fix(negp ? -rlo : rlo);
1818  }
1819}
1820
1821void C_ccall
1822C_u_2_integer_times(C_word c, C_word self, C_word k, C_word x, C_word y)
1823{
1824  if (x & C_FIXNUM_BIT) {
1825    if (y & C_FIXNUM_BIT) {
1826      C_word *a = C_alloc(C_SIZEOF_BIGNUM(2));
1827      C_kontinue(k, C_a_u_i_2_fixnum_times(&a, 2, x, y));
1828    } else {
1829      C_word tmp = x; /* swap to ensure x is a bignum and y a fixnum */
1830      x = y;
1831      y = tmp;
1832    }
1833  }
1834  /* Here, we know for sure that X is a bignum */
1835  if (y == C_fix(0)) {
1836    C_kontinue(k, C_fix(0));
1837  } else if (y == C_fix(1)) {
1838    C_kontinue(k, x);
1839  } else if (y == C_fix(-1)) {
1840    C_u_integer_negate(3, (C_word)NULL, k, x);
1841  } else if (y & C_FIXNUM_BIT) { /* Any other fixnum */
1842    C_word absy = (y & C_INT_SIGN_BIT) ? -C_unfix(y) : C_unfix(y),
1843           negp = C_mk_bool((y & C_INT_SIGN_BIT) ?
1844                            !C_bignum_negativep(x) :
1845                            C_bignum_negativep(x));
1846 
1847    if (C_fitsinbignumhalfdigitp(absy) ||
1848        (((C_uword)1 << (C_ilen(absy)-1)) == absy && C_fitsinfixnump(absy))) {
1849      C_word size, k2, *a = C_alloc(C_SIZEOF_CLOSURE(4));
1850      k2 = C_closure(&a, 4, (C_word)integer_times_2, k, x, C_fix(absy));
1851      size = C_fix(C_bignum_size(x) + 1); /* Needs _at most_ one more digit */
1852      C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
1853    } else {
1854      C_word *a = C_alloc(C_SIZEOF_FIX_BIGNUM);
1855      y = C_a_u_i_fix_to_big(&a, y);
1856      bignum_times_bignum_unsigned(k, x, y, negp);
1857    }
1858  } else {
1859    C_word negp = C_bignum_negativep(x) ?
1860                  !C_bignum_negativep(y) :
1861                  C_bignum_negativep(y);
1862    bignum_times_bignum_unsigned(k, x, y, C_mk_bool(negp));
1863  }
1864}
1865
1866static void
1867integer_times_2(C_word c, C_word self, C_word new_big)
1868{
1869  C_word k = C_block_item(self, 1),
1870         old_bigx = C_block_item(self, 2),
1871         absy = C_unfix(C_block_item(self, 3));
1872  C_uword *digits = C_bignum_digits(new_big),
1873          *end_digit = digits + C_bignum_size(old_bigx);
1874  int shift;
1875
1876  bignum_digits_destructive_copy(new_big, old_bigx);
1877
1878  /* Scale up, and sanitise the result. */
1879  shift = C_ilen(absy) - 1;
1880  if (((C_uword)1 << shift) == absy) { /* Power of two? */
1881    *end_digit = bignum_digits_destructive_shift_left(digits, end_digit, shift);
1882  } else {
1883    *end_digit =
1884      bignum_digits_destructive_scale_up_with_carry(digits, end_digit, absy, 0);
1885  }
1886  C_kontinue(k, C_bignum_simplify(new_big));
1887}
1888
1889static void
1890bignum_times_bignum_unsigned(C_word k, C_word x, C_word y, C_word negp)
1891{
1892  if (C_bignum_size(y) < C_bignum_size(x)) { /* Ensure size(x) <= size(y) */
1893    C_word z = x;
1894    x = y;
1895    y = z;
1896  }
1897
1898  if (C_bignum_size(x) < C_KARATSUBA_THRESHOLD) {  /* Gradebook */
1899    C_word kab[C_SIZEOF_CLOSURE(4)], *ka = kab, k2, size;
1900    k2 = C_closure(&ka, 4, (C_word)bignum_times_bignum_unsigned_2, k, x, y);
1901    size = C_fix(C_bignum_size(x) + C_bignum_size(y));
1902    C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_TRUE);
1903  } else {
1904    try_extended_number("numbers#@bignum-2-times-karatsuba", 3, k, x, y);
1905  }
1906}
1907
1908static C_regparm void
1909bignum_digits_multiply(C_word x, C_word y, C_word result)
1910{
1911  C_uword product,
1912          *xd = C_bignum_digits(x),
1913          *yd = C_bignum_digits(y),
1914          *rd = C_bignum_digits(result);
1915  C_uhword carry, yj;
1916  /* Lengths in halfwords */
1917  int i, j, length_x = C_bignum_size(x) * 2, length_y = C_bignum_size(y) * 2;
1918
1919  /* From Hacker's Delight, Figure 8-1 (top part) */
1920  for (j = 0; j < length_y; ++j) {
1921    yj = C_uhword_ref(yd, j);
1922    if (yj == 0) continue;
1923    carry = 0;
1924    for (i = 0; i < length_x; ++i) {
1925      product = (C_uword)C_uhword_ref(xd, i) * yj +
1926                (C_uword)C_uhword_ref(rd, i + j) + carry;
1927      C_uhword_set(rd, i + j, product);
1928      carry = C_BIGNUM_DIGIT_HI_HALF(product);
1929    }
1930    C_uhword_set(rd, j + length_x, carry);
1931  }
1932}
1933
1934static void
1935bignum_times_bignum_unsigned_2(C_word c, C_word self, C_word result)
1936{
1937  C_word k = C_block_item(self, 1),
1938         x = C_block_item(self, 2),
1939         y = C_block_item(self, 3);
1940
1941  bignum_digits_multiply(x, y, result);
1942  C_kontinue(k, C_bignum_simplify(result));
1943}
1944
1945void C_ccall
1946C_digits_to_integer(C_word c, C_word self, C_word k, C_word str,
1947                    C_word start, C_word end, C_word radix, C_word negp)
1948{
1949  assert((C_unfix(radix) > 1) && C_fitsinbignumhalfdigitp(C_unfix(radix)));
1950 
1951  if (start == end) {
1952    C_kontinue(k, C_SCHEME_FALSE);
1953  } else {
1954    C_word kab[C_SIZEOF_CLOSURE(6)], *ka = kab, k2, size;
1955    size_t nbits;
1956    k2 = C_closure(&ka, 6, (C_word)digits_to_integer_2, k, str, start, end, radix);
1957
1958    nbits = (C_unfix(end) - C_unfix(start)) * C_ilen(C_unfix(radix)-1);
1959    size = C_fix(C_BIGNUM_BITS_TO_DIGITS(nbits));
1960    C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
1961  }
1962}
1963
1964C_inline int hex_char_to_digit(int ch)
1965{
1966  if (ch == (int)'#') return 0; /* Hash characters in numbers are mapped to 0 */
1967  else if (ch >= (int)'a') return ch - (int)'a' + 10; /* lower hex */
1968  else if (ch >= (int)'A') return ch - (int)'A' + 10; /* upper hex */
1969  else return ch - (int)'0'; /* decimal (OR INVALID; handled elsewhere) */
1970}
1971
1972static void
1973digits_to_integer_2(C_word c, C_word self, C_word result)
1974{
1975  C_word k = C_block_item(self, 1),
1976         str = C_block_item(self, 2),
1977         start = C_unfix(C_block_item(self, 3)),
1978         end = C_unfix(C_block_item(self, 4)),
1979         radix = C_unfix(C_block_item(self, 5));
1980  char *s = C_c_string(str);
1981
1982  C_kontinue(k, str_to_bignum(result, s + start, s + end, radix));
1983}
1984
1985/* Write from digit character stream to bignum.  Bignum does not need
1986 * to be initialised.  Returns the bignum, or a fixnum.  Assumes the
1987 * string contains only digits that fit within radix (checked by
1988 * string->number).
1989 */
1990static C_regparm C_word
1991str_to_bignum(C_word bignum, char *str, char *str_end, int radix)
1992{
1993  int radix_shift, str_digit;
1994  C_uword *digits = C_bignum_digits(bignum),
1995          *end_digits = digits + C_bignum_size(bignum), big_digit = 0;
1996
1997  /* Below, we try to save up as much as possible in big_digit, and
1998   * only when it exceeds what we would be able to multiply easily, we
1999   * scale up the bignum and add what we saved up.
2000   */
2001  radix_shift = C_ilen(radix) - 1;
2002  if (((C_uword)1 << radix_shift) == radix) { /* Power of two? */
2003    int n = 0; /* Number of bits read so far into current big digit */
2004
2005    /* Read from least to most significant digit to avoid shifting or scaling */
2006    while (str_end > str) {
2007      str_digit = hex_char_to_digit((int)*--str_end);
2008
2009      big_digit |= (C_uword)str_digit << n;
2010      n += radix_shift;
2011
2012      if (n >= C_BIGNUM_DIGIT_LENGTH) {
2013        n -= C_BIGNUM_DIGIT_LENGTH;
2014        *digits++ = big_digit;
2015        big_digit = str_digit >> (radix_shift - n);
2016      }
2017    }
2018    assert(n < C_BIGNUM_DIGIT_LENGTH);
2019    /* If radix isn't an exact divisor of digit length, write final digit */
2020    if (n > 0) *digits++ = big_digit;
2021    assert(digits == end_digits);
2022  } else {                        /* Not a power of two */
2023    C_uword *last_digit = digits, factor;  /* bignum starts as zero */
2024
2025    do {
2026      factor = radix;
2027      while (str < str_end && C_fitsinbignumhalfdigitp(factor)) {
2028        str_digit = hex_char_to_digit((int)*str++);
2029        factor *= radix;
2030        big_digit = radix * big_digit + str_digit;
2031      }
2032
2033      big_digit = bignum_digits_destructive_scale_up_with_carry(
2034                   digits, last_digit, factor / radix, big_digit);
2035
2036      if (big_digit) {
2037        (*last_digit++) = big_digit; /* Move end */
2038        big_digit = 0;
2039      }
2040    } while (str < str_end);
2041
2042    /* Set remaining digits to zero so bignum_simplify can do its work */
2043    assert(last_digit <= end_digits);
2044    while (last_digit < end_digits) *last_digit++ = 0;
2045  }
2046
2047  return C_bignum_simplify(bignum);
2048}
2049
2050/* TODO: Copied from runtime.c */
2051# define STRING_BUFFER_SIZE   4096
2052
2053static C_TLS C_char buffer[ STRING_BUFFER_SIZE ];
2054static char *to_n_nary(C_uword num, C_uword base, int negp, int as_flonum)
2055{
2056  static char *digits = "0123456789abcdef";
2057  char *p;
2058  C_uword shift = C_ilen(base) - 1;
2059  int mask = (1 << shift) - 1;
2060  if (as_flonum) {
2061    buffer[68] = '\0';
2062    buffer[67] = '0';
2063    buffer[66] = '.';
2064  } else {
2065    buffer[66] = '\0';
2066  }
2067  p = buffer + 66;
2068  if (mask == base - 1) {
2069    do {
2070      *(--p) = digits [ num & mask ];
2071      num >>= shift;
2072    } while (num);
2073  } else {
2074    do {
2075      *(--p) = digits [ num % base ];
2076      num /= base;
2077    } while (num);
2078  }
2079  if (negp) *(--p) = '-';
2080  return p;
2081}
2082
2083void C_ccall C_basic_number_to_string(C_word c, C_word closure, C_word k, C_word num, ...)
2084{
2085  C_word radix;
2086
2087  if(c == 3) {
2088    radix = C_fix(10);
2089  } else if(c == 4) {
2090    va_list v;
2091
2092    va_start(v, num);
2093    radix = va_arg(v, C_word);
2094    va_end(v);
2095   
2096    if(!(radix & C_FIXNUM_BIT))
2097      barf(C_BAD_ARGUMENT_TYPE_BAD_BASE_ERROR, "number->string", radix);
2098  } else {
2099    C_bad_argc(c, 3);
2100  }
2101
2102  if(num & C_FIXNUM_BIT) {
2103    C_u_fixnum_to_string(4, (C_word)NULL, k, num, radix);
2104  } else if (C_immediatep(num)) {
2105    barf(C_BAD_ARGUMENT_TYPE_ERROR, "number->string", num);
2106  } else if(C_block_header(num) == C_FLONUM_TAG) {
2107    C_u_flonum_to_string(4, (C_word)NULL, k, num, radix);
2108  } else if (C_IS_BIGNUM_TYPE(num)) {
2109    C_u_integer_to_string(4, (C_word)NULL, k, num, radix);
2110  } else {
2111    try_extended_number("numbers#@extended-number->string", 3, k, num, radix);
2112  }
2113}
2114
2115void C_ccall
2116C_u_fixnum_to_string(C_word c, C_word self, C_word k, C_word num, C_word radix)
2117{
2118  C_char *p;
2119  C_word *a, neg = num & C_INT_SIGN_BIT ? 1 : 0;
2120
2121  radix = C_unfix(radix);
2122  if (radix < 2 || radix > 16) {
2123    barf(C_BAD_ARGUMENT_TYPE_BAD_BASE_ERROR, "number->string", radix);
2124  }
2125
2126  num = neg ? -C_unfix(num) : C_unfix(num);
2127  p = to_n_nary(num, radix, neg, 0);
2128
2129  num = C_strlen(p);
2130  a = C_alloc((C_bytestowords(num) + 1));
2131  C_kontinue(k, C_string(&a, num, p));
2132}
2133
2134void C_ccall
2135C_u_flonum_to_string(C_word c, C_word self, C_word k, C_word num, C_word radix)
2136{
2137  C_word *a;
2138  C_char *p;
2139  double f;
2140
2141  radix = C_unfix(radix);
2142  f = C_flonum_magnitude(num);
2143
2144  /* XXX TODO: Should inexacts be printable in other bases than 10?
2145   * Perhaps output a string starting with #i?
2146   * Right now something like (number->string 1e40 16) results in
2147   * a string that can't be read back using string->number.
2148   */
2149  if((radix < 2) || (radix > 16)){
2150    barf(C_BAD_ARGUMENT_TYPE_BAD_BASE_ERROR, "number->string", C_fix(radix));
2151  }
2152
2153  if(C_fits_in_unsigned_int_p(num) == C_SCHEME_TRUE) { /* Use fast int code */
2154    if(f < 0) {
2155      p = to_n_nary((C_uword)-f, radix, 1, 1);
2156    } else {
2157      p = to_n_nary((C_uword)f, radix, 0, 1);
2158    }
2159  } else if(C_isnan(f)) {
2160    p = "+nan.0";
2161  } else if(C_isinf(f)) {
2162    p = f > 0 ? "+inf.0" : "-inf.0";
2163  } else { /* Doesn't fit an unsigned int and not "special"; use system libc */
2164    C_snprintf(buffer, STRING_BUFFER_SIZE, C_text("%.*g"),
2165               /* XXX: flonum_print_precision */
2166               (int)C_unfix(C_get_print_precision()), f);
2167    buffer[STRING_BUFFER_SIZE-1] = '\0';
2168
2169    if((p = C_strpbrk(buffer, C_text(".eE"))) == NULL) {
2170      /* Already checked for these, so shouldn't happen */
2171      assert(*buffer != 'i'); /* "inf" */
2172      assert(*buffer != 'n'); /* "nan" */
2173      /* Ensure integral flonums w/o expt are always terminated by .0 */
2174#if defined(HAVE_STRLCAT) || !defined(C_strcat)
2175      C_strlcat(buffer, C_text(".0"), sizeof(buffer));
2176#else
2177      C_strcat(buffer, C_text(".0"));
2178#endif
2179    }
2180    p = buffer;
2181  }
2182
2183  radix = C_strlen(p);
2184  a = C_alloc((C_bytestowords(radix) + 1));
2185  radix = C_string(&a, radix, p);
2186  C_kontinue(k, radix);
2187}
2188
2189/* Naming is a little inconsistent, but looks saner.  We're not R-O-B-O-T-S! */
2190void C_ccall
2191C_u_integer_to_string(C_word c, C_word self, C_word k, C_word num, C_word radix)
2192{
2193  if (num & C_FIXNUM_BIT) {
2194    C_u_fixnum_to_string(4, (C_word)NULL, k, num, radix);
2195  } else {
2196    int len, radix_shift;
2197    size_t nbits;
2198
2199    if ((C_unfix(radix) < 2) || (C_unfix(radix) > 16)) {
2200      barf(C_BAD_ARGUMENT_TYPE_BAD_BASE_ERROR, "number->string", radix);
2201    }
2202
2203    /* Approximation of the number of radix digits we'll need.  We try
2204     * to be as precise as possible to avoid memmove overhead at the end
2205     * of the non-powers of two part of the conversion procedure, which
2206     * we may need to do because we write strings back-to-front, and
2207     * pointers must be aligned (even for byte blocks).
2208     */
2209    len = C_bignum_size(num)-1;
2210
2211    nbits  = (size_t)len * C_BIGNUM_DIGIT_LENGTH;
2212    nbits += C_ilen(C_bignum_digits(num)[len]);
2213
2214    len = C_ilen(C_unfix(radix))-1;
2215    len = (nbits + len - 1) / len;
2216    len += C_bignum_negativep(num) ? 1 : 0; /* Add space for negative sign */
2217   
2218    radix_shift = C_ilen(C_unfix(radix)) - 1;
2219    if (len > C_RECURSIVE_TO_STRING_THRESHOLD &&
2220        /* The power of two fast path is much faster than recursion */
2221        ((C_uword)1 << radix_shift) != C_unfix(radix)) {
2222      try_extended_number("numbers#@integer->string/recursive",
2223                          4, k, num, radix, C_fix(len));
2224    } else {
2225      C_word k2, *ka;
2226      ka = C_alloc(C_SIZEOF_CLOSURE(4));
2227      k2 = C_closure(&ka, 4, (C_word)bignum_to_str_2, k, num, radix);
2228      C_allocate_vector(6, (C_word)NULL, k2, C_fix(len),
2229                        /* Byte vec, no initialization, no align at 8 bytes */
2230                        C_SCHEME_TRUE, C_SCHEME_FALSE, C_SCHEME_FALSE);
2231    }
2232  }
2233}
2234
2235static void
2236bignum_to_str_2(C_word c, C_word self, C_word string)
2237{
2238  static char *characters = "0123456789abcdef";
2239  C_word k = C_block_item(self, 1),
2240         bignum = C_block_item(self, 2),
2241         radix = C_unfix(C_block_item(self, 3));
2242  char *buf = C_c_string(string), *index = buf + C_header_size(string) - 1;
2243  int radix_shift, negp = (C_bignum_negativep(bignum) ? 1 : 0);
2244
2245  radix_shift = C_ilen(radix) - 1;
2246  if (((C_uword)1 << radix_shift) == radix) { /* Power of two? */
2247    int radix_mask = radix - 1, big_digit_len = 0, radix_digit;
2248    C_uword *scan, *end, big_digit = 0;
2249
2250    scan = C_bignum_digits(bignum);
2251    end = scan + C_bignum_size(bignum);
2252
2253    while (scan < end) {
2254      /* If radix isn't an exact divisor of digit length, handle overlap */
2255      if (big_digit_len == 0) {
2256        big_digit = *scan++;
2257        big_digit_len = C_BIGNUM_DIGIT_LENGTH;
2258      } else {
2259        assert(index >= buf);
2260        radix_digit = big_digit;
2261        big_digit = *scan++;
2262        radix_digit |= (big_digit << big_digit_len) & radix_mask;
2263        big_digit >>= (radix_shift - big_digit_len);
2264        big_digit_len = C_BIGNUM_DIGIT_LENGTH - big_digit_len;
2265      }
2266
2267      while(big_digit_len >= radix_shift && index >= buf) {
2268        radix_digit = big_digit & radix_mask;
2269        *index-- = characters[radix_digit];
2270        big_digit >>= radix_shift;
2271        big_digit_len -= radix_shift;
2272      }
2273    }
2274
2275    assert(big_digit < radix);
2276
2277    /* Final digit (like overlap at start of while loop) */
2278    if (big_digit) *index-- = characters[big_digit];
2279
2280    if (negp) {
2281      /* Loop above might've overwritten sign position with a zero */
2282      if (*(index+1) == '0') *(index+1) = '-';
2283      else *index-- = '-';
2284    }
2285
2286    /* Length calculation is always precise for radix powers of two. */
2287    assert(index == buf-1);
2288  } else {
2289    C_uword base, *start, *scan, big_digit;
2290    C_word working_copy;
2291    int steps, i;
2292
2293    working_copy = allocate_tmp_bignum(C_fix(C_bignum_size(bignum)),
2294                                       C_mk_bool(negp), C_SCHEME_FALSE);
2295    bignum_digits_destructive_copy(working_copy, bignum);
2296
2297    start = C_bignum_digits(working_copy);
2298
2299    scan = start + C_bignum_size(bignum);
2300    /* Calculate the largest power of radix that fits a halfdigit:
2301     * steps = log10(2^halfdigit_bits), base = 10^steps
2302     */
2303    for(steps = 0, base = radix; C_fitsinbignumhalfdigitp(base); base *= radix)
2304      steps++;
2305
2306    base /= radix; /* Back down: we overshot in the loop */
2307
2308    while (scan > start) {
2309      big_digit = bignum_digits_destructive_scale_down(start, scan, base);
2310
2311      if (*(scan-1) == 0) scan--; /* Adjust if we exhausted the highest digit */
2312
2313      for(i = 0; i < steps && index >= buf; ++i) {
2314        C_word tmp = big_digit / radix;
2315        *index-- = characters[big_digit - (tmp*radix)]; /* big_digit % radix */
2316        big_digit = tmp;
2317      }
2318    }
2319    assert(index >= buf-1);
2320    free_tmp_bignum(working_copy);
2321
2322    /* Move index onto first nonzero digit.  We're writing a bignum
2323       here: it can't consist of only zeroes. */
2324    while(*++index == '0');
2325 
2326    if (negp) *--index = '-';
2327 
2328    /* Shorten with distance between start and index. */
2329    if (buf != index) {
2330      i = C_header_size(string) - (index - buf);
2331      C_memmove(buf, index, i); /* Move start of number to beginning. */
2332      C_block_header(string) = C_STRING_TYPE | i; /* Mutate strlength. */
2333    }
2334  }
2335
2336  C_kontinue(k, string);
2337}
2338
2339C_regparm double C_bignum_to_double(C_word bignum)
2340{
2341  double accumulator = 0;
2342  C_uword *start = C_bignum_digits(bignum),
2343          *scan = start + C_bignum_size(bignum);
2344  while (start < scan) {
2345    accumulator *= (C_word)1 << C_BIGNUM_HALF_DIGIT_LENGTH;
2346    accumulator *= (C_word)1 << C_BIGNUM_HALF_DIGIT_LENGTH;
2347    accumulator += (*--scan);
2348  }
2349  return(C_bignum_negativep(bignum) ? -accumulator : accumulator);
2350}
2351
2352static void
2353fabs_frexp_to_digits(C_uword exp, double sign, C_uword *start, C_uword *scan)
2354{
2355  C_uword digit, odd_bits = exp % C_BIGNUM_DIGIT_LENGTH;
2356
2357  assert(C_isfinite(sign));
2358  assert(0.5 <= sign && sign < 1); /* Guaranteed by frexp() and fabs() */
2359  assert((scan - start) == C_BIGNUM_BITS_TO_DIGITS(exp));
2360 
2361  if (odd_bits > 0) { /* Handle most significant digit first */
2362    sign *= (C_uword)1 << odd_bits;
2363    digit = (C_uword)sign;
2364    (*--scan) = digit;
2365    sign -= (double)digit;
2366  }
2367
2368  while (start < scan && sign > 0) {
2369    sign *= pow(2.0, C_BIGNUM_DIGIT_LENGTH);
2370    digit = (C_uword)sign;
2371    (*--scan) = digit;
2372    sign -= (double)digit;
2373  }
2374
2375  /* Finish up by clearing any remaining, lower, digits */
2376  while (start < scan)
2377    (*--scan) = 0;
2378}
2379
2380static C_word
2381flo_to_tmp_bignum(C_word x)
2382{
2383  /* TODO: allocating and initialising the bignum is pointless if we
2384   * already know the number of limbs in the comparand.  In fact,
2385   * bignum_cmp will first check the number of limbs and *then*
2386   * compare.  Instead, we can check beforehand and check the limbs
2387   * directly against the generated limbs, without allocating at all!
2388   */
2389  C_word tmp_big, negp = C_mk_bool(C_flonum_magnitude(x) < 0.0);
2390  int exponent;
2391  double significand = frexp(C_flonum_magnitude(x), &exponent);
2392
2393  assert(C_u_i_fpintegerp_fixed(x));
2394
2395  if (exponent <= 0) {
2396    tmp_big = allocate_tmp_bignum(C_fix(0), C_SCHEME_FALSE, C_SCHEME_FALSE);
2397  } else if (exponent == 1) { /* TODO: check significand * 2^exp fits fixnum? */
2398    /* Don't use fix_to_big to simplify caller code: it can just free this */
2399    tmp_big = allocate_tmp_bignum(C_fix(1), negp, C_SCHEME_FALSE);
2400    C_bignum_digits(tmp_big)[0] = 1;
2401  } else {
2402    C_uword size, *start, *end;
2403
2404    size = C_fix(C_BIGNUM_BITS_TO_DIGITS(exponent));
2405
2406    tmp_big = allocate_tmp_bignum(size, negp, C_SCHEME_FALSE);
2407    start = C_bignum_digits(tmp_big);
2408    end = start + C_bignum_size(tmp_big);
2409
2410    fabs_frexp_to_digits(exponent, fabs(significand), start, end);
2411  }
2412  return tmp_big;
2413}
2414
2415void C_ccall
2416C_u_flo_to_int(C_word c, C_word self, C_word k, C_word x)
2417{
2418  int exponent;
2419  double significand = frexp(C_flonum_magnitude(x), &exponent);
2420
2421  assert(C_truep(C_u_i_fpintegerp_fixed(x)));
2422
2423  if (exponent <= 0) {
2424    C_kontinue(k, C_fix(0));
2425  } else if (exponent == 1) { /* TODO: check significand * 2^exp fits fixnum? */
2426    C_kontinue(k, significand < 0.0 ? C_fix(-1) : C_fix(1));
2427  } else {
2428    C_word kab[C_SIZEOF_CLOSURE(4) + C_SIZEOF_FLONUM], *ka = kab, k2, size,
2429           negp = C_mk_bool(C_flonum_magnitude(x) < 0.0),
2430           sign = C_flonum(&ka, fabs(significand));
2431
2432    k2 = C_closure(&ka, 4, (C_word)flo_to_int_2, k, C_fix(exponent), sign);
2433
2434    size = C_fix(C_BIGNUM_BITS_TO_DIGITS(exponent));
2435    C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
2436  }
2437}
2438
2439static void
2440flo_to_int_2(C_word c, C_word self, C_word result)
2441{
2442  C_word k = C_block_item(self, 1);
2443  C_uword exponent = C_unfix(C_block_item(self, 2)),
2444          *start = C_bignum_digits(result),
2445          *scan = start + C_bignum_size(result);
2446  double significand = C_flonum_magnitude(C_block_item(self, 3));
2447
2448  fabs_frexp_to_digits(exponent, significand, start, scan);
2449  C_kontinue(k, C_bignum_simplify(result));
2450}
2451
2452C_word C_ccall
2453C_u_i_integer_length(C_word x)
2454{
2455  if (x & C_FIXNUM_BIT) {
2456    return C_u_i_fixnum_length(x);
2457  } else {
2458    C_uword result = (C_bignum_size(x) - 1) * C_BIGNUM_DIGIT_LENGTH,
2459            *last_digit = C_bignum_digits(x) + C_bignum_size(x) - 1,
2460            last_digit_length = C_ilen(*last_digit);
2461
2462    /* If *only* the highest bit is set, negating will give one less bit */
2463    if (C_bignum_negativep(x) &&
2464        *last_digit == ((C_uword)1 << (last_digit_length-1))) {
2465      C_uword *startx = C_bignum_digits(x);
2466      while (startx < last_digit && *startx == 0) ++startx;
2467      if (startx == last_digit) result--;
2468    }
2469    return C_fix(result + last_digit_length);
2470  }
2471}
2472
2473void C_ccall /* x is any exact integer but y is _always_ a fixnum */
2474C_u_integer_shift(C_word c, C_word self, C_word k, C_word x, C_word y)
2475{
2476  C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab;
2477
2478  y = C_unfix(y);
2479  if (y == 0 || x == C_fix(0)) { /* Done (no shift) */
2480    C_kontinue(k, x);
2481  } else if (x & C_FIXNUM_BIT) {
2482    if (y < 0) {
2483      /* Don't shift more than a word's length (that's undefined in C!) */
2484      if (-y < C_WORD_SIZE) {
2485        C_kontinue(k, C_fix(C_unfix(x) >> -y));
2486      } else {
2487        C_kontinue(k, (x < 0) ? C_fix(-1) : C_fix(0));
2488      }
2489    } else if (y > 0 && y < C_WORD_SIZE-2 &&
2490               /* After shifting, the length still fits a fixnum */
2491               (C_ilen(C_unfix(x)) + y) < C_WORD_SIZE-2) {
2492      C_kontinue(k, C_fix(C_unfix(x) << y));
2493    } else {
2494      x = C_a_u_i_fix_to_big(&a, x);
2495    }
2496  }
2497
2498  {
2499    C_word ab[C_SIZEOF_CLOSURE(6)], *a = ab,
2500           k2, size, negp, digit_offset, bit_offset;
2501
2502    negp = C_mk_bool(C_bignum_negativep(x));
2503 
2504    if (y > 0) { /* y is guaranteed not to be 0 here */
2505      digit_offset = y / C_BIGNUM_DIGIT_LENGTH;
2506      bit_offset =   y % C_BIGNUM_DIGIT_LENGTH;
2507
2508      k2 = C_closure(&a, 6, (C_word)bignum_actual_shift, k,
2509                     x, C_SCHEME_TRUE, C_fix(digit_offset), C_fix(bit_offset));
2510      size = C_fix(C_bignum_size(x) + digit_offset + 1);
2511      C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
2512    } else if (-y >= C_bignum_size(x) * (C_word)C_BIGNUM_DIGIT_LENGTH) {
2513      /* All bits are shifted out, just return 0 or -1 */
2514      C_kontinue(k, C_truep(negp) ? C_fix(-1) : C_fix(0));
2515    } else {
2516      digit_offset = -y / C_BIGNUM_DIGIT_LENGTH;
2517      bit_offset =   -y % C_BIGNUM_DIGIT_LENGTH;
2518   
2519      k2 = C_closure(&a, 6, (C_word)bignum_actual_shift, k,
2520                     x, C_SCHEME_FALSE, C_fix(digit_offset), C_fix(bit_offset));
2521
2522      size = C_fix(C_bignum_size(x) - digit_offset);
2523      C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
2524    }
2525  }
2526}
2527
2528C_inline C_word maybe_negate_bignum_for_bitwise_op(C_word x, C_word size)
2529{
2530  C_word nx = C_SCHEME_FALSE, xsize;
2531  if (C_bignum_negativep(x)) {
2532    nx = allocate_tmp_bignum(C_fix(size), C_SCHEME_FALSE, C_SCHEME_FALSE);
2533    xsize = C_bignum_size(x);
2534    /* Copy up until requested size, and init any remaining upper digits */
2535    C_memcpy(C_bignum_digits(nx), C_bignum_digits(x),
2536             C_wordstobytes(nmin(size, xsize)));
2537    if (size > xsize)
2538      C_memset(C_bignum_digits(nx)+xsize, 0, C_wordstobytes(size-xsize));
2539    bignum_digits_destructive_negate(nx);
2540  }
2541  return nx;
2542}
2543
2544static void
2545bignum_actual_shift(C_word c, C_word self, C_word result)
2546{
2547  C_word k = C_block_item(self, 1),
2548         x = C_block_item(self, 2),
2549         shift_left = C_truep(C_block_item(self, 3)),
2550         digit_offset = C_unfix(C_block_item(self, 4)),
2551         bit_offset = C_unfix(C_block_item(self, 5));
2552  C_uword *startr = C_bignum_digits(result),
2553          *startx = C_bignum_digits(x),
2554          *endx = startx + C_bignum_size(x),
2555          *endr = startr + C_bignum_size(result);
2556
2557  if (shift_left) {
2558    /* Initialize only the lower digits we're skipping and the MSD */
2559    C_memset(startr, 0, C_wordstobytes(digit_offset));
2560    *(endr-1) = 0;
2561    startr += digit_offset;
2562    /* Can't use bignum_digits_destructive_copy because it assumes
2563     * we want to copy from the start.
2564     */
2565    C_memcpy(startr, startx, C_wordstobytes(endx-startx));
2566    if(bit_offset > 0)
2567      bignum_digits_destructive_shift_left(startr, endr, bit_offset);
2568  } else {
2569    C_word nx, size = C_bignum_size(x) + 1;
2570    if (C_truep(nx = maybe_negate_bignum_for_bitwise_op(x, size))) {
2571      startx = C_bignum_digits(nx); /* update startx; x and endx are unused */
2572    }
2573
2574    startx += digit_offset;
2575    /* Can't use bignum_digits_destructive_copy because that assumes
2576     * target is at least as big as source.
2577     */
2578    C_memcpy(startr, startx, C_wordstobytes(endr-startr));
2579    if(bit_offset > 0)
2580      bignum_digits_destructive_shift_right(startr,endr,bit_offset,C_truep(nx));
2581
2582    if (C_truep(nx)) {
2583      free_tmp_bignum(nx);
2584      bignum_digits_destructive_negate(result);
2585    }
2586  }
2587  C_kontinue(k, C_bignum_simplify(result));
2588}
2589
2590/* This is currently only used by Karatsuba multiplication and
2591 * Burnikel-Ziegler division.  It is not intended as a public API!
2592 */
2593void C_ccall
2594C_u_bignum_extract_digits(C_word c, C_word self, C_word k, C_word x, C_word start, C_word end)
2595{
2596  if (x & C_FIXNUM_BIT) { /* Needed? */
2597    if (C_unfix(start) == 0 && (end == C_SCHEME_FALSE || C_unfix(end) > 0))
2598      C_kontinue(k, x);
2599    else
2600      C_kontinue(k, C_fix(0));
2601  } else {
2602    C_word negp, size;
2603
2604    negp = C_mk_bool(C_bignum_negativep(x)); /* Always false */
2605
2606    start = C_unfix(start);
2607    /* We might get passed larger values than actually fits; pad w/ zeroes */
2608    if (end == C_SCHEME_FALSE) end = C_bignum_size(x);
2609    else end = nmin(C_unfix(end), C_bignum_size(x));
2610    assert(start >= 0);
2611
2612    size = end - start;
2613
2614    if (size == 0 || start >= C_bignum_size(x)) {
2615      C_kontinue(k, C_fix(0));
2616    } else {
2617      C_word k2, kab[C_SIZEOF_CLOSURE(5)], *ka = kab;
2618      k2 = C_closure(&ka, 5, (C_word)bignum_actual_extraction,
2619                     k, x, C_fix(start), C_fix(end));
2620      C_allocate_bignum(5, (C_word)NULL, k2, C_fix(size), negp, C_SCHEME_FALSE);
2621    }
2622  }
2623}
2624
2625static void
2626bignum_actual_extraction(C_word c, C_word self, C_word result)
2627{
2628  C_word k = C_block_item(self, 1),
2629         x = C_block_item(self, 2),
2630         start = C_unfix(C_block_item(self, 3)),
2631         end = C_unfix(C_block_item(self, 4));
2632  C_uword *result_digits = C_bignum_digits(result),
2633          *x_digits = C_bignum_digits(x);
2634
2635  /* Can't use bignum_digits_destructive_copy because that assumes
2636   * target is at least as big as source.
2637   */
2638  C_memcpy(result_digits, x_digits + start, C_wordstobytes(end-start));
2639  C_kontinue(k, C_bignum_simplify(result));
2640}
2641
2642C_regparm C_word C_ccall C_u_i_integer_randomize(C_word seed)
2643{
2644  /* TODO: Rename C_randomize to C_u_i_fixnum_randomize */
2645  if (seed & C_FIXNUM_BIT) {
2646    return C_randomize(seed);
2647  } else {
2648    /*
2649     * This random number generator is very simple. Probably too simple...
2650     */
2651    C_uword *scan = C_bignum_digits(seed),
2652            *end = scan + C_bignum_size(seed), iseed = 0;
2653
2654    /* What a cheap way to initialize the random generator. I feel dirty! */
2655    while (scan < end)
2656      iseed ^= *scan++;
2657
2658    srand(iseed);
2659    return C_SCHEME_UNDEFINED;
2660  }
2661}
2662
2663void C_ccall
2664C_u_integer_random(C_word c, C_word self, C_word k, C_word max)
2665{
2666  /* TODO: for consistency C_random_fixnum should be called C_u_i_fixnum_random */
2667  if (max & C_FIXNUM_BIT) {
2668    C_kontinue(k, C_random_fixnum(max));
2669  } else {
2670    C_word k2, kab[C_SIZEOF_CLOSURE(4)], *ka = kab, size,
2671           max_len, max_bits, max_top_digit, d, negp;
2672
2673    max_len = C_bignum_size(max);
2674    max_top_digit = d = C_bignum_digits(max)[max_len - 1];
2675 
2676    max_bits = (max_len - 1) * C_BIGNUM_DIGIT_LENGTH;
2677    while(d) {
2678      max_bits++;
2679      d >>= 1;
2680    }
2681    /* Subtract/add one because we don't want zero to be over-represented */
2682    size = ((double)rand())/(RAND_MAX + 1.0) * (double)(max_bits - 1);
2683    size = C_fix(C_BIGNUM_BITS_TO_DIGITS(size + 1));
2684
2685    negp = C_mk_bool(C_bignum_negativep(max));
2686    k2 = C_closure(&ka, 4, (C_word)bignum_random_2, k, C_fix(max_top_digit), C_fix(max_len));
2687    C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
2688  }
2689}
2690
2691static void
2692bignum_random_2(C_word c, C_word self, C_word result)
2693{
2694  C_word k = C_block_item(self, 1),
2695         max_top_digit = C_unfix(C_block_item(self, 2)),
2696         max_len = C_unfix(C_block_item(self, 3));
2697  C_uword *scan = C_bignum_digits(result),
2698          *end = scan + C_bignum_size(result); /* Go to just before the end. */
2699
2700  while(scan < end)
2701    *scan++ = ((double)rand())/(RAND_MAX + 1.0) * pow(2.0, C_BIGNUM_DIGIT_LENGTH);
2702  /*
2703   * Last word is special when length is max_len: It must be less than
2704   * max's most significant digit, instead of 2^{digitlen}.
2705   */
2706  if (max_len == C_bignum_size(result))
2707    *scan = ((double)rand())/(RAND_MAX + 1.0) * (double)max_top_digit;
2708  else
2709    *scan = ((double)rand())/(RAND_MAX + 1.0) * pow(2.0, C_BIGNUM_DIGIT_LENGTH);
2710
2711  C_kontinue(k, C_bignum_simplify(result));
2712}
2713
2714C_word C_ccall
2715C_u_i_integer_bit_setp(C_word n, C_word i)
2716{
2717  if (!(i & C_FIXNUM_BIT)) { /* A bit silly, but might be useful */
2718    return C_u_i_integer_negativep(n);
2719  } else if (i & C_INT_SIGN_BIT) {
2720    barf(C_BAD_ARGUMENT_TYPE_NO_UINTEGER_ERROR, "bit-set?", n, i);
2721  } else {
2722    i = C_unfix(i);
2723    if (n & C_FIXNUM_BIT) {
2724      if (i >= C_WORD_SIZE) return C_mk_bool(n & C_INT_SIGN_BIT);
2725      else return C_mk_bool((C_unfix(n) & ((C_word)1 << i)) != 0);
2726    } else {
2727      C_word nn, d;
2728      d = i / C_BIGNUM_DIGIT_LENGTH;
2729      if (d >= C_bignum_size(n)) return C_mk_bool(C_bignum_negativep(n));
2730
2731      if (C_truep(nn = maybe_negate_bignum_for_bitwise_op(n, d))) n = nn;
2732
2733      i %= C_BIGNUM_DIGIT_LENGTH;
2734      d = C_mk_bool((C_bignum_digits(n)[d] & (C_uword)1 << i) != 0);
2735      if (C_truep(nn)) free_tmp_bignum(nn);
2736      return d;
2737    }
2738  }
2739}
2740
2741void C_ccall
2742C_u_2_integer_bitwise_and(C_word c, C_word self, C_word k, C_word x, C_word y)
2743{
2744  if ((x & y) & C_FIXNUM_BIT) {
2745    C_kontinue(k, C_u_fixnum_and(x, y));
2746  } else {
2747    C_word kab[C_SIZEOF_FIX_BIGNUM*2], *ka = kab, negp, size, k2;
2748    if (x & C_FIXNUM_BIT) x = C_a_u_i_fix_to_big(&ka, x);
2749    if (y & C_FIXNUM_BIT) y = C_a_u_i_fix_to_big(&ka, y);
2750
2751    negp = C_mk_bool(C_bignum_negativep(x) && C_bignum_negativep(y));
2752    /* Allow negative 1-bits to propagate */
2753    if (C_bignum_negativep(x) || C_bignum_negativep(y))
2754      size = nmax(C_bignum_size(x), C_bignum_size(y)) + 1;
2755    else
2756      size = nmin(C_bignum_size(x), C_bignum_size(y));
2757
2758    ka = C_alloc(C_SIZEOF_CLOSURE(4)); /* Why faster than static alloc? */
2759    k2 = C_closure(&ka, 4, (C_word)bignum_bitwise_and_2, k, x, y);
2760    C_allocate_bignum(5, (C_word)NULL, k2, C_fix(size), negp, C_SCHEME_FALSE);
2761  }
2762}
2763
2764static void
2765bignum_bitwise_and_2(C_word c, C_word self, C_word result)
2766{
2767  C_word k = C_block_item(self, 1),
2768         x = C_block_item(self, 2),
2769         y = C_block_item(self, 3),
2770         size = C_bignum_size(result), nx, ny;
2771  C_uword *scanr = C_bignum_digits(result),
2772          *endr = scanr + C_bignum_size(result),
2773          *scans1, *ends1, *scans2;
2774
2775  if (C_truep(nx = maybe_negate_bignum_for_bitwise_op(x, size))) x = nx;
2776  if (C_truep(ny = maybe_negate_bignum_for_bitwise_op(y, size))) y = ny;
2777
2778  if (C_bignum_size(x) < C_bignum_size(y)) {
2779    scans1 = C_bignum_digits(x); ends1 = scans1 + C_bignum_size(x);
2780    scans2 = C_bignum_digits(y);
2781  } else {
2782    scans1 = C_bignum_digits(y); ends1 = scans1 + C_bignum_size(y);
2783    scans2 = C_bignum_digits(x);
2784  }
2785
2786  while (scans1 < ends1) *scanr++ = *scans1++ & *scans2++;
2787  C_memset(scanr, 0, C_wordstobytes(endr - scanr));
2788
2789  if (C_truep(nx)) free_tmp_bignum(nx);
2790  if (C_truep(ny)) free_tmp_bignum(ny);
2791  if (C_bignum_negativep(result)) bignum_digits_destructive_negate(result);
2792
2793  C_kontinue(k, C_bignum_simplify(result));
2794}
2795
2796void C_ccall
2797C_u_2_integer_bitwise_ior(C_word c, C_word self, C_word k, C_word x, C_word y)
2798{
2799  if ((x & y) & C_FIXNUM_BIT) {
2800    C_kontinue(k, C_u_fixnum_or(x, y));
2801  } else {
2802    C_word kab[C_SIZEOF_FIX_BIGNUM*2], *ka = kab, negp, size, k2;
2803    if (x & C_FIXNUM_BIT) x = C_a_u_i_fix_to_big(&ka, x);
2804    if (y & C_FIXNUM_BIT) y = C_a_u_i_fix_to_big(&ka, y);
2805
2806    ka = C_alloc(C_SIZEOF_CLOSURE(4)); /* Why faster than static alloc? */
2807    k2 = C_closure(&ka, 4, (C_word)bignum_bitwise_ior_2, k, x, y);
2808    size = C_fix(nmax(C_bignum_size(x), C_bignum_size(y)) + 1);
2809    negp = C_mk_bool(C_bignum_negativep(x) || C_bignum_negativep(y));
2810    C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
2811  }
2812}
2813
2814static void
2815bignum_bitwise_ior_2(C_word c, C_word self, C_word result)
2816{
2817  C_word k = C_block_item(self, 1),
2818         x = C_block_item(self, 2),
2819         y = C_block_item(self, 3),
2820         size = C_bignum_size(result), nx, ny;
2821  C_uword *scanr = C_bignum_digits(result),
2822          *endr = scanr + C_bignum_size(result),
2823          *scans1, *ends1, *scans2, *ends2;
2824
2825  if (C_truep(nx = maybe_negate_bignum_for_bitwise_op(x, size))) x = nx;
2826  if (C_truep(ny = maybe_negate_bignum_for_bitwise_op(y, size))) y = ny;
2827
2828  if (C_bignum_size(x) < C_bignum_size(y)) {
2829    scans1 = C_bignum_digits(x); ends1 = scans1 + C_bignum_size(x);
2830    scans2 = C_bignum_digits(y); ends2 = scans2 + C_bignum_size(y);
2831  } else {
2832    scans1 = C_bignum_digits(y); ends1 = scans1 + C_bignum_size(y);
2833    scans2 = C_bignum_digits(x); ends2 = scans2 + C_bignum_size(x);
2834  }
2835
2836  while (scans1 < ends1) *scanr++ = *scans1++ | *scans2++;
2837  while (scans2 < ends2) *scanr++ = *scans2++;
2838  if (scanr < endr) *scanr++ = 0; /* Only done when result is positive */
2839  assert(scanr == endr);
2840
2841  if (C_truep(nx)) free_tmp_bignum(nx);
2842  if (C_truep(ny)) free_tmp_bignum(ny);
2843  if (C_bignum_negativep(result)) bignum_digits_destructive_negate(result);
2844
2845  C_kontinue(k, C_bignum_simplify(result));
2846}
2847
2848void C_ccall
2849C_u_2_integer_bitwise_xor(C_word c, C_word self, C_word k, C_word x, C_word y)
2850{
2851  if ((x & y) & C_FIXNUM_BIT) {
2852    C_kontinue(k, C_fixnum_xor(x, y));
2853  } else {
2854    C_word kab[C_SIZEOF_FIX_BIGNUM*2], *ka = kab, size, k2, negp;
2855    if (x & C_FIXNUM_BIT) x = C_a_u_i_fix_to_big(&ka, x);
2856    if (y & C_FIXNUM_BIT) y = C_a_u_i_fix_to_big(&ka, y);
2857
2858    ka = C_alloc(C_SIZEOF_CLOSURE(4)); /* Why faster than static alloc? */
2859    k2 = C_closure(&ka, 4, (C_word)bignum_bitwise_xor_2, k, x, y);
2860    size = C_fix(nmax(C_bignum_size(x), C_bignum_size(y)) + 1);
2861    negp = C_mk_bool(C_bignum_negativep(x) != C_bignum_negativep(y));
2862    C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
2863  }
2864}
2865
2866static void
2867bignum_bitwise_xor_2(C_word c, C_word self, C_word result)
2868{
2869  C_word k = C_block_item(self, 1),
2870         x = C_block_item(self, 2),
2871         y = C_block_item(self, 3),
2872         size = C_bignum_size(result), nx, ny;
2873  C_uword *scanr = C_bignum_digits(result),
2874          *endr = scanr + C_bignum_size(result),
2875          *scans1, *ends1, *scans2, *ends2;
2876
2877  if (C_truep(nx = maybe_negate_bignum_for_bitwise_op(x, size))) x = nx;
2878  if (C_truep(ny = maybe_negate_bignum_for_bitwise_op(y, size))) y = ny;
2879
2880  if (C_bignum_size(x) < C_bignum_size(y)) {
2881    scans1 = C_bignum_digits(x); ends1 = scans1 + C_bignum_size(x);
2882    scans2 = C_bignum_digits(y); ends2 = scans2 + C_bignum_size(y);
2883  } else {
2884    scans1 = C_bignum_digits(y); ends1 = scans1 + C_bignum_size(y);
2885    scans2 = C_bignum_digits(x); ends2 = scans2 + C_bignum_size(x);
2886  }
2887
2888  while (scans1 < ends1) *scanr++ = *scans1++ ^ *scans2++;
2889  while (scans2 < ends2) *scanr++ = *scans2++;
2890  if (scanr < endr) *scanr++ = 0; /* Only done when result is positive */
2891  assert(scanr == endr);
2892
2893  if (C_truep(nx)) free_tmp_bignum(nx);
2894  if (C_truep(ny)) free_tmp_bignum(ny);
2895  if (C_bignum_negativep(result)) bignum_digits_destructive_negate(result);
2896
2897  C_kontinue(k, C_bignum_simplify(result));
2898}
2899
2900static void bignum_digits_destructive_negate(C_word result)
2901{
2902  C_uword *scan, *end, digit, sum;
2903
2904  scan = C_bignum_digits(result);
2905  end = scan + C_bignum_size(result);
2906
2907  do {
2908    digit = ~*scan;
2909    sum = digit + 1;
2910    *scan++ = sum;
2911  } while (sum == 0 && scan < end);
2912
2913  for (; scan < end; scan++) {
2914    *scan = ~*scan;
2915  }
2916}
2917
2918/* This is ugly but really cleans up the code below */
2919#define RETURN_Q_AND_OR_R(calc_q, calc_r)                 \
2920  if (C_truep(C_and(return_q, return_r))) {               \
2921    C_values(4, C_SCHEME_UNDEFINED, k, calc_q, calc_r);   \
2922  } else if (C_truep(return_r)) {                         \
2923    C_kontinue(k, calc_r);                                \
2924  } else {                                                \
2925    C_kontinue(k, calc_q);                                \
2926  }
2927
2928/* Lossy; we could be in "quotient&remainder" or "modulo" */
2929#define DIVREM_LOC ((C_truep(C_and(return_q, return_r))) ? "/" :        \
2930                    (C_truep(return_q) ? "quotient" : "remainder"))
2931
2932static C_regparm void
2933basic_divrem(C_word c, C_word self, C_word k, C_word x, C_word y, C_word return_q, C_word return_r)
2934{
2935  if (x & C_FIXNUM_BIT) {
2936    if (y & C_FIXNUM_BIT) {
2937      C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab;
2938      if (y == C_fix(0)) numbers_div_by_zero_error(DIVREM_LOC);
2939
2940      RETURN_Q_AND_OR_R(C_a_u_i_fixnum_quotient_checked(&a, 2, x, y),
2941                        C_u_i_fixnum_remainder_checked(x, y));
2942    } else if (C_immediatep(y)) {
2943      barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, DIVREM_LOC, y);
2944    } else if (C_block_header(y) == C_FLONUM_TAG) {
2945      C_word ab[C_SIZEOF_FLONUM*3], *a = ab;
2946      if (C_flonum_magnitude(y) == 0.0) numbers_div_by_zero_error(DIVREM_LOC);
2947
2948      x = C_a_i_fix_to_flo(&a, 1, x);
2949      RETURN_Q_AND_OR_R(C_a_i_flonum_actual_quotient_checked(&a, 2, x, y),
2950                        C_a_i_flonum_remainder_checked(&a, 2, x, y));
2951    } else if (C_IS_BIGNUM_TYPE(y)) {
2952      integer_divrem(6, (C_word)NULL, k, x, y, return_q, return_r);
2953    } else {
2954      barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, DIVREM_LOC, y);
2955    }
2956  } else if (C_immediatep(x)) {
2957    barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, DIVREM_LOC, x);
2958  } else if (C_block_header(x) == C_FLONUM_TAG) {
2959    if (!C_truep(C_u_i_fpintegerp_fixed(x))) {
2960      barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, DIVREM_LOC, x);
2961    } else if (y & C_FIXNUM_BIT) {
2962      C_word ab[C_SIZEOF_FLONUM*3], *a = ab;
2963      if (y == C_fix(0)) numbers_div_by_zero_error(DIVREM_LOC);
2964
2965      y = C_a_i_fix_to_flo(&a, 1, y);
2966      RETURN_Q_AND_OR_R(C_a_i_flonum_actual_quotient_checked(&a, 2, x, y),
2967                        C_a_i_flonum_remainder_checked(&a, 2, x, y));
2968    } else if (C_immediatep(y)) {
2969      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, DIVREM_LOC, y);
2970    } else if (C_block_header(y) == C_FLONUM_TAG) {
2971      C_word ab[C_SIZEOF_FLONUM*2], *a = ab;
2972      if (C_flonum_magnitude(y) == 0.0) numbers_div_by_zero_error(DIVREM_LOC);
2973
2974      RETURN_Q_AND_OR_R(C_a_i_flonum_actual_quotient_checked(&a, 2, x, y),
2975                        C_a_i_flonum_remainder_checked(&a, 2, x, y));
2976    } else if (C_IS_BIGNUM_TYPE(y)) {
2977      C_word k2, ab[C_SIZEOF_CLOSURE(3)], *a = ab;
2978      x = flo_to_tmp_bignum(x);
2979      k2 = C_closure(&a, 3, (C_word)divrem_intflo_2, k, x);
2980      integer_divrem(6, (C_word)NULL, k2, x, y, return_q, return_r);
2981    } else {
2982      barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, DIVREM_LOC, y);
2983    }
2984  } else if (C_IS_BIGNUM_TYPE(x)) {
2985    if (y & C_FIXNUM_BIT) {
2986      integer_divrem(6, (C_word)NULL, k, x, y, return_q, return_r);
2987    } else if (C_immediatep(y)) {
2988      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, DIVREM_LOC, y);
2989    } else if (C_block_header(y) == C_FLONUM_TAG) {
2990      if (!C_truep(C_u_i_fpintegerp_fixed(y))) {
2991        barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, DIVREM_LOC, y);
2992      } else if (C_flonum_magnitude(y) == 0.0) {
2993        numbers_div_by_zero_error(DIVREM_LOC);
2994      } else {
2995        C_word k2, ab[C_SIZEOF_CLOSURE(3)], *a = ab;
2996        y = flo_to_tmp_bignum(y);
2997        k2 = C_closure(&a, 3, (C_word)divrem_intflo_2, k, y);
2998        integer_divrem(6, (C_word)NULL, k2, x, y, return_q, return_r);
2999      }
3000    } else if (C_IS_BIGNUM_TYPE(y)) {
3001      bignum_divrem(6, (C_word)NULL, k, x, y, return_q, return_r);
3002    } else {
3003      barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, DIVREM_LOC, y);
3004    }
3005  } else {
3006    barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, DIVREM_LOC, x);
3007  }
3008}
3009
3010static void divrem_intflo_2(C_word c, C_word self, ...)
3011{
3012  C_word k = C_block_item(self, 1), x, y;
3013  va_list v;
3014
3015  free_tmp_bignum(C_block_item(self, 2));
3016 
3017  va_start(v, self);
3018  if (c == 2) {
3019    C_word ab[C_SIZEOF_FLONUM], *a = ab;
3020    x = va_arg(v, C_word);
3021    va_end(v);
3022    if (x & C_FIXNUM_BIT) x = C_a_i_fix_to_flo(&a, 1, x);
3023    else x = C_a_u_i_big_to_flo(&a, 1, x);
3024    C_kontinue(k, x);
3025  } else { /* c == 3 */
3026    C_word ab[C_SIZEOF_FLONUM*2], *a = ab;
3027    x = va_arg(v, C_word);
3028    y = va_arg(v, C_word);
3029    va_end(v);
3030
3031    if (x & C_FIXNUM_BIT) x = C_a_i_fix_to_flo(&a, 1, x);
3032    else x = C_a_u_i_big_to_flo(&a, 1, x);
3033    if (y & C_FIXNUM_BIT) y = C_a_i_fix_to_flo(&a, 1, y);
3034    else y = C_a_u_i_big_to_flo(&a, 1, y);
3035    C_values(4, C_SCHEME_UNDEFINED, k, x, y);
3036  }
3037}
3038
3039static void bignum_divrem_fixnum_2(C_word c, C_word self, C_word negated_big)
3040{
3041   C_word k = C_block_item(self, 1),
3042          return_q = C_block_item(self, 2),
3043          return_r = C_block_item(self, 3);
3044   RETURN_Q_AND_OR_R(negated_big, C_fix(0));
3045}
3046
3047static C_regparm void
3048integer_divrem(C_word c, C_word self, C_word k, C_word x, C_word y, C_word return_q, C_word return_r)
3049{
3050  if (!(y & C_FIXNUM_BIT)) { /* y is bignum. */
3051    if (x & C_FIXNUM_BIT) {
3052      /* abs(x) < abs(y), so it will always be [0, x] except for this case: */
3053      if (x == C_fix(C_MOST_NEGATIVE_FIXNUM) &&
3054          C_bignum_negated_fitsinfixnump(y)) {
3055        RETURN_Q_AND_OR_R(C_fix(-1), C_fix(0));
3056      } else {
3057        RETURN_Q_AND_OR_R(C_fix(0), x);
3058      }
3059    } else {
3060      bignum_divrem(6, (C_word)NULL, k, x, y, return_q, return_r);
3061    }
3062  } else if (x & C_FIXNUM_BIT) { /* both x and y are fixnum. */
3063    C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab;
3064    if (y == C_fix(0)) numbers_div_by_zero_error(DIVREM_LOC);
3065
3066    RETURN_Q_AND_OR_R(C_a_u_i_fixnum_quotient_checked(&a, 2, x, y),
3067                      C_u_i_fixnum_remainder_checked(x, y));
3068  } else { /* x is bignum, y is fixnum. */
3069    C_word absy = (y & C_INT_SIGN_BIT) ? -C_unfix(y) : C_unfix(y);
3070
3071    if (y == C_fix(1)) {
3072      RETURN_Q_AND_OR_R(x, C_fix(0));
3073    } else if (y == C_fix(0)) {
3074      numbers_div_by_zero_error(DIVREM_LOC);
3075    } else if (y == C_fix(-1)) {
3076      C_word *ka, k2;
3077      ka = C_alloc(C_SIZEOF_CLOSURE(4));     
3078      k2 = C_closure(&ka, 4, (C_word)bignum_divrem_fixnum_2,
3079                     k, return_q, return_r);
3080      C_u_integer_negate(3, (C_word)NULL, k2, x);
3081    } else if (C_fitsinbignumhalfdigitp(absy) ||
3082               ((((C_uword)1 << (C_ilen(absy)-1)) == absy) &&
3083                C_fitsinfixnump(absy))) {
3084      if (C_truep(return_q)) {
3085        C_word q_negp = C_mk_bool((y & C_INT_SIGN_BIT) ?
3086                                  !(C_bignum_negativep(x)) :
3087                                  C_bignum_negativep(x)),
3088                r_negp = C_mk_bool(C_bignum_negativep(x)),
3089                *ka, k2, size;
3090        ka = C_alloc(C_SIZEOF_CLOSURE(7));
3091        size = C_fix(C_bignum_size(x));
3092        k2 = C_closure(&ka, 7,
3093                       (C_word)bignum_destructive_divide_unsigned_small,
3094                       k, x, C_fix(absy),
3095                       return_q, return_r, r_negp);
3096        C_allocate_bignum(5, (C_word)NULL, k2, size, q_negp, C_SCHEME_FALSE);
3097      } else {
3098        C_word rem;
3099        C_uword next_power = (C_uword)1 << (C_ilen(absy)-1);
3100
3101        if (next_power == absy) { /* Is absy a power of two? */
3102          rem = *(C_bignum_digits(x)) & (next_power - 1);
3103        } else { /* Too bad, we have to do some real work */
3104          rem = bignum_remainder_unsigned_halfdigit(x, absy);
3105        }
3106        C_kontinue(k, C_bignum_negativep(x) ? C_fix(-rem) : C_fix(rem));
3107      }
3108    } else {                    /* Just divide it as two bignums */
3109      C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab;
3110      bignum_divrem(6, (C_word)NULL, k, x, C_a_u_i_fix_to_big(&a, y),
3111                    return_q, return_r);
3112    }
3113  }
3114}
3115
3116static C_regparm void
3117bignum_divrem(C_word c, C_word self, C_word k, C_word x, C_word y, C_word return_q, C_word return_r)
3118{
3119  C_word q_negp = C_mk_bool(C_bignum_negativep(y) ?
3120                            !C_bignum_negativep(x) :
3121                            C_bignum_negativep(x)),
3122         r_negp = C_mk_bool(C_bignum_negativep(x)), size;
3123
3124  switch(bignum_cmp_unsigned(x, y)) {
3125  case 0:
3126    RETURN_Q_AND_OR_R(C_truep(q_negp) ? C_fix(-1) : C_fix(1), C_fix(0));
3127  case -1:
3128    RETURN_Q_AND_OR_R(C_fix(0), x);
3129  case 1:
3130  default:
3131    size = C_bignum_size(x) - C_bignum_size(y);
3132    if (size > C_BURNIKEL_ZIEGLER_DIFF_THRESHOLD) {
3133      try_extended_number("numbers#@bignum-2-divrem-burnikel-ziegler", 5,
3134                          k, x, y, return_q, return_r);
3135    } else if (C_truep(return_q)) {
3136      C_word kab[C_SIZEOF_CLOSURE(6)], *ka = kab, k2;
3137      k2 = C_closure(&ka, 6, (C_word)bignum_divide_2_unsigned, k,
3138                     x, y, return_r, r_negp);
3139      size = C_fix(C_bignum_size(x) + 1 - C_bignum_size(y));
3140      C_allocate_bignum(5, (C_word)NULL, k2, size, q_negp, C_SCHEME_FALSE);
3141    } else { /* We can skip bignum_divide_2_unsigned if we need no quotient */
3142      C_word kab[C_SIZEOF_CLOSURE(7)], *ka = kab, k2;
3143      k2 = C_closure(&ka, 7, (C_word)bignum_divide_2_unsigned_2, k,
3144                     x, y, return_q, return_r, C_SCHEME_UNDEFINED);
3145      size = C_fix(C_bignum_size(x) + 1); /* May need to be normalized */
3146      C_allocate_bignum(5, (C_word)NULL, k2, size, r_negp, C_SCHEME_FALSE);
3147    }
3148  }
3149}
3150
3151static C_word
3152bignum_remainder_unsigned_halfdigit(C_word num, C_word den)
3153{
3154  C_uword *start = C_bignum_digits(num),
3155          *scan = start + C_bignum_size(num),
3156          rem = 0, two_digits;
3157
3158  assert((den > 1) && (C_fitsinbignumhalfdigitp(den)));
3159  while (start < scan) {
3160    two_digits = (*--scan);
3161    rem = C_BIGNUM_DIGIT_COMBINE(rem, C_BIGNUM_DIGIT_HI_HALF(two_digits)) % den;
3162    rem = C_BIGNUM_DIGIT_COMBINE(rem, C_BIGNUM_DIGIT_LO_HALF(two_digits)) % den;
3163  }
3164  return rem;
3165}
3166
3167/* External interface for the above internal divrem functions */
3168void C_ccall
3169C_basic_divrem(C_word c, C_word self, C_word k, C_word x, C_word y)
3170{
3171  if (c != 4) C_bad_argc_2(c, 4, self);
3172  basic_divrem(6, (C_word)NULL, k, x, y, C_SCHEME_TRUE, C_SCHEME_TRUE);
3173}
3174
3175void C_ccall
3176C_u_integer_divrem(C_word c, C_word self, C_word k, C_word x, C_word y)
3177{
3178  integer_divrem(6, (C_word)NULL, k, x, y, C_SCHEME_TRUE, C_SCHEME_TRUE);
3179}
3180
3181void C_ccall
3182C_basic_remainder(C_word c, C_word self, C_word k, C_word x, C_word y)
3183{
3184  if (c != 4) C_bad_argc_2(c, 4, self);
3185  basic_divrem(6, (C_word)NULL, k, x, y, C_SCHEME_FALSE, C_SCHEME_TRUE);
3186}
3187
3188void C_ccall
3189C_u_integer_remainder(C_word c, C_word self, C_word k, C_word x, C_word y)
3190{
3191  integer_divrem(6, (C_word)NULL, k, x, y, C_SCHEME_FALSE, C_SCHEME_TRUE);
3192}
3193
3194void C_ccall
3195C_basic_quotient(C_word c, C_word self, C_word k, C_word x, C_word y)
3196{
3197  if (c != 4) C_bad_argc_2(c, 4, self);
3198  basic_divrem(6, (C_word)NULL, k, x, y, C_SCHEME_TRUE, C_SCHEME_FALSE);
3199}
3200
3201void C_ccall
3202C_u_integer_quotient(C_word c, C_word self, C_word k, C_word x, C_word y)
3203{
3204  integer_divrem(6, (C_word)NULL, k, x, y, C_SCHEME_TRUE, C_SCHEME_FALSE);
3205}
3206
3207/* "small" is either a number that fits a halfdigit, or a power of two */
3208static void
3209bignum_destructive_divide_unsigned_small(C_word c, C_word self, C_word quotient)
3210{
3211  C_word k = C_block_item(self, 1),
3212         numerator = C_block_item(self, 2),
3213         denominator = C_unfix(C_block_item(self, 3)),
3214         /* return_quotient = C_block_item(self, 4), */
3215         return_remainder = C_block_item(self, 5),
3216         remainder_negp = C_block_item(self, 6);
3217  C_uword *start = C_bignum_digits(quotient),
3218          *end = start + C_bignum_size(quotient),
3219          remainder;
3220  int shift_amount;
3221
3222  bignum_digits_destructive_copy(quotient, numerator);
3223
3224  shift_amount = C_ilen(denominator)-1;
3225  if (((C_uword)1 << shift_amount) == denominator) { /* Power of two?  Shift! */
3226    remainder = bignum_digits_destructive_shift_right(start,end,shift_amount,0);
3227    assert(C_ufitsinfixnump(remainder));
3228  } else {
3229    remainder = bignum_digits_destructive_scale_down(start, end, denominator);
3230    assert(C_fitsinbignumhalfdigitp(remainder));
3231  }
3232
3233  quotient = C_bignum_simplify(quotient);
3234
3235  if (C_truep(return_remainder)) {
3236    remainder = C_truep(remainder_negp) ? -remainder : remainder;
3237    C_values(4, C_SCHEME_UNDEFINED, k, quotient, C_fix(remainder));
3238  } else {
3239    C_kontinue(k, quotient);
3240  }
3241}
3242
3243
3244/* Full bignum division */
3245
3246static void
3247bignum_divide_2_unsigned(C_word c, C_word self, C_word quotient)
3248{
3249  C_word k = C_block_item(self, 1),
3250         x = C_block_item(self, 2),
3251         y = C_block_item(self, 3),
3252         size = C_fix(C_bignum_size(x) + 1),
3253         return_r = C_block_item(self, 4),
3254         r_negp = C_block_item(self, 5),
3255         kab[C_SIZEOF_CLOSURE(7)], *ka = kab, k2;
3256
3257  k2 = C_closure(&ka, 7, (C_word)bignum_divide_2_unsigned_2, k,
3258                 x, y, C_SCHEME_TRUE, return_r, quotient);
3259  C_allocate_bignum(5, (C_word)NULL, k2, size, r_negp, C_SCHEME_FALSE);
3260}
3261
3262/* For help understanding this algorithm, see:
3263   Knuth, Donald E., "The Art of Computer Programming",
3264   volume 2, "Seminumerical Algorithms"
3265   section 4.3.1, "Multiple-Precision Arithmetic".
3266
3267   [Yeah, that's a nice book but that particular section is not
3268   helpful at all, which is also pointed out by P. Brinch Hansen's
3269   "Multiple-Length Division Revisited: A Tour Of The Minefield".
3270   That's a more down-to-earth step-by-step explanation of the
3271   algorithm.  Add to this the C implementation in Hacker's Delight
3272   (section 9-2, p141--142) and you may be able to grok this...
3273   ...barely, if you're as math-challenged as I am -- sjamaan]
3274*/
3275
3276static void
3277bignum_divide_2_unsigned_2(C_word c, C_word self, C_word remainder)
3278{
3279  C_word k = C_block_item(self, 1),
3280         numerator = C_block_item(self, 2),
3281         denominator = C_block_item(self, 3),
3282         return_quotient = C_block_item(self, 4),
3283         return_remainder = C_block_item(self, 5),
3284         quotient = C_block_item(self, 6),
3285         length = C_bignum_size(denominator);
3286  C_uword d1 = *(C_bignum_digits(denominator) + length - 1),
3287          *startr = C_bignum_digits(remainder),
3288          *endr = startr + C_bignum_size(remainder);
3289  int shift;
3290
3291  shift = C_BIGNUM_DIGIT_LENGTH - C_ilen(d1); /* nlz */
3292
3293  /* We have to work on halfdigits, so we shift out only the necessary
3294   * amount in order fill out that halfdigit (base is halved).
3295   * This trick is shamelessly stolen from Gauche :)
3296   * See below for part 2 of the trick.
3297   */
3298  if (shift >= C_BIGNUM_HALF_DIGIT_LENGTH)
3299    shift -= C_BIGNUM_HALF_DIGIT_LENGTH;
3300
3301  /* Code below won't always set high halfdigit of quotient, so do it here. */
3302  if (quotient != C_SCHEME_UNDEFINED)
3303    C_bignum_digits(quotient)[C_bignum_size(quotient)-1] = 0;
3304
3305  bignum_digits_destructive_copy(remainder, numerator);
3306  *(endr-1) = 0;    /* Ensure most significant digit is initialised */
3307  if (shift == 0) { /* Already normalized */
3308    bignum_destructive_divide_normalized(remainder, denominator, quotient);
3309  } else { /* Requires normalisation; allocate scratch denominator for this */
3310    C_uword *startnd;
3311    C_word ndenom;
3312
3313    bignum_digits_destructive_shift_left(startr, endr, shift);
3314
3315    ndenom = allocate_tmp_bignum(C_fix(length), C_SCHEME_FALSE, C_SCHEME_FALSE);
3316    startnd = C_bignum_digits(ndenom);
3317    bignum_digits_destructive_copy(ndenom, denominator);
3318    bignum_digits_destructive_shift_left(startnd, startnd+length, shift);
3319
3320    bignum_destructive_divide_normalized(remainder, ndenom, quotient);
3321    if (C_truep(return_remainder)) /* Otherwise, don't bother shifting back */
3322      bignum_digits_destructive_shift_right(startr, endr, shift, 0);
3323
3324    free_tmp_bignum(ndenom);
3325  }
3326
3327  if (C_truep(return_remainder)) {
3328    if (C_truep(return_quotient)) {
3329      C_values(4, C_SCHEME_UNDEFINED, k,
3330               C_bignum_simplify(quotient), C_bignum_simplify(remainder));
3331    } else {
3332      C_kontinue(k, C_bignum_simplify(remainder));
3333    }
3334  } else {
3335    assert(C_truep(return_quotient));
3336    C_kontinue(k, C_bignum_simplify(quotient));
3337  }
3338}
3339
3340static void
3341bignum_destructive_divide_normalized(C_word big_u, C_word big_v, C_word big_q)
3342{
3343  C_uword *v = C_bignum_digits(big_v),
3344          *u = C_bignum_digits(big_u),
3345          *q = big_q == C_SCHEME_UNDEFINED ? NULL : C_bignum_digits(big_q),
3346           p,               /* product of estimated quotient & "denominator" */
3347           hat, qhat, rhat, /* estimated quotient and remainder digit */
3348           vn_1, vn_2;      /* "cached" values v[n-1], v[n-2] */
3349  C_word t, k;              /* Two helpers: temp/final remainder and "borrow" */
3350  /* We use plain ints here, which theoretically may not be enough on
3351   * 64-bit for an insanely huge number, but it is a _lot_ faster.
3352   */
3353  int n = C_bignum_size(big_v) * 2,       /* in halfwords */
3354      m = (C_bignum_size(big_u) * 2) - 2; /* Correct for extra digit */
3355  int i, j;                /* loop  vars */
3356
3357  /* Part 2 of Gauche's aforementioned trick: */
3358  if (C_uhword_ref(v, n-1) == 0) n--;
3359
3360  /* These won't change during the loop, but are used in every step. */
3361  vn_1 = C_uhword_ref(v, n-1);
3362  vn_2 = C_uhword_ref(v, n-2);
3363
3364  /* See also Hacker's Delight, Figure 9-1.  This is almost exactly that. */
3365  for (j = m - n; j >= 0; j--) {
3366    hat = C_BIGNUM_DIGIT_COMBINE(C_uhword_ref(u, j+n), C_uhword_ref(u, j+n-1));
3367    if (hat == 0) {
3368      if (q != NULL) C_uhword_set(q, j, 0);
3369      continue;
3370    }
3371    qhat = hat / vn_1;
3372    rhat = hat % vn_1;
3373
3374    /* Two whiles is faster than one big check with an OR.  Thanks, Gauche! */
3375    while(qhat >= (1L << C_BIGNUM_HALF_DIGIT_LENGTH)) { qhat--; rhat += vn_1; }
3376    while(qhat * vn_2 > C_BIGNUM_DIGIT_COMBINE(rhat, C_uhword_ref(u, j+n-2))
3377          && rhat < (1L << C_BIGNUM_HALF_DIGIT_LENGTH)) {
3378      qhat--;
3379      rhat += vn_1;
3380    }
3381
3382    /* Multiply and subtract */
3383    k = 0;
3384    for (i = 0; i < n; i++) {
3385      p = qhat * C_uhword_ref(v, i);
3386      t = C_uhword_ref(u, i+j) - k - C_BIGNUM_DIGIT_LO_HALF(p);
3387      C_uhword_set(u, i+j, t);
3388      k = C_BIGNUM_DIGIT_HI_HALF(p) - (t >> C_BIGNUM_HALF_DIGIT_LENGTH);
3389    }
3390    t = C_uhword_ref(u,j+n) - k;
3391    C_uhword_set(u, j+n, t);
3392
3393    if (t < 0) {                /* Subtracted too much? */
3394      qhat--;
3395      k = 0;
3396      for (i = 0; i < n; i++) {
3397        t = (C_uword)C_uhword_ref(u, i+j) + C_uhword_ref(v, i) + k;
3398        C_uhword_set(u, i+j, t);
3399        k = t >> C_BIGNUM_HALF_DIGIT_LENGTH;
3400      }
3401      C_uhword_set(u, j+n, (C_uhword_ref(u, j+n) + k));
3402    }
3403    if (q != NULL) C_uhword_set(q, j, qhat);
3404  } /* end j */
3405}
Note: See TracBrowser for help on using the repository browser.