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

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

numbers: First attempt at converting to argvector. Start with wrappers for C_values and allocation continuation closures. Compiles cleanly, but crashes and burns (as expected)

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