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

Last change on this file since 33148 was 33148, checked in by sjamaan, 3 years ago

numbers: Remove yet another pathological case in Burnikel/Ziegler?. We can now revert to the original implementation; reason behind the original slowness was that you shouldn't divide recursively if Y is too small

File size: 115.5 KB
Line 
1/* numbers-c.c
2 *
3 * Copyright 2010-2016 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 CONT_PROCN(divrem_intflo_2, c, 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    av[ 0 ] = C_block_item(err, 0);
203    /* No continuation is passed: '##sys#error-hook' may not return: */
204    av[ 1 ] = C_SCHEME_UNDEFINED;
205    av[ 2 ] = C_fix(code);
206    av[ 3 ] = C_intern2(&a, loc); /* loc is never NULL here, unlike in core */
207
208    va_start(v, loc);
209    for(i = 0; i < c; ++i)
210      av[ i + 4 ] = va_arg(v, C_word);
211    va_end(v);
212
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 CPS_PROC2(C_2_basic_plus, c, self, k, x, y)
510{
511  CPS_BODY2(c, self, k, x, y);
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      CPS_CALL(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      CPS_CALL(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      CPS_CALL(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 CPS_PROC2(C_u_2_integer_plus, c, self, k, x, y)
560{
561  CPS_BODY2(c, self, k, x, y);
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 CPS_PROC2(C_2_basic_minus, c, self, k, x, y)
601{
602  CPS_BODY2(c, self, k, x, y);
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      CPS_CALL(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      CPS_CALL(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      CPS_CALL(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 CPS_PROC2(C_u_2_integer_minus, c, self, k, x, y)
651{
652  CPS_BODY2(c, self, k, x, y);
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 CPS_PROC1(C_basic_abs, c, self, k, x)
877{
878  CPS_BODY1(c, self, k, x);
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    CPS_CALL(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 CPS_PROC1(C_u_integer_abs, c, self, k, x)
897{
898  CPS_BODY1(c, self, k, x);
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    CPS_CALL(C_u_integer_negate, 3, (C_word)NULL, k, x);
904  } else {
905    C_kontinue(k, x);
906  }
907}
908
909void C_ccall CPS_PROC1(C_basic_signum, c, self, k, x)
910{
911  CPS_BODY1(c, self, k, x);
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 CPS_PROC1(C_basic_negate, c, self, k, x)
987{
988  CPS_BODY1(c, self, k, x);
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    CPS_CALL(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 CPS_PROC1(C_u_integer_negate, c, self, k, x)
1005{
1006  CPS_BODY1(c, self, k, x);
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 CPS_PROCN(C_numbers_nequalp, c, self, k)
1044{
1045  CPS_BODYN(c, self, k);
1046  C_word x, y, result;
1047#ifdef ARGVECTOR_CHICKEN
1048  int i = 2;
1049#else
1050  va_list v;
1051#endif
1052
1053  if (c < 4) C_bad_argc_2(c, 4, self);
1054
1055  c -= 2; 
1056#ifdef ARGVECTOR_CHICKEN
1057  x = __av[i++];
1058#else
1059  va_start(v, k);
1060  x = va_arg(v, C_word);
1061#endif
1062
1063  while(--c) {
1064#ifdef ARGVECTOR_CHICKEN
1065    y = __av[i++];
1066#else
1067    y = va_arg(v, C_word);
1068#endif
1069    result = C_i_2_basic_equalp(x, y);
1070    if (result == C_SCHEME_FALSE) break;
1071  }
1072
1073#ifndef ARGVECTOR_CHICKEN
1074  va_end(v);
1075#endif
1076  C_kontinue(k, result);
1077}
1078
1079/* Compare two numbers as ratnums.  Either may be rat-, fix- or bignums */
1080static C_word rat_cmp(C_word x, C_word y)
1081{
1082  C_word ab[C_SIZEOF_FIX_BIGNUM*4], *a = ab, x1, x2, y1, y2,
1083         s, t, ssize, tsize, result, negp;
1084  C_uword *scan;
1085
1086  /* Check for 1 or 0; if x or y is this, the other must be the ratnum */
1087  if (x == C_fix(0)) {        /* Only the sign of y1 matters */
1088    return basic_cmp(x, C_block_item(y, 1), "ratcmp", 0);
1089  } else if (x == C_fix(1)) { /* x1*y1 <> x2*y2 --> y2 <> y1 | x1/x2 = 1/1 */
1090    return basic_cmp(C_block_item(y, 2), C_block_item(y, 1), "ratcmp", 0);
1091  } else if (y == C_fix(0)) { /* Only the sign of x1 matters */
1092    return basic_cmp(C_block_item(x, 1), y, "ratcmp", 0);
1093  } else if (y == C_fix(1)) { /* x1*y1 <> x2*y2 --> x1 <> x2 | y1/y2 = 1/1 */
1094    return basic_cmp(C_block_item(x, 1), C_block_item(x, 2), "ratcmp", 0);
1095  }
1096
1097  /* Extract components x=x1/x2 and y=y1/y2 */
1098  if (x & C_FIXNUM_BIT || C_IS_BIGNUM_TYPE(x)) x1 = x, x2 = C_fix(1);
1099  else x1 = C_block_item(x, 1), x2 = C_block_item(x, 2);
1100
1101  if (y & C_FIXNUM_BIT || C_IS_BIGNUM_TYPE(y)) y1 = y, y2 = C_fix(1);
1102  else y1 = C_block_item(y, 1), y2 = C_block_item(y, 2);
1103
1104  /* We only want to deal with bignums (this is tricky enough) */
1105  if (x1 & C_FIXNUM_BIT) x1 = C_a_u_i_fix_to_big(&a, x1);
1106  if (x2 & C_FIXNUM_BIT) x2 = C_a_u_i_fix_to_big(&a, x2);
1107  if (y1 & C_FIXNUM_BIT) y1 = C_a_u_i_fix_to_big(&a, y1);
1108  if (y2 & C_FIXNUM_BIT) y2 = C_a_u_i_fix_to_big(&a, y2);
1109
1110  /* We multiply using schoolbook method, so this will be very slow in
1111   * extreme cases.  This is a tradeoff we make so that comparisons
1112   * are inlineable, which makes a big difference for the common case.
1113   */
1114  ssize = C_bignum_size(x1) + C_bignum_size(y2);
1115  negp = C_mk_bool(C_bignum_negativep(x1));
1116  s = allocate_tmp_bignum(C_fix(ssize), negp, C_SCHEME_TRUE);
1117  bignum_digits_multiply(x1, y2, s); /* Swap args if x1 < y2? */
1118
1119  tsize = C_bignum_size(y1) + C_bignum_size(x2);
1120  negp = C_mk_bool(C_bignum_negativep(y1));
1121  t = allocate_tmp_bignum(C_fix(tsize), negp, C_SCHEME_TRUE);
1122  bignum_digits_multiply(y1, x2, t); /* Swap args if y1 < x2? */
1123
1124  /* Shorten the numbers if needed */
1125  for (scan = C_bignum_digits(s)+ssize-1; *scan == 0; scan--) ssize--;
1126  C_bignum_mutate_size(s, ssize);
1127  for (scan = C_bignum_digits(t)+tsize-1; *scan == 0; scan--) tsize--;
1128  C_bignum_mutate_size(t, tsize);
1129
1130  result = C_u_i_bignum_cmp(s, t);
1131
1132  free_tmp_bignum(t);
1133  free_tmp_bignum(s);
1134  return result;
1135}
1136
1137/* The primitive comparison operator.  eqp should be 1 if we're only
1138 * interested in equality testing (can speed things up and in case of
1139 * compnums, equality checking is the only available operation).  This
1140 * may return #f, in case there is no answer (for NaNs) or as a quick
1141 * and dirty non-zero answer when eqp is true.  Ugly but effective :)
1142 */
1143static C_word basic_cmp(C_word x, C_word y, char *loc, int eqp)
1144{
1145  if (x & C_FIXNUM_BIT) {
1146    if (y & C_FIXNUM_BIT) {
1147      return C_fix((x < y) ? -1 : ((x > y) ? 1 : 0));
1148    } else if (C_immediatep(y)) {
1149      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y);
1150    } else if (C_block_header(y) == C_FLONUM_TAG) {
1151      return int_flo_cmp(x, y);
1152    } else if (C_IS_BIGNUM_TYPE(y)) {
1153      C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab;
1154      return C_u_i_bignum_cmp(C_a_u_i_fix_to_big(&a, x), y);
1155    } else if (C_header_bits(y) == C_STRUCTURE_TYPE &&
1156               C_block_item(y, 0) == ratnum_type_tag) {
1157      if (eqp) return C_SCHEME_FALSE;
1158      else return rat_cmp(x, y);
1159    } else if (C_block_item(y, 0) == compnum_type_tag) {
1160      if (eqp) return C_SCHEME_FALSE;
1161      else barf(C_BAD_ARGUMENT_TYPE_COMPLEX_NO_ORDERING_ERROR, loc, y);
1162    } else {
1163      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y);
1164    }
1165  } else if (C_immediatep(x)) {
1166    barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, x);
1167  } else if (C_block_header(x) == C_FLONUM_TAG) {
1168    if (y & C_FIXNUM_BIT) {
1169      return flo_int_cmp(x, y);
1170    } else if (C_immediatep(y)) {
1171      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y);
1172    } else if (C_block_header(y) == C_FLONUM_TAG) {
1173      double a = C_flonum_magnitude(x), b = C_flonum_magnitude(y);
1174      if (C_isnan(a) || C_isnan(b)) return C_SCHEME_FALSE; /* "mu" */
1175      else return C_fix((a < b) ? -1 : ((a > b) ? 1 : 0));
1176    } else if (C_IS_BIGNUM_TYPE(y)) {
1177      return flo_int_cmp(x, y);
1178    } else if (C_header_bits(y) == C_STRUCTURE_TYPE) {
1179      if (C_block_item(y, 0) == ratnum_type_tag) {
1180        return flo_rat_cmp(x, y);
1181      } else if (C_block_item(y, 0) == compnum_type_tag) {
1182        if (eqp) return C_SCHEME_FALSE;
1183        else barf(C_BAD_ARGUMENT_TYPE_COMPLEX_NO_ORDERING_ERROR, loc, y);
1184      } else {
1185        barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y);
1186      }
1187    } else {
1188      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y);
1189    }
1190  } else if (C_IS_BIGNUM_TYPE(x)) {
1191    if (y & C_FIXNUM_BIT) {
1192      C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab;
1193      return C_u_i_bignum_cmp(x, C_a_u_i_fix_to_big(&a, y));
1194    } else if (C_immediatep(y)) {
1195      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y);
1196    } else if (C_block_header(y) == C_FLONUM_TAG) {
1197      return int_flo_cmp(x, y);
1198    } else if (C_IS_BIGNUM_TYPE(y)) {
1199      return C_u_i_bignum_cmp(x, y);
1200    } else if (C_header_bits(y) == C_STRUCTURE_TYPE &&
1201               C_block_item(y, 0) == ratnum_type_tag) {
1202      if (eqp) return C_SCHEME_FALSE;
1203      else return rat_cmp(x, y);
1204    } else if (C_block_item(y, 0) == compnum_type_tag) {
1205      if (eqp) return C_SCHEME_FALSE;
1206      else barf(C_BAD_ARGUMENT_TYPE_COMPLEX_NO_ORDERING_ERROR, loc, y);
1207    } else {
1208      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y);
1209    }
1210  } else if (C_header_bits(x) == C_STRUCTURE_TYPE &&
1211             (C_block_item(x, 0) == ratnum_type_tag)) {
1212    if (y & C_FIXNUM_BIT) {
1213      if (eqp) return C_SCHEME_FALSE;
1214      else return rat_cmp(x, y);
1215    } else if (C_immediatep(y)) {
1216      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y);
1217    } else if (C_block_header(y) == C_FLONUM_TAG) {
1218      return rat_flo_cmp(x, y);
1219    } else if (C_IS_BIGNUM_TYPE(y)) {
1220      if (eqp) return C_SCHEME_FALSE;
1221      else return rat_cmp(x, y);
1222    } else if (C_header_bits(y) == C_STRUCTURE_TYPE &&
1223               (C_block_item(y, 0) == ratnum_type_tag)) {
1224      if (eqp) {
1225        return C_and(C_and(C_u_i_2_integer_equalp(C_block_item(x, 1),
1226                                                  C_block_item(y, 1)),
1227                           C_u_i_2_integer_equalp(C_block_item(x, 2),
1228                                                  C_block_item(y, 2))),
1229                     C_fix(0));
1230      } else {
1231        return rat_cmp(x, y);
1232      }
1233    } else if (C_header_bits(y) == C_STRUCTURE_TYPE &&
1234               (C_block_item(y, 0) == compnum_type_tag)) {
1235      if (eqp) return C_SCHEME_FALSE;
1236      else barf(C_BAD_ARGUMENT_TYPE_COMPLEX_NO_ORDERING_ERROR, loc, y);
1237    } else {
1238      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y);
1239    }
1240  } else if (C_header_bits(x) == C_STRUCTURE_TYPE &&
1241             (C_block_item(x, 0) == compnum_type_tag)) {
1242    if (!eqp) {
1243      barf(C_BAD_ARGUMENT_TYPE_COMPLEX_NO_ORDERING_ERROR, loc, x);
1244    } else if (y & C_FIXNUM_BIT) {
1245      return C_SCHEME_FALSE;
1246    } else if (C_immediatep(y)) {
1247      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y);
1248    } else if (C_block_header(y) == C_FLONUM_TAG ||
1249               C_IS_BIGNUM_TYPE(y) ||
1250               (C_header_bits(y) == C_STRUCTURE_TYPE &&
1251                (C_block_item(y, 0) == ratnum_type_tag))) {
1252      return C_SCHEME_FALSE;
1253    } else if (C_header_bits(y) == C_STRUCTURE_TYPE &&
1254               (C_block_item(y, 0) == compnum_type_tag)) {
1255      return C_and(C_and(C_i_2_basic_equalp(C_block_item(x, 1),
1256                                            C_block_item(y, 1)),
1257                         C_i_2_basic_equalp(C_block_item(x, 2),
1258                                            C_block_item(y, 2))),
1259                   C_fix(0));
1260    } else {
1261      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, y);
1262    }
1263  } else {
1264    barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, loc, x);
1265  }
1266}
1267
1268C_regparm C_word C_fcall
1269C_i_2_basic_equalp(C_word x, C_word y)
1270{
1271   return C_mk_bool(basic_cmp(x, y, "=", 1) == C_fix(0));
1272}
1273
1274C_word C_ccall C_u_i_2_integer_equalp(C_word x, C_word y)
1275{
1276  if (x & C_FIXNUM_BIT)
1277    return C_mk_bool(x == y);
1278  else if (y & C_FIXNUM_BIT)
1279    return C_SCHEME_FALSE;
1280  else
1281    return C_mk_bool(C_u_i_bignum_cmp(x, y) == C_fix(0));
1282}
1283
1284/* TODO: Rename to C_lessp */
1285void C_ccall CPS_PROCN(C_numbers_lessp, c, self, k)
1286{
1287  CPS_BODYN(c, self, k);
1288  C_word x, y, result;
1289#ifdef ARGVECTOR_CHICKEN
1290  int i = 2;
1291#else
1292  va_list v;
1293#endif
1294 
1295  if (c < 4) C_bad_argc_2(c, 4, self);
1296
1297  c -= 2; 
1298#ifdef ARGVECTOR_CHICKEN
1299  x = __av[i++];
1300#else
1301  va_start(v, k);
1302  x = va_arg(v, C_word);
1303#endif
1304
1305  while(--c) {
1306#ifdef ARGVECTOR_CHICKEN
1307    y = __av[i++];
1308#else
1309    y = va_arg(v, C_word);
1310#endif
1311    result = C_i_2_basic_lessp(x, y);
1312    if (result == C_SCHEME_FALSE) break;
1313    x = y;
1314  }
1315
1316#ifndef ARGVECTOR_CHICKEN
1317  va_end(v);
1318#endif
1319  C_kontinue(k, result);
1320}
1321
1322/* TODO: Rename to C_less_or_equal_p */
1323void C_ccall CPS_PROCN(C_numbers_less_or_equal_p, c, self, k)
1324{
1325  CPS_BODYN(c, self, k);
1326  C_word x, y, result;
1327#ifdef ARGVECTOR_CHICKEN
1328  int i = 2;
1329#else
1330  va_list v;
1331#endif
1332
1333  if (c < 4) C_bad_argc_2(c, 4, self);
1334
1335  c -= 2; 
1336#ifdef ARGVECTOR_CHICKEN
1337  x = __av[i++];
1338#else
1339  va_start(v, k);
1340  x = va_arg(v, C_word);
1341#endif
1342
1343  while(--c) {
1344#ifdef ARGVECTOR_CHICKEN
1345    y = __av[i++];
1346#else
1347    y = va_arg(v, C_word);
1348#endif
1349    result = C_i_2_basic_less_or_equalp(x, y);
1350    if (result == C_SCHEME_FALSE) break;
1351    x = y;
1352  }
1353
1354#ifndef ARGVECTOR_CHICKEN
1355  va_end(v);
1356#endif
1357  C_kontinue(k, result);
1358}
1359
1360C_regparm C_word C_fcall C_i_2_basic_lessp(C_word x, C_word y)
1361{
1362   return C_mk_bool(basic_cmp(x, y, "<", 0) == C_fix(-1));
1363}
1364
1365C_regparm C_word C_fcall C_i_2_basic_less_or_equalp(C_word x, C_word y)
1366{
1367   C_word res = basic_cmp(x, y, "<=", 0);
1368   return C_mk_bool(res == C_fix(0) || res == C_fix(-1));
1369}
1370
1371C_word C_ccall
1372C_u_i_2_integer_lessp(C_word x, C_word y)
1373{
1374  if (x & C_FIXNUM_BIT) {
1375    if (y & C_FIXNUM_BIT) {
1376      return C_mk_bool(C_unfix(x) < C_unfix(y));
1377    } else {
1378      return C_mk_nbool(C_bignum_negativep(y));
1379    }
1380  } else if (y & C_FIXNUM_BIT) {
1381    return C_mk_bool(C_bignum_negativep(x));
1382  } else {
1383    return C_mk_bool(C_u_i_bignum_cmp(x, y) == C_fix(-1));
1384  }
1385}
1386
1387C_word C_ccall
1388C_u_i_2_integer_less_or_equalp(C_word x, C_word y)
1389{
1390  if (x & C_FIXNUM_BIT) {
1391    if (y & C_FIXNUM_BIT) {
1392      return C_mk_bool(C_unfix(x) <= C_unfix(y));
1393    } else {
1394      return C_mk_nbool(C_bignum_negativep(y));
1395    }
1396  } else if (y & C_FIXNUM_BIT) {
1397    return C_mk_bool(C_bignum_negativep(x));
1398  } else {
1399    C_word res = C_u_i_bignum_cmp(x, y);
1400    return C_mk_bool(res == C_fix(0) || res == C_fix(-1));
1401  }
1402}
1403
1404/* TODO: Rename to C_greater_p */
1405void C_ccall CPS_PROCN(C_numbers_greaterp, c, self, k)
1406{
1407  CPS_BODYN(c, self, k);
1408  C_word x, y, result;
1409#ifdef ARGVECTOR_CHICKEN
1410  int i = 2;
1411#else
1412  va_list v;
1413#endif
1414
1415  if (c < 4) C_bad_argc_2(c, 4, self);
1416
1417  c -= 2; 
1418#ifdef ARGVECTOR_CHICKEN
1419  x = __av[i++];
1420#else
1421  va_start(v, k);
1422  x = va_arg(v, C_word);
1423#endif
1424
1425  while(--c) {
1426#ifdef ARGVECTOR_CHICKEN
1427    y = __av[i++];
1428#else
1429    y = va_arg(v, C_word);
1430#endif
1431    result = C_i_2_basic_greaterp(x, y);
1432    if (result == C_SCHEME_FALSE) break;
1433    x = y;
1434  }
1435
1436#ifndef ARGVECTOR_CHICKEN
1437  va_end(v);
1438#endif
1439  C_kontinue(k, result);
1440}
1441
1442/* TODO: Rename to C_greater_or_equal_p */
1443void C_ccall CPS_PROCN(C_numbers_greater_or_equal_p, c, self, k)
1444{
1445  CPS_BODYN(c, self, k);
1446  C_word x, y, result;
1447#ifdef ARGVECTOR_CHICKEN
1448  int i = 2;
1449#else
1450  va_list v;
1451#endif
1452
1453  if (c < 4) C_bad_argc_2(c, 4, self);
1454
1455  c -= 2; 
1456#ifdef ARGVECTOR_CHICKEN
1457  x = __av[i++];
1458#else
1459  va_start(v, k);
1460  x = va_arg(v, C_word);
1461#endif
1462
1463  while(--c) {
1464#ifdef ARGVECTOR_CHICKEN
1465    y = __av[i++];
1466#else
1467    y = va_arg(v, C_word);
1468#endif
1469    result = C_i_2_basic_greater_or_equalp(x, y);
1470    if (result == C_SCHEME_FALSE) break;
1471    x = y;
1472  }
1473
1474#ifndef ARGVECTOR_CHICKEN
1475  va_end(v);
1476#endif
1477  C_kontinue(k, result);
1478}
1479
1480C_regparm C_word C_fcall C_i_2_basic_greaterp(C_word x, C_word y)
1481{
1482   return C_mk_bool(basic_cmp(x, y, ">", 0) == C_fix(1));
1483}
1484
1485C_regparm C_word C_fcall C_i_2_basic_greater_or_equalp(C_word x, C_word y)
1486{
1487   C_word res = basic_cmp(x, y, ">=", 0);
1488   return C_mk_bool(res == C_fix(0) || res == C_fix(1));
1489}
1490
1491C_word C_ccall
1492C_u_i_2_integer_greaterp(C_word x, C_word y)
1493{
1494  if (x & C_FIXNUM_BIT) {
1495    if (y & C_FIXNUM_BIT) {
1496      return C_mk_bool(C_unfix(x) > C_unfix(y));
1497    } else {
1498      return C_mk_bool(C_bignum_negativep(y));
1499    }
1500  } else if (y & C_FIXNUM_BIT) {
1501    return C_mk_nbool(C_bignum_negativep(x));
1502  } else {
1503    return C_mk_bool(C_u_i_bignum_cmp(x, y) == C_fix(1));
1504  }
1505}
1506
1507C_word C_ccall
1508C_u_i_2_integer_greater_or_equalp(C_word x, C_word y)
1509{
1510  if (x & C_FIXNUM_BIT) {
1511    if (y & C_FIXNUM_BIT) {
1512      return C_mk_bool(C_unfix(x) >= C_unfix(y));
1513    } else {
1514      return C_mk_bool(C_bignum_negativep(y));
1515    }
1516  } else if (y & C_FIXNUM_BIT) {
1517    return C_mk_nbool(C_bignum_negativep(x));
1518  } else {
1519    C_word res = C_u_i_bignum_cmp(x, y);
1520    return C_mk_bool(res == C_fix(0) || res == C_fix(1));
1521  }
1522}
1523
1524/* This is a bit weird: We have to compare flonums as bignums due to
1525 * precision loss on 64-bit platforms.  For simplicity, we convert
1526 * fixnums to bignums here.
1527 */
1528static C_word int_flo_cmp(C_word intnum, C_word flonum)
1529{
1530  C_word ab[C_SIZEOF_FIX_BIGNUM + C_SIZEOF_FLONUM], *a = ab, x, y, res;
1531  double i, f;
1532
1533  f = C_flonum_magnitude(flonum);
1534
1535  if (C_isnan(f)) {
1536    return C_SCHEME_FALSE; /* "mu" */
1537  } else if (C_isinf(f)) {
1538    return C_fix((f > 0.0) ? -1 : 1); /* x is smaller if f is +inf.0 */
1539  } else {
1540    f = modf(f, &i);
1541
1542    x = (intnum & C_FIXNUM_BIT) ? C_a_u_i_fix_to_big(&a, intnum) : intnum;
1543    y = flo_to_tmp_bignum(C_flonum(&a, i));
1544
1545    res = C_u_i_bignum_cmp(x, y);
1546    free_tmp_bignum(y);
1547
1548    if (res == C_fix(0)) /* Use fraction to break tie. If f > 0, x is smaller */
1549      return C_fix((f > 0.0) ? -1 : ((f < 0.0) ? 1 : 0));
1550    else
1551      return res;
1552  }
1553}
1554
1555/* For convenience (ie, to reduce the degree of mindfuck) */
1556static C_word flo_int_cmp(C_word flonum, C_word intnum)
1557{
1558  C_word res = int_flo_cmp(intnum, flonum);
1559  switch(res) {
1560  case C_fix(1): return C_fix(-1);
1561  case C_fix(-1): return C_fix(1);
1562  default: return res; /* Can be either C_fix(0) or C_SCHEME_FALSE(!) */
1563  }
1564}
1565
1566/* This code is completely braindead, but at least it allows us to do
1567 * inline comparisons!
1568 */
1569static C_word rat_flo_cmp(C_word ratnum, C_word flonum)
1570{
1571  C_word ab[C_SIZEOF_FIX_BIGNUM * 2 + C_SIZEOF_FLONUM], *a = ab,
1572         num, denom, ibig, res, nscaled, iscaled, negp;
1573  C_uword *scan;
1574  int shift_amount, ilen, nlen;
1575  double i, f;
1576
1577  f = C_flonum_magnitude(flonum);
1578
1579  if (C_isnan(f)) {
1580    return C_SCHEME_FALSE; /* "mu" */
1581  } else if (C_isinf(f)) {
1582    return C_fix((f > 0.0) ? -1 : 1); /* x is smaller if f is +inf.0 */
1583  } else {
1584    /* Scale up the floating-point number to become a whole integer,
1585     * and remember power of two (# of bits) to shift the numerator.
1586     */
1587    shift_amount = 0;
1588
1589    /* TODO: This doesn't work for denormalized flonums! */
1590    while (modf(f, &i) != 0.0) {
1591      f = ldexp(f, 1);
1592      shift_amount++;
1593    }
1594
1595    i = f; /* TODO: split i and f so it'll work for denormalized flonums */
1596
1597    num = C_block_item(ratnum, 1);
1598    num = (num & C_FIXNUM_BIT) ? C_a_u_i_fix_to_big(&a, num) : num;
1599
1600    if (C_bignum_negativep(num) && i >= 0.0) { /* Save time if signs differ */
1601      return C_fix(-1);
1602    } else if (!C_bignum_negativep(num) && i <= 0.0) { /* num is never 0 */
1603      return C_fix(1);
1604    } else {
1605      negp = C_mk_bool(C_bignum_negativep(num));
1606
1607      denom = C_block_item(ratnum, 2);
1608      denom = (denom & C_FIXNUM_BIT) ? C_a_u_i_fix_to_big(&a, denom) : denom;
1609
1610      ibig = flo_to_tmp_bignum(C_flonum(&a, i));
1611
1612      nlen = C_bignum_size(num) + C_bignum_size(denom);
1613      ilen = C_bignum_size(ibig) + C_bignum_size(denom);
1614
1615      /* Now, multiply the scaled flonum by the denominator, so we can
1616       * compare it directly to the scaled numerator.  Unfortunately,
1617       * this won't use Karatsuba multiplication, so for large numbers
1618       * it will be slower than it could be if comparisons were done
1619       * in CPS context.
1620       */
1621      iscaled = allocate_tmp_bignum(C_fix(ilen), negp, C_SCHEME_TRUE);
1622      bignum_digits_multiply(denom, ibig, iscaled); /* Swap args if i < d? */
1623      free_tmp_bignum(ibig);
1624
1625      nlen += C_BIGNUM_BITS_TO_DIGITS(shift_amount);
1626      nscaled = allocate_tmp_bignum(C_fix(nlen), negp, C_SCHEME_TRUE);
1627
1628      scan = C_bignum_digits(nscaled) + shift_amount / C_BIGNUM_DIGIT_LENGTH;
1629      C_memcpy(scan, C_bignum_digits(num), C_wordstobytes(C_bignum_size(num)));
1630      shift_amount = shift_amount % C_BIGNUM_DIGIT_LENGTH;
1631      if(shift_amount > 0) {
1632        bignum_digits_destructive_shift_left(
1633         scan, C_bignum_digits(nscaled) + nlen, shift_amount);
1634      }
1635
1636      /* Shorten the numbers if needed */
1637      for (scan = C_bignum_digits(iscaled)+ilen-1; *scan == 0; scan--) ilen--;
1638      C_bignum_mutate_size(iscaled, ilen);
1639      for (scan = C_bignum_digits(nscaled)+nlen-1; *scan == 0; scan--) nlen--;
1640      C_bignum_mutate_size(nscaled, nlen);
1641
1642      /* Finally, we're ready to compare them! */
1643      res = C_u_i_bignum_cmp(nscaled, iscaled);
1644      free_tmp_bignum(nscaled);
1645      free_tmp_bignum(iscaled);
1646
1647      return res;
1648    }
1649  }
1650}
1651
1652static C_word flo_rat_cmp(C_word flonum, C_word ratnum)
1653{
1654  C_word res = rat_flo_cmp(ratnum, flonum);
1655  switch(res) {
1656  case C_fix(1): return C_fix(-1);
1657  case C_fix(-1): return C_fix(1);
1658  default: return res; /* Can be either C_fix(0) or C_SCHEME_FALSE(!) */
1659  }
1660}
1661
1662C_word
1663C_u_i_bignum_cmp(C_word x, C_word y)
1664{
1665  if (C_bignum_negativep(x)) {
1666    if (C_bignum_negativep(y)) { /* Largest negative number is smallest */
1667      return C_fix(bignum_cmp_unsigned(y, x));
1668    } else {
1669      return C_fix(-1);
1670    }
1671  } else {
1672    if (C_bignum_negativep(y)) {
1673      return C_fix(1);
1674    } else {
1675      return C_fix(bignum_cmp_unsigned(x, y));
1676    }
1677  }
1678}
1679
1680/* NOTE: If C_allocate_bignum is to be callable from Scheme, it needs
1681 * to be converted to argvector.  But that's a very big change.
1682 */
1683void C_ccall
1684C_allocate_bignum(C_word c, C_word self, C_word k, C_word size, C_word negp, C_word initp)
1685{
1686#ifdef ARGVECTOR_CHICKEN
1687  C_word kab[C_SIZEOF_CLOSURE(3)], *ka = kab, av[6];
1688
1689  av[ 0 ] = (C_word)NULL;   /* No "self" closure */
1690  av[ 1 ] = C_closure(&ka, 3, (C_word)allocate_bignum_2, k, negp);
1691  av[ 2 ] = C_bytes(C_fixnum_plus(size, C_fix(1))); /* Add header */
1692  av[ 3 ] = C_SCHEME_TRUE;  /* Byte vector */
1693  av[ 4 ] = C_and(initp, C_make_character('\0'));
1694  av[ 5 ] = C_SCHEME_FALSE; /* Don't align at 8 bytes */
1695  C_allocate_vector(6, av);
1696#else
1697  C_word kab[C_SIZEOF_CLOSURE(3)], *ka = kab, k2, init;
1698  k2 = C_closure(&ka, 3, (C_word)allocate_bignum_2, k, negp);
1699
1700  init = C_and(initp, C_make_character('\0'));
1701  C_allocate_vector(6, (C_word)NULL, k2,
1702                    C_bytes(C_fixnum_plus(size, C_fix(1))), /* Add header */
1703                    /* Byte vec, initialization, align at 8 bytes */
1704                    C_SCHEME_TRUE, init, C_SCHEME_FALSE);
1705#endif
1706}
1707
1708static void CONT_PROC(allocate_bignum_2, c, self, bigvec)
1709{
1710  CONT_BODY(self, bigvec);
1711  C_word ab[C_SIZEOF_STRUCTURE(2)], *a = ab, bignum,
1712         k = C_block_item(self, 1),
1713         negp = C_truep(C_block_item(self, 2)) ? 1 : 0,
1714         tagvec = CHICKEN_gc_root_ref(tags);
1715
1716  C_set_block_item(bigvec, 0, negp);
1717
1718  bignum = C_a_i_record2(&a, 2, C_block_item(tagvec, BIG_TAG), bigvec);
1719  C_kontinue(k, bignum);
1720}
1721
1722static C_word
1723allocate_tmp_bignum(C_word size, C_word negp, C_word initp)
1724{
1725  C_word *mem = malloc(C_wordstobytes(C_SIZEOF_BIGNUM(C_unfix(size)))),
1726          bigvec = (C_word)(mem + C_SIZEOF_STRUCTURE(2)),
1727          tagvec = CHICKEN_gc_root_ref(tags);
1728  if (mem == NULL) abort();     /* TODO: panic */
1729 
1730  C_block_header_init(bigvec, (C_STRING_TYPE | C_wordstobytes(C_unfix(size)+1)));
1731  C_set_block_item(bigvec, 0, C_truep(negp));
1732
1733  if (C_truep(initp)) {
1734    C_memset(((C_uword *)C_data_pointer(bigvec))+1,
1735             0, C_wordstobytes(C_unfix(size)));
1736  }
1737
1738  return C_a_i_record2(&mem, 2, C_block_item(tagvec, BIG_TAG), bigvec);
1739}
1740
1741/* Simplification: scan trailing zeroes, then return a fixnum if the
1742 * value fits, or trim the bignum's length. */
1743C_word C_ccall
1744C_bignum_simplify(C_word big)
1745{
1746  C_uword *start = C_bignum_digits(big),
1747          *last_digit = start + C_bignum_size(big) - 1,
1748          *scan = last_digit, tmp;
1749  int length;
1750
1751  while (scan >= start && *scan == 0)
1752    scan--;
1753  length = scan - start + 1;
1754 
1755  switch(length) {
1756  case 0:
1757    return C_fix(0);
1758  case 1:
1759    tmp = *start;
1760    if (C_bignum_negativep(big) ?
1761        !(tmp & C_INT_SIGN_BIT) && C_fitsinfixnump(-(C_word)tmp) :
1762        C_ufitsinfixnump(tmp))
1763      return C_bignum_negativep(big) ? C_fix(-(C_word)tmp) : C_fix(tmp);
1764    /* FALLTHROUGH */
1765  default:
1766    if (scan < last_digit) C_bignum_mutate_size(big, length);
1767    return big;
1768  }
1769}
1770
1771static C_uword
1772bignum_digits_destructive_scale_up_with_carry(C_uword *start, C_uword *end, C_uword factor, C_uword carry)
1773{
1774  C_uword digit, p;
1775
1776  assert(C_fitsinbignumhalfdigitp(carry));
1777  assert(C_fitsinbignumhalfdigitp(factor));
1778
1779  /* See fixnum_times.  Substitute xlo = factor, xhi = 0, y = digit
1780   * and simplify the result to reduce variable usage.
1781   */
1782  while (start < end) {
1783    digit = (*start);
1784
1785    p = factor * C_BIGNUM_DIGIT_LO_HALF(digit) + carry;
1786    carry = C_BIGNUM_DIGIT_LO_HALF(p);
1787
1788    p = factor * C_BIGNUM_DIGIT_HI_HALF(digit) + C_BIGNUM_DIGIT_HI_HALF(p);
1789    (*start++) = C_BIGNUM_DIGIT_COMBINE(C_BIGNUM_DIGIT_LO_HALF(p), carry);
1790    carry = C_BIGNUM_DIGIT_HI_HALF(p);
1791  }
1792  return carry;
1793}
1794
1795static C_uword
1796bignum_digits_destructive_scale_down(C_uword *start, C_uword *end, C_uword denominator)
1797{
1798  C_uword digit, k = 0;
1799  C_uhword q_j_hi, q_j_lo;
1800
1801  /* Single digit divisor case from Hacker's Delight, Figure 9-1,
1802   * adapted to modify u[] in-place instead of writing to q[].
1803   */
1804  while (start < end) {
1805    digit = (*--end);
1806
1807    k = C_BIGNUM_DIGIT_COMBINE(k, C_BIGNUM_DIGIT_HI_HALF(digit)); /* j */
1808    q_j_hi = k / denominator;
1809    k -= q_j_hi * denominator;
1810
1811    k = C_BIGNUM_DIGIT_COMBINE(k, C_BIGNUM_DIGIT_LO_HALF(digit)); /* j-1 */
1812    q_j_lo = k / denominator;
1813    k -= q_j_lo * denominator;
1814   
1815    *end = C_BIGNUM_DIGIT_COMBINE(q_j_hi, q_j_lo);
1816  }
1817  return k;
1818}
1819
1820static C_uword
1821bignum_digits_destructive_shift_right(C_uword *start, C_uword *end, int shift_right, int negp)
1822{
1823  int shift_left = C_BIGNUM_DIGIT_LENGTH - shift_right;
1824  C_uword digit, carry = negp ? ((~(C_uword)0) << shift_left) : 0;
1825
1826  assert(shift_right < C_BIGNUM_DIGIT_LENGTH);
1827
1828  while (start < end) {
1829    digit = *(--end);
1830    *end = (digit >> shift_right) | carry;
1831    carry = digit << shift_left;
1832  }
1833  return carry >> shift_left; /* The bits that were shifted out to the right */
1834}
1835
1836static C_uword
1837bignum_digits_destructive_shift_left(C_uword *start, C_uword *end, int shift_left)
1838{
1839  C_uword carry = 0, digit;
1840  int shift_right = C_BIGNUM_DIGIT_LENGTH - shift_left;
1841
1842  assert(shift_left < C_BIGNUM_DIGIT_LENGTH);
1843
1844  while (start < end) {
1845    digit = *start;
1846    (*start++) = (digit << shift_left) | carry;
1847    carry = digit >> shift_right;
1848  }
1849  return carry;  /* This would end up as most significant digit if it fit */
1850}
1851
1852void C_ccall CPS_PROC2(C_2_basic_times, c, self, k, x, y)
1853{
1854  CPS_BODY2(c, self, k, x, y);
1855  if (x & C_FIXNUM_BIT) {
1856    if (y & C_FIXNUM_BIT) {
1857      C_word *a = C_alloc(C_SIZEOF_BIGNUM(2));
1858      C_kontinue(k, C_a_u_i_2_fixnum_times(&a, 2, x, y));
1859    } else if (C_immediatep(y)) {
1860      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", y);
1861    } else if (C_block_header(y) == C_FLONUM_TAG) {
1862      C_word *a = C_alloc(C_SIZEOF_FLONUM);
1863      C_kontinue(k, C_flonum(&a, (double)C_unfix(x) * C_flonum_magnitude(y)));
1864    } else if (C_IS_BIGNUM_TYPE(y)) {
1865      CPS_CALL(C_u_2_integer_times, 4, (C_word)NULL, k, x, y);
1866    } else {
1867      try_extended_number("numbers#@extended-2-times", 3, k, x, y);
1868    }
1869  } else if (C_immediatep(x)) {
1870    barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", x);
1871  } else if (C_block_header(x) == C_FLONUM_TAG) {
1872    C_word *a = C_alloc(C_SIZEOF_FLONUM);
1873    if (y & C_FIXNUM_BIT) {
1874      C_kontinue(k, C_flonum(&a, C_flonum_magnitude(x) * (double)C_unfix(y)));
1875    } else if (C_immediatep(y)) {
1876      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", y);
1877    } else if (C_block_header(y) == C_FLONUM_TAG) {
1878      C_kontinue(k, C_a_i_flonum_times(&a, 2, x, y));
1879    } else if (C_IS_BIGNUM_TYPE(y)) {
1880      C_kontinue(k, C_flonum(&a, C_flonum_magnitude(x)*C_bignum_to_double(y)));
1881    } else {
1882      try_extended_number("numbers#@extended-2-times", 3, k, x, y);
1883    }
1884  } else if (C_IS_BIGNUM_TYPE(x)) {
1885    if (y & C_FIXNUM_BIT) {
1886      CPS_CALL(C_u_2_integer_times, 4, (C_word)NULL, k, x, y);
1887    } else if (C_immediatep(y)) {
1888      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, "*", x);
1889    } else if (C_block_header(y) == C_FLONUM_TAG) {
1890      C_word *a = C_alloc(C_SIZEOF_FLONUM);
1891      C_kontinue(k, C_flonum(&a, C_bignum_to_double(x)*C_flonum_magnitude(y)));
1892    } else if (C_IS_BIGNUM_TYPE(y)) {
1893      CPS_CALL(C_u_2_integer_times, 4, (C_word)NULL, k, x, y);
1894    } else {
1895      try_extended_number("numbers#@extended-2-times", 3, k, x, y);
1896    }
1897  } else {
1898    try_extended_number("numbers#@extended-2-times", 3, k, x, y);
1899  }
1900}
1901
1902/* Needs SIZEOF_BIGNUM(2) */
1903C_regparm C_word C_fcall
1904C_a_u_i_2_fixnum_times(C_word **ptr, C_word n, C_word x, C_word y)
1905{
1906  C_uword negp, xhi, xlo, yhi, ylo, p, rhi, rlo;
1907
1908  negp = ((x & C_INT_SIGN_BIT) ? !(y & C_INT_SIGN_BIT) : (y & C_INT_SIGN_BIT));
1909  x = (x & C_INT_SIGN_BIT) ? -C_unfix(x) : C_unfix(x);
1910  y = (y & C_INT_SIGN_BIT) ? -C_unfix(y) : C_unfix(y);
1911
1912  xhi = C_BIGNUM_DIGIT_HI_HALF(x); xlo = C_BIGNUM_DIGIT_LO_HALF(x);
1913  yhi = C_BIGNUM_DIGIT_HI_HALF(y); ylo = C_BIGNUM_DIGIT_LO_HALF(y);
1914 
1915  /* This is simply bignum_digits_multiply unrolled for 2x2 halfdigits */
1916  p = xlo * ylo;
1917  rlo = C_BIGNUM_DIGIT_LO_HALF(p);
1918
1919  p = xhi * ylo + C_BIGNUM_DIGIT_HI_HALF(p);
1920  rhi = C_BIGNUM_DIGIT_HI_HALF(p);
1921
1922  p = xlo * yhi + C_BIGNUM_DIGIT_LO_HALF(p);
1923  rlo = C_BIGNUM_DIGIT_COMBINE(C_BIGNUM_DIGIT_LO_HALF(p), rlo);
1924
1925  rhi = xhi * yhi + C_BIGNUM_DIGIT_HI_HALF(p) + rhi;
1926
1927  if (rhi) {
1928    return C_bignum2(ptr, negp != 0, rlo, rhi);
1929  } else if (negp ?
1930             ((rlo & C_INT_SIGN_BIT) || !C_fitsinfixnump(-(C_word)rlo)) :
1931             !C_ufitsinfixnump(rlo)) {
1932    return C_bignum1(ptr, negp != 0, rlo);
1933  } else {
1934    return C_fix(negp ? -rlo : rlo);
1935  }
1936}
1937
1938void C_ccall CPS_PROC2(C_u_2_integer_times, c, self, k, x, y)
1939{
1940  CPS_BODY2(c, self, k, x, y);
1941  if (x & C_FIXNUM_BIT) {
1942    if (y & C_FIXNUM_BIT) {
1943      C_word *a = C_alloc(C_SIZEOF_BIGNUM(2));
1944      C_kontinue(k, C_a_u_i_2_fixnum_times(&a, 2, x, y));
1945    } else {
1946      C_word tmp = x; /* swap to ensure x is a bignum and y a fixnum */
1947      x = y;
1948      y = tmp;
1949    }
1950  }
1951  /* Here, we know for sure that X is a bignum */
1952  if (y == C_fix(0)) {
1953    C_kontinue(k, C_fix(0));
1954  } else if (y == C_fix(1)) {
1955    C_kontinue(k, x);
1956  } else if (y == C_fix(-1)) {
1957    CPS_CALL(C_u_integer_negate, 3, (C_word)NULL, k, x);
1958  } else if (y & C_FIXNUM_BIT) { /* Any other fixnum */
1959    C_word absy = (y & C_INT_SIGN_BIT) ? -C_unfix(y) : C_unfix(y),
1960           negp = C_mk_bool((y & C_INT_SIGN_BIT) ?
1961                            !C_bignum_negativep(x) :
1962                            C_bignum_negativep(x));
1963 
1964    if (C_fitsinbignumhalfdigitp(absy) ||
1965        (((C_uword)1 << (C_ilen(absy)-1)) == absy && C_fitsinfixnump(absy))) {
1966      C_word size, k2, *a = C_alloc(C_SIZEOF_CLOSURE(4));
1967      k2 = C_closure(&a, 4, (C_word)integer_times_2, k, x, C_fix(absy));
1968      size = C_fix(C_bignum_size(x) + 1); /* Needs _at most_ one more digit */
1969      C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
1970    } else {
1971      C_word *a = C_alloc(C_SIZEOF_FIX_BIGNUM);
1972      y = C_a_u_i_fix_to_big(&a, y);
1973      bignum_times_bignum_unsigned(k, x, y, negp);
1974    }
1975  } else {
1976    C_word negp = C_bignum_negativep(x) ?
1977                  !C_bignum_negativep(y) :
1978                  C_bignum_negativep(y);
1979    bignum_times_bignum_unsigned(k, x, y, C_mk_bool(negp));
1980  }
1981}
1982
1983static void CONT_PROC(integer_times_2, c, self, new_big)
1984{
1985  CONT_BODY(self, new_big);
1986  C_word k = C_block_item(self, 1),
1987         old_bigx = C_block_item(self, 2),
1988         absy = C_unfix(C_block_item(self, 3));
1989  C_uword *digits = C_bignum_digits(new_big),
1990          *end_digit = digits + C_bignum_size(old_bigx);
1991  int shift;
1992
1993  bignum_digits_destructive_copy(new_big, old_bigx);
1994
1995  /* Scale up, and sanitise the result. */
1996  shift = C_ilen(absy) - 1;
1997  if (((C_uword)1 << shift) == absy) { /* Power of two? */
1998    *end_digit = bignum_digits_destructive_shift_left(digits, end_digit, shift);
1999  } else {
2000    *end_digit =
2001      bignum_digits_destructive_scale_up_with_carry(digits, end_digit, absy, 0);
2002  }
2003  C_kontinue(k, C_bignum_simplify(new_big));
2004}
2005
2006static void
2007bignum_times_bignum_unsigned(C_word k, C_word x, C_word y, C_word negp)
2008{
2009  if (C_bignum_size(y) < C_bignum_size(x)) { /* Ensure size(x) <= size(y) */
2010    C_word z = x;
2011    x = y;
2012    y = z;
2013  }
2014
2015  if (C_bignum_size(x) < C_KARATSUBA_THRESHOLD) {  /* Gradebook */
2016    C_word kab[C_SIZEOF_CLOSURE(4)], *ka = kab, k2, size;
2017    k2 = C_closure(&ka, 4, (C_word)bignum_times_bignum_unsigned_2, k, x, y);
2018    size = C_fix(C_bignum_size(x) + C_bignum_size(y));
2019    C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_TRUE);
2020  } else {
2021    try_extended_number("numbers#@bignum-2-times-karatsuba", 3, k, x, y);
2022  }
2023}
2024
2025static C_regparm void
2026bignum_digits_multiply(C_word x, C_word y, C_word result)
2027{
2028  C_uword product,
2029          *xd = C_bignum_digits(x),
2030          *yd = C_bignum_digits(y),
2031          *rd = C_bignum_digits(result);
2032  C_uhword carry, yj;
2033  /* Lengths in halfwords */
2034  int i, j, length_x = C_bignum_size(x) * 2, length_y = C_bignum_size(y) * 2;
2035
2036  /* From Hacker's Delight, Figure 8-1 (top part) */
2037  for (j = 0; j < length_y; ++j) {
2038    yj = C_uhword_ref(yd, j);
2039    if (yj == 0) continue;
2040    carry = 0;
2041    for (i = 0; i < length_x; ++i) {
2042      product = (C_uword)C_uhword_ref(xd, i) * yj +
2043                (C_uword)C_uhword_ref(rd, i + j) + carry;
2044      C_uhword_set(rd, i + j, product);
2045      carry = C_BIGNUM_DIGIT_HI_HALF(product);
2046    }
2047    C_uhword_set(rd, j + length_x, carry);
2048  }
2049}
2050
2051static void CONT_PROC(bignum_times_bignum_unsigned_2, c, self, result)
2052{
2053  CONT_BODY(self, result);
2054  C_word k = C_block_item(self, 1),
2055         x = C_block_item(self, 2),
2056         y = C_block_item(self, 3);
2057
2058  bignum_digits_multiply(x, y, result);
2059  C_kontinue(k, C_bignum_simplify(result));
2060}
2061
2062void C_ccall
2063CPS_PROC5(C_digits_to_integer, c, self, k, str, start, end, radix, negp)
2064{
2065  CPS_BODY5(c, self, k, str, start, end, radix, negp);
2066  assert((C_unfix(radix) > 1) && C_fitsinbignumhalfdigitp(C_unfix(radix)));
2067 
2068  if (start == end) {
2069    C_kontinue(k, C_SCHEME_FALSE);
2070  } else {
2071    C_word kab[C_SIZEOF_CLOSURE(6)], *ka = kab, k2, size;
2072    size_t nbits;
2073    k2 = C_closure(&ka, 6, (C_word)digits_to_integer_2, k, str, start, end, radix);
2074
2075    nbits = (C_unfix(end) - C_unfix(start)) * C_ilen(C_unfix(radix)-1);
2076    size = C_fix(C_BIGNUM_BITS_TO_DIGITS(nbits));
2077    C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
2078  }
2079}
2080
2081C_inline int hex_char_to_digit(int ch)
2082{
2083  if (ch == (int)'#') return 0; /* Hash characters in numbers are mapped to 0 */
2084  else if (ch >= (int)'a') return ch - (int)'a' + 10; /* lower hex */
2085  else if (ch >= (int)'A') return ch - (int)'A' + 10; /* upper hex */
2086  else return ch - (int)'0'; /* decimal (OR INVALID; handled elsewhere) */
2087}
2088
2089static void CONT_PROC(digits_to_integer_2, c, self, result)
2090{
2091  CONT_BODY(self, result);
2092  C_word k = C_block_item(self, 1),
2093         str = C_block_item(self, 2),
2094         start = C_unfix(C_block_item(self, 3)),
2095         end = C_unfix(C_block_item(self, 4)),
2096         radix = C_unfix(C_block_item(self, 5));
2097  char *s = C_c_string(str);
2098
2099  C_kontinue(k, str_to_bignum(result, s + start, s + end, radix));
2100}
2101
2102/* Write from digit character stream to bignum.  Bignum does not need
2103 * to be initialised.  Returns the bignum, or a fixnum.  Assumes the
2104 * string contains only digits that fit within radix (checked by
2105 * string->number).
2106 */
2107static C_regparm C_word
2108str_to_bignum(C_word bignum, char *str, char *str_end, int radix)
2109{
2110  int radix_shift, str_digit;
2111  C_uword *digits = C_bignum_digits(bignum),
2112          *end_digits = digits + C_bignum_size(bignum), big_digit = 0;
2113
2114  /* Below, we try to save up as much as possible in big_digit, and
2115   * only when it exceeds what we would be able to multiply easily, we
2116   * scale up the bignum and add what we saved up.
2117   */
2118  radix_shift = C_ilen(radix) - 1;
2119  if (((C_uword)1 << radix_shift) == radix) { /* Power of two? */
2120    int n = 0; /* Number of bits read so far into current big digit */
2121
2122    /* Read from least to most significant digit to avoid shifting or scaling */
2123    while (str_end > str) {
2124      str_digit = hex_char_to_digit((int)*--str_end);
2125
2126      big_digit |= (C_uword)str_digit << n;
2127      n += radix_shift;
2128
2129      if (n >= C_BIGNUM_DIGIT_LENGTH) {
2130        n -= C_BIGNUM_DIGIT_LENGTH;
2131        *digits++ = big_digit;
2132        big_digit = str_digit >> (radix_shift - n);
2133      }
2134    }
2135    assert(n < C_BIGNUM_DIGIT_LENGTH);
2136    /* If radix isn't an exact divisor of digit length, write final digit */
2137    if (n > 0) *digits++ = big_digit;
2138    assert(digits == end_digits);
2139  } else {                        /* Not a power of two */
2140    C_uword *last_digit = digits, factor;  /* bignum starts as zero */
2141
2142    do {
2143      factor = radix;
2144      while (str < str_end && C_fitsinbignumhalfdigitp(factor)) {
2145        str_digit = hex_char_to_digit((int)*str++);
2146        factor *= radix;
2147        big_digit = radix * big_digit + str_digit;
2148      }
2149
2150      big_digit = bignum_digits_destructive_scale_up_with_carry(
2151                   digits, last_digit, factor / radix, big_digit);
2152
2153      if (big_digit) {
2154        (*last_digit++) = big_digit; /* Move end */
2155        big_digit = 0;
2156      }
2157    } while (str < str_end);
2158
2159    /* Set remaining digits to zero so bignum_simplify can do its work */
2160    assert(last_digit <= end_digits);
2161    while (last_digit < end_digits) *last_digit++ = 0;
2162  }
2163
2164  return C_bignum_simplify(bignum);
2165}
2166
2167/* TODO: Copied from runtime.c */
2168# define STRING_BUFFER_SIZE   4096
2169
2170static C_TLS C_char buffer[ STRING_BUFFER_SIZE ];
2171static char *to_n_nary(C_uword num, C_uword base, int negp, int as_flonum)
2172{
2173  static char *digits = "0123456789abcdef";
2174  char *p;
2175  C_uword shift = C_ilen(base) - 1;
2176  int mask = (1 << shift) - 1;
2177  if (as_flonum) {
2178    buffer[68] = '\0';
2179    buffer[67] = '0';
2180    buffer[66] = '.';
2181  } else {
2182    buffer[66] = '\0';
2183  }
2184  p = buffer + 66;
2185  if (mask == base - 1) {
2186    do {
2187      *(--p) = digits [ num & mask ];
2188      num >>= shift;
2189    } while (num);
2190  } else {
2191    do {
2192      *(--p) = digits [ num % base ];
2193      num /= base;
2194    } while (num);
2195  }
2196  if (negp) *(--p) = '-';
2197  return p;
2198}
2199
2200void C_ccall CPS_PROCNPLUS1(C_basic_number_to_string, c, closure, k, num)
2201{
2202  CPS_BODYNPLUS1(c, closure, k, num);
2203  C_word radix;
2204
2205  if(c == 3) {
2206    radix = C_fix(10);
2207  } else if(c == 4) {
2208#ifdef ARGVECTOR_CHICKEN
2209    radix = __av[3];
2210#else
2211    va_list v;
2212
2213    va_start(v, num);
2214    radix = va_arg(v, C_word);
2215    va_end(v);
2216#endif
2217   
2218    if(!(radix & C_FIXNUM_BIT))
2219      barf(C_BAD_ARGUMENT_TYPE_BAD_BASE_ERROR, "number->string", radix);
2220  } else {
2221    C_bad_argc(c, 3);
2222  }
2223
2224  if(num & C_FIXNUM_BIT) {
2225    CPS_CALL(C_u_fixnum_to_string, 4, (C_word)NULL, k, num, radix);
2226  } else if (C_immediatep(num)) {
2227    barf(C_BAD_ARGUMENT_TYPE_ERROR, "number->string", num);
2228  } else if(C_block_header(num) == C_FLONUM_TAG) {
2229    CPS_CALL(C_u_flonum_to_string, 4, (C_word)NULL, k, num, radix);
2230  } else if (C_IS_BIGNUM_TYPE(num)) {
2231    CPS_CALL(C_u_integer_to_string, 4, (C_word)NULL, k, num, radix);
2232  } else {
2233    try_extended_number("numbers#@extended-number->string", 3, k, num, radix);
2234  }
2235}
2236
2237void C_ccall CPS_PROC2(C_u_fixnum_to_string, c, self, k, num, radix)
2238{
2239  CPS_BODY2(c, self, k, num, radix);
2240  C_char *p;
2241  C_word *a, neg = num & C_INT_SIGN_BIT ? 1 : 0;
2242
2243  radix = C_unfix(radix);
2244  if (radix < 2 || radix > 16) {
2245    barf(C_BAD_ARGUMENT_TYPE_BAD_BASE_ERROR, "number->string", radix);
2246  }
2247
2248  num = neg ? -C_unfix(num) : C_unfix(num);
2249  p = to_n_nary(num, radix, neg, 0);
2250
2251  num = C_strlen(p);
2252  a = C_alloc((C_bytestowords(num) + 1));
2253  C_kontinue(k, C_string(&a, num, p));
2254}
2255
2256void C_ccall CPS_PROC2(C_u_flonum_to_string, c, self, k, num, radix)
2257{
2258  CPS_BODY2(c, self, k, num, radix);
2259  C_word *a;
2260  C_char *p;
2261  double f;
2262
2263  radix = C_unfix(radix);
2264  f = C_flonum_magnitude(num);
2265
2266  /* XXX TODO: Should inexacts be printable in other bases than 10?
2267   * Perhaps output a string starting with #i?
2268   * Right now something like (number->string 1e40 16) results in
2269   * a string that can't be read back using string->number.
2270   */
2271  if((radix < 2) || (radix > 16)){
2272    barf(C_BAD_ARGUMENT_TYPE_BAD_BASE_ERROR, "number->string", C_fix(radix));
2273  }
2274
2275  if(C_fits_in_unsigned_int_p(num) == C_SCHEME_TRUE) { /* Use fast int code */
2276    if(f < 0) {
2277      p = to_n_nary((C_uword)-f, radix, 1, 1);
2278    } else {
2279      p = to_n_nary((C_uword)f, radix, 0, 1);
2280    }
2281  } else if(C_isnan(f)) {
2282    p = "+nan.0";
2283  } else if(C_isinf(f)) {
2284    p = f > 0 ? "+inf.0" : "-inf.0";
2285  } else { /* Doesn't fit an unsigned int and not "special"; use system libc */
2286    C_snprintf(buffer, STRING_BUFFER_SIZE, C_text("%.*g"),
2287               /* XXX: flonum_print_precision */
2288               (int)C_unfix(C_get_print_precision()), f);
2289    buffer[STRING_BUFFER_SIZE-1] = '\0';
2290
2291    if((p = C_strpbrk(buffer, C_text(".eE"))) == NULL) {
2292      /* Already checked for these, so shouldn't happen */
2293      assert(*buffer != 'i'); /* "inf" */
2294      assert(*buffer != 'n'); /* "nan" */
2295      /* Ensure integral flonums w/o expt are always terminated by .0 */
2296#if defined(HAVE_STRLCAT) || !defined(C_strcat)
2297      C_strlcat(buffer, C_text(".0"), sizeof(buffer));
2298#else
2299      C_strcat(buffer, C_text(".0"));
2300#endif
2301    }
2302    p = buffer;
2303  }
2304
2305  radix = C_strlen(p);
2306  a = C_alloc((C_bytestowords(radix) + 1));
2307  radix = C_string(&a, radix, p);
2308  C_kontinue(k, radix);
2309}
2310
2311/* Naming is a little inconsistent, but looks saner.  We're not R-O-B-O-T-S! */
2312void C_ccall CPS_PROC2(C_u_integer_to_string, c, self, k, num, radix)
2313{
2314  CPS_BODY2(c, self, k, num, radix);
2315  if (num & C_FIXNUM_BIT) {
2316    CPS_CALL(C_u_fixnum_to_string, 4, (C_word)NULL, k, num, radix);
2317  } else {
2318    int len, radix_shift;
2319    size_t nbits;
2320
2321    if ((C_unfix(radix) < 2) || (C_unfix(radix) > 16)) {
2322      barf(C_BAD_ARGUMENT_TYPE_BAD_BASE_ERROR, "number->string", radix);
2323    }
2324
2325    /* Approximation of the number of radix digits we'll need.  We try
2326     * to be as precise as possible to avoid memmove overhead at the end
2327     * of the non-powers of two part of the conversion procedure, which
2328     * we may need to do because we write strings back-to-front, and
2329     * pointers must be aligned (even for byte blocks).
2330     */
2331    len = C_bignum_size(num)-1;
2332
2333    nbits  = (size_t)len * C_BIGNUM_DIGIT_LENGTH;
2334    nbits += C_ilen(C_bignum_digits(num)[len]);
2335
2336    len = C_ilen(C_unfix(radix))-1;
2337    len = (nbits + len - 1) / len;
2338    len += C_bignum_negativep(num) ? 1 : 0; /* Add space for negative sign */
2339   
2340    radix_shift = C_ilen(C_unfix(radix)) - 1;
2341    if (len > C_RECURSIVE_TO_STRING_THRESHOLD &&
2342        /* The power of two fast path is much faster than recursion */
2343        ((C_uword)1 << radix_shift) != C_unfix(radix)) {
2344      try_extended_number("numbers#@integer->string/recursive",
2345                          4, k, num, radix, C_fix(len));
2346    } else {
2347#ifdef ARGVECTOR_CHICKEN
2348      C_word kab[C_SIZEOF_CLOSURE(4)], *ka = kab, av[6];
2349
2350      av[ 0 ] = (C_word)NULL;   /* No "self" closure */
2351      av[ 1 ] = C_closure(&ka, 4, (C_word)bignum_to_str_2, k, num, radix);
2352      av[ 2 ] = C_fix(len);
2353      av[ 3 ] = C_SCHEME_TRUE; /* Byte vector */
2354      av[ 4 ] = C_SCHEME_FALSE; /* No initialization */
2355      av[ 5 ] = C_SCHEME_FALSE; /* Don't align at 8 bytes */
2356      C_allocate_vector(6, av);
2357#else
2358      C_word k2, *ka;
2359      ka = C_alloc(C_SIZEOF_CLOSURE(4));
2360      k2 = C_closure(&ka, 4, (C_word)bignum_to_str_2, k, num, radix);
2361      C_allocate_vector(6, (C_word)NULL, k2, C_fix(len),
2362                        /* Byte vec, no initialization, no align at 8 bytes */
2363                        C_SCHEME_TRUE, C_SCHEME_FALSE, C_SCHEME_FALSE);
2364#endif
2365    }
2366  }
2367}
2368
2369static void CONT_PROC(bignum_to_str_2, c, self, string)
2370{
2371  static char *characters = "0123456789abcdef";
2372  CONT_BODY(self, string);
2373  C_word k = C_block_item(self, 1),
2374         bignum = C_block_item(self, 2),
2375         radix = C_unfix(C_block_item(self, 3));
2376  char *buf = C_c_string(string), *index = buf + C_header_size(string) - 1;
2377  int radix_shift, negp = (C_bignum_negativep(bignum) ? 1 : 0);
2378
2379  radix_shift = C_ilen(radix) - 1;
2380  if (((C_uword)1 << radix_shift) == radix) { /* Power of two? */
2381    int radix_mask = radix - 1, big_digit_len = 0, radix_digit;
2382    C_uword *scan, *end, big_digit = 0;
2383
2384    scan = C_bignum_digits(bignum);
2385    end = scan + C_bignum_size(bignum);
2386
2387    while (scan < end) {
2388      /* If radix isn't an exact divisor of digit length, handle overlap */
2389      if (big_digit_len == 0) {
2390        big_digit = *scan++;
2391        big_digit_len = C_BIGNUM_DIGIT_LENGTH;
2392      } else {
2393        assert(index >= buf);
2394        radix_digit = big_digit;
2395        big_digit = *scan++;
2396        radix_digit |= (big_digit << big_digit_len) & radix_mask;
2397        *index-- = characters[radix_digit];
2398        big_digit >>= (radix_shift - big_digit_len);
2399        big_digit_len = C_BIGNUM_DIGIT_LENGTH - (radix_shift - big_digit_len);
2400      }
2401
2402      while(big_digit_len >= radix_shift && index >= buf) {
2403        radix_digit = big_digit & radix_mask;
2404        *index-- = characters[radix_digit];
2405        big_digit >>= radix_shift;
2406        big_digit_len -= radix_shift;
2407      }
2408    }
2409
2410    assert(big_digit < radix);
2411
2412    /* Final digit (like overlap at start of while loop) */
2413    if (big_digit) *index-- = characters[big_digit];
2414
2415    if (negp) {
2416      /* Loop above might've overwritten sign position with a zero */
2417      if (*(index+1) == '0') *(index+1) = '-';
2418      else *index-- = '-';
2419    }
2420
2421    /* Length calculation is always precise for radix powers of two. */
2422    assert(index == buf-1);
2423  } else {
2424    C_uword base, *start, *scan, big_digit;
2425    C_word working_copy;
2426    int steps, i;
2427
2428    working_copy = allocate_tmp_bignum(C_fix(C_bignum_size(bignum)),
2429                                       C_mk_bool(negp), C_SCHEME_FALSE);
2430    bignum_digits_destructive_copy(working_copy, bignum);
2431
2432    start = C_bignum_digits(working_copy);
2433
2434    scan = start + C_bignum_size(bignum);
2435    /* Calculate the largest power of radix that fits a halfdigit:
2436     * steps = log10(2^halfdigit_bits), base = 10^steps
2437     */
2438    for(steps = 0, base = radix; C_fitsinbignumhalfdigitp(base); base *= radix)
2439      steps++;
2440
2441    base /= radix; /* Back down: we overshot in the loop */
2442
2443    while (scan > start) {
2444      big_digit = bignum_digits_destructive_scale_down(start, scan, base);
2445
2446      if (*(scan-1) == 0) scan--; /* Adjust if we exhausted the highest digit */
2447
2448      for(i = 0; i < steps && index >= buf; ++i) {
2449        C_word tmp = big_digit / radix;
2450        *index-- = characters[big_digit - (tmp*radix)]; /* big_digit % radix */
2451        big_digit = tmp;
2452      }
2453    }
2454    assert(index >= buf-1);
2455    free_tmp_bignum(working_copy);
2456
2457    /* Move index onto first nonzero digit.  We're writing a bignum
2458       here: it can't consist of only zeroes. */
2459    while(*++index == '0');
2460 
2461    if (negp) *--index = '-';
2462 
2463    /* Shorten with distance between start and index. */
2464    if (buf != index) {
2465      i = C_header_size(string) - (index - buf);
2466      C_memmove(buf, index, i); /* Move start of number to beginning. */
2467      C_block_header(string) = C_STRING_TYPE | i; /* Mutate strlength. */
2468    }
2469  }
2470
2471  C_kontinue(k, string);
2472}
2473
2474C_regparm double C_bignum_to_double(C_word bignum)
2475{
2476  double accumulator = 0;
2477  C_uword *start = C_bignum_digits(bignum),
2478          *scan = start + C_bignum_size(bignum);
2479  while (start < scan) {
2480    accumulator *= (C_word)1 << C_BIGNUM_HALF_DIGIT_LENGTH;
2481    accumulator *= (C_word)1 << C_BIGNUM_HALF_DIGIT_LENGTH;
2482    accumulator += (*--scan);
2483  }
2484  return(C_bignum_negativep(bignum) ? -accumulator : accumulator);
2485}
2486
2487static void
2488fabs_frexp_to_digits(C_uword exp, double sign, C_uword *start, C_uword *scan)
2489{
2490  C_uword digit, odd_bits = exp % C_BIGNUM_DIGIT_LENGTH;
2491
2492  assert(C_isfinite(sign));
2493  assert(0.5 <= sign && sign < 1); /* Guaranteed by frexp() and fabs() */
2494  assert((scan - start) == C_BIGNUM_BITS_TO_DIGITS(exp));
2495 
2496  if (odd_bits > 0) { /* Handle most significant digit first */
2497    sign *= (C_uword)1 << odd_bits;
2498    digit = (C_uword)sign;
2499    (*--scan) = digit;
2500    sign -= (double)digit;
2501  }
2502
2503  while (start < scan && sign > 0) {
2504    sign *= pow(2.0, C_BIGNUM_DIGIT_LENGTH);
2505    digit = (C_uword)sign;
2506    (*--scan) = digit;
2507    sign -= (double)digit;
2508  }
2509
2510  /* Finish up by clearing any remaining, lower, digits */
2511  while (start < scan)
2512    (*--scan) = 0;
2513}
2514
2515static C_word
2516flo_to_tmp_bignum(C_word x)
2517{
2518  /* TODO: allocating and initialising the bignum is pointless if we
2519   * already know the number of limbs in the comparand.  In fact,
2520   * bignum_cmp will first check the number of limbs and *then*
2521   * compare.  Instead, we can check beforehand and check the limbs
2522   * directly against the generated limbs, without allocating at all!
2523   */
2524  C_word tmp_big, negp = C_mk_bool(C_flonum_magnitude(x) < 0.0);
2525  int exponent;
2526  double significand = frexp(C_flonum_magnitude(x), &exponent);
2527
2528  assert(C_u_i_fpintegerp_fixed(x));
2529
2530  if (exponent <= 0) {
2531    tmp_big = allocate_tmp_bignum(C_fix(0), C_SCHEME_FALSE, C_SCHEME_FALSE);
2532  } else if (exponent == 1) { /* TODO: check significand * 2^exp fits fixnum? */
2533    /* Don't use fix_to_big to simplify caller code: it can just free this */
2534    tmp_big = allocate_tmp_bignum(C_fix(1), negp, C_SCHEME_FALSE);
2535    C_bignum_digits(tmp_big)[0] = 1;
2536  } else {
2537    C_uword size, *start, *end;
2538
2539    size = C_fix(C_BIGNUM_BITS_TO_DIGITS(exponent));
2540
2541    tmp_big = allocate_tmp_bignum(size, negp, C_SCHEME_FALSE);
2542    start = C_bignum_digits(tmp_big);
2543    end = start + C_bignum_size(tmp_big);
2544
2545    fabs_frexp_to_digits(exponent, fabs(significand), start, end);
2546  }
2547  return tmp_big;
2548}
2549
2550void C_ccall CPS_PROC1(C_u_flo_to_int, c, self, k, x)
2551{
2552  CPS_BODY1(c, self, k, x);
2553  int exponent;
2554  double significand = frexp(C_flonum_magnitude(x), &exponent);
2555
2556  assert(C_truep(C_u_i_fpintegerp_fixed(x)));
2557
2558  if (exponent <= 0) {
2559    C_kontinue(k, C_fix(0));
2560  } else if (exponent == 1) { /* TODO: check significand * 2^exp fits fixnum? */
2561    C_kontinue(k, significand < 0.0 ? C_fix(-1) : C_fix(1));
2562  } else {
2563    C_word kab[C_SIZEOF_CLOSURE(4) + C_SIZEOF_FLONUM], *ka = kab, k2, size,
2564           negp = C_mk_bool(C_flonum_magnitude(x) < 0.0),
2565           sign = C_flonum(&ka, fabs(significand));
2566
2567    k2 = C_closure(&ka, 4, (C_word)flo_to_int_2, k, C_fix(exponent), sign);
2568
2569    size = C_fix(C_BIGNUM_BITS_TO_DIGITS(exponent));
2570    C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
2571  }
2572}
2573
2574static void CONT_PROC(flo_to_int_2, c, self, result)
2575{
2576  CONT_BODY(self, result);
2577  C_word k = C_block_item(self, 1);
2578  C_uword exponent = C_unfix(C_block_item(self, 2)),
2579          *start = C_bignum_digits(result),
2580          *scan = start + C_bignum_size(result);
2581  double significand = C_flonum_magnitude(C_block_item(self, 3));
2582
2583  fabs_frexp_to_digits(exponent, significand, start, scan);
2584  C_kontinue(k, C_bignum_simplify(result));
2585}
2586
2587C_word C_ccall
2588C_u_i_integer_length(C_word x)
2589{
2590  if (x & C_FIXNUM_BIT) {
2591    return C_u_i_fixnum_length(x);
2592  } else {
2593    C_uword result = (C_bignum_size(x) - 1) * C_BIGNUM_DIGIT_LENGTH,
2594            *last_digit = C_bignum_digits(x) + C_bignum_size(x) - 1,
2595            last_digit_length = C_ilen(*last_digit);
2596
2597    /* If *only* the highest bit is set, negating will give one less bit */
2598    if (C_bignum_negativep(x) &&
2599        *last_digit == ((C_uword)1 << (last_digit_length-1))) {
2600      C_uword *startx = C_bignum_digits(x);
2601      while (startx < last_digit && *startx == 0) ++startx;
2602      if (startx == last_digit) result--;
2603    }
2604    return C_fix(result + last_digit_length);
2605  }
2606}
2607
2608/* x is any exact integer but y is _always_ a fixnum */
2609void C_ccall CPS_PROC2(C_u_integer_shift, c, self, k, x, y)
2610{
2611  CPS_BODY2(c, self, k, x, y);
2612  C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab;
2613
2614  y = C_unfix(y);
2615  if (y == 0 || x == C_fix(0)) { /* Done (no shift) */
2616    C_kontinue(k, x);
2617  } else if (x & C_FIXNUM_BIT) {
2618    if (y < 0) {
2619      /* Don't shift more than a word's length (that's undefined in C!) */
2620      if (-y < C_WORD_SIZE) {
2621        C_kontinue(k, C_fix(C_unfix(x) >> -y));
2622      } else {
2623        C_kontinue(k, (x < 0) ? C_fix(-1) : C_fix(0));
2624      }
2625    } else if (y > 0 && y < C_WORD_SIZE-2 &&
2626               /* After shifting, the length still fits a fixnum */
2627               (C_ilen(C_unfix(x)) + y) < C_WORD_SIZE-2) {
2628      C_kontinue(k, C_fix(C_unfix(x) << y));
2629    } else {
2630      x = C_a_u_i_fix_to_big(&a, x);
2631    }
2632  }
2633
2634  {
2635    C_word ab[C_SIZEOF_CLOSURE(6)], *a = ab,
2636           k2, size, negp, digit_offset, bit_offset;
2637
2638    negp = C_mk_bool(C_bignum_negativep(x));
2639 
2640    if (y > 0) { /* y is guaranteed not to be 0 here */
2641      digit_offset = y / C_BIGNUM_DIGIT_LENGTH;
2642      bit_offset =   y % C_BIGNUM_DIGIT_LENGTH;
2643
2644      k2 = C_closure(&a, 6, (C_word)bignum_actual_shift, k,
2645                     x, C_SCHEME_TRUE, C_fix(digit_offset), C_fix(bit_offset));
2646      size = C_fix(C_bignum_size(x) + digit_offset + 1);
2647      C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
2648    } else if (-y >= C_bignum_size(x) * (C_word)C_BIGNUM_DIGIT_LENGTH) {
2649      /* All bits are shifted out, just return 0 or -1 */
2650      C_kontinue(k, C_truep(negp) ? C_fix(-1) : C_fix(0));
2651    } else {
2652      digit_offset = -y / C_BIGNUM_DIGIT_LENGTH;
2653      bit_offset =   -y % C_BIGNUM_DIGIT_LENGTH;
2654   
2655      k2 = C_closure(&a, 6, (C_word)bignum_actual_shift, k,
2656                     x, C_SCHEME_FALSE, C_fix(digit_offset), C_fix(bit_offset));
2657
2658      size = C_fix(C_bignum_size(x) - digit_offset);
2659      C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
2660    }
2661  }
2662}
2663
2664C_inline C_word maybe_negate_bignum_for_bitwise_op(C_word x, C_word size)
2665{
2666  C_word nx = C_SCHEME_FALSE, xsize;
2667  if (C_bignum_negativep(x)) {
2668    nx = allocate_tmp_bignum(C_fix(size), C_SCHEME_FALSE, C_SCHEME_FALSE);
2669    xsize = C_bignum_size(x);
2670    /* Copy up until requested size, and init any remaining upper digits */
2671    C_memcpy(C_bignum_digits(nx), C_bignum_digits(x),
2672             C_wordstobytes(nmin(size, xsize)));
2673    if (size > xsize)
2674      C_memset(C_bignum_digits(nx)+xsize, 0, C_wordstobytes(size-xsize));
2675    bignum_digits_destructive_negate(nx);
2676  }
2677  return nx;
2678}
2679
2680static void CONT_PROC(bignum_actual_shift, c, self, result)
2681{
2682  CONT_BODY(self, result);
2683  C_word k = C_block_item(self, 1),
2684         x = C_block_item(self, 2),
2685         shift_left = C_truep(C_block_item(self, 3)),
2686         digit_offset = C_unfix(C_block_item(self, 4)),
2687         bit_offset = C_unfix(C_block_item(self, 5));
2688  C_uword *startr = C_bignum_digits(result),
2689          *startx = C_bignum_digits(x),
2690          *endx = startx + C_bignum_size(x),
2691          *endr = startr + C_bignum_size(result);
2692
2693  if (shift_left) {
2694    /* Initialize only the lower digits we're skipping and the MSD */
2695    C_memset(startr, 0, C_wordstobytes(digit_offset));
2696    *(endr-1) = 0;
2697    startr += digit_offset;
2698    /* Can't use bignum_digits_destructive_copy because it assumes
2699     * we want to copy from the start.
2700     */
2701    C_memcpy(startr, startx, C_wordstobytes(endx-startx));
2702    if(bit_offset > 0)
2703      bignum_digits_destructive_shift_left(startr, endr, bit_offset);
2704  } else {
2705    C_word nx, size = C_bignum_size(x) + 1;
2706    if (C_truep(nx = maybe_negate_bignum_for_bitwise_op(x, size))) {
2707      startx = C_bignum_digits(nx); /* update startx; x and endx are unused */
2708    }
2709
2710    startx += digit_offset;
2711    /* Can't use bignum_digits_destructive_copy because that assumes
2712     * target is at least as big as source.
2713     */
2714    C_memcpy(startr, startx, C_wordstobytes(endr-startr));
2715    if(bit_offset > 0)
2716      bignum_digits_destructive_shift_right(startr,endr,bit_offset,C_truep(nx));
2717
2718    if (C_truep(nx)) {
2719      free_tmp_bignum(nx);
2720      bignum_digits_destructive_negate(result);
2721    }
2722  }
2723  C_kontinue(k, C_bignum_simplify(result));
2724}
2725
2726/* This is currently only used by Karatsuba multiplication and
2727 * Burnikel-Ziegler division.  It is not intended as a public API!
2728 */
2729void C_ccall
2730CPS_PROC3(C_u_bignum_extract_digits, c, self, k, x, start, end)
2731{
2732  CPS_BODY3(c, self, k, x, start, end);
2733  if (x & C_FIXNUM_BIT) { /* Needed? */
2734    if (C_unfix(start) == 0 && (end == C_SCHEME_FALSE || C_unfix(end) > 0)) {
2735      C_kontinue(k, x);
2736    } else {
2737      C_kontinue(k, C_fix(0));
2738    }
2739  } else {
2740    C_word negp, size;
2741
2742    negp = C_mk_bool(C_bignum_negativep(x)); /* Always false */
2743
2744    start = C_unfix(start);
2745    /* We might get passed larger values than actually fits; pad w/ zeroes */
2746    if (end == C_SCHEME_FALSE) end = C_bignum_size(x);
2747    else end = nmin(C_unfix(end), C_bignum_size(x));
2748    assert(start >= 0);
2749
2750    size = end - start;
2751
2752    if (size == 0 || start >= C_bignum_size(x)) {
2753      C_kontinue(k, C_fix(0));
2754    } else {
2755      C_word k2, kab[C_SIZEOF_CLOSURE(5)], *ka = kab;
2756      k2 = C_closure(&ka, 5, (C_word)bignum_actual_extraction,
2757                     k, x, C_fix(start), C_fix(end));
2758      C_allocate_bignum(5, (C_word)NULL, k2, C_fix(size), negp, C_SCHEME_FALSE);
2759    }
2760  }
2761}
2762
2763static void CONT_PROC(bignum_actual_extraction, c, self, result)
2764{
2765  CONT_BODY(self, result);
2766  C_word k = C_block_item(self, 1),
2767         x = C_block_item(self, 2),
2768         start = C_unfix(C_block_item(self, 3)),
2769         end = C_unfix(C_block_item(self, 4));
2770  C_uword *result_digits = C_bignum_digits(result),
2771          *x_digits = C_bignum_digits(x);
2772
2773  /* Can't use bignum_digits_destructive_copy because that assumes
2774   * target is at least as big as source.
2775   */
2776  C_memcpy(result_digits, x_digits + start, C_wordstobytes(end-start));
2777  C_kontinue(k, C_bignum_simplify(result));
2778}
2779
2780C_regparm C_word C_ccall C_u_i_integer_randomize(C_word seed)
2781{
2782  /* TODO: Rename C_randomize to C_u_i_fixnum_randomize */
2783  if (seed & C_FIXNUM_BIT) {
2784    return C_randomize(seed);
2785  } else {
2786    /*
2787     * This random number generator is very simple. Probably too simple...
2788     */
2789    C_uword *scan = C_bignum_digits(seed),
2790            *end = scan + C_bignum_size(seed), iseed = 0;
2791
2792    /* What a cheap way to initialize the random generator. I feel dirty! */
2793    while (scan < end)
2794      iseed ^= *scan++;
2795
2796    srand(iseed);
2797    return C_SCHEME_UNDEFINED;
2798  }
2799}
2800
2801void C_ccall CPS_PROC1(C_u_integer_random, c, self, k, max)
2802{
2803  CPS_BODY1(c, self, k, max);
2804  /* TODO: for consistency C_random_fixnum should be called C_u_i_fixnum_random */
2805  if (max & C_FIXNUM_BIT) {
2806    C_kontinue(k, C_random_fixnum(max));
2807  } else {
2808    C_word k2, kab[C_SIZEOF_CLOSURE(4)], *ka = kab, size,
2809           max_len, max_bits, max_top_digit, d, negp;
2810
2811    max_len = C_bignum_size(max);
2812    max_top_digit = d = C_bignum_digits(max)[max_len - 1];
2813 
2814    max_bits = (max_len - 1) * C_BIGNUM_DIGIT_LENGTH;
2815    while(d) {
2816      max_bits++;
2817      d >>= 1;
2818    }
2819    /* Subtract/add one because we don't want zero to be over-represented */
2820    size = ((double)rand())/(RAND_MAX + 1.0) * (double)(max_bits - 1);
2821    size = C_fix(C_BIGNUM_BITS_TO_DIGITS(size + 1));
2822
2823    negp = C_mk_bool(C_bignum_negativep(max));
2824    k2 = C_closure(&ka, 4, (C_word)bignum_random_2, k, C_fix(max_top_digit), C_fix(max_len));
2825    C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
2826  }
2827}
2828
2829static void CONT_PROC(bignum_random_2, c, self, result)
2830{
2831  CONT_BODY(self, result);
2832  C_word k = C_block_item(self, 1),
2833         max_top_digit = C_unfix(C_block_item(self, 2)),
2834         max_len = C_unfix(C_block_item(self, 3));
2835  C_uword *scan = C_bignum_digits(result),
2836          *end = scan + C_bignum_size(result); /* Go to just before the end. */
2837
2838  while(scan < end)
2839    *scan++ = ((double)rand())/(RAND_MAX + 1.0) * pow(2.0, C_BIGNUM_DIGIT_LENGTH);
2840  /*
2841   * Last word is special when length is max_len: It must be less than
2842   * max's most significant digit, instead of 2^{digitlen}.
2843   */
2844  if (max_len == C_bignum_size(result))
2845    *scan = ((double)rand())/(RAND_MAX + 1.0) * (double)max_top_digit;
2846  else
2847    *scan = ((double)rand())/(RAND_MAX + 1.0) * pow(2.0, C_BIGNUM_DIGIT_LENGTH);
2848
2849  C_kontinue(k, C_bignum_simplify(result));
2850}
2851
2852C_word C_ccall
2853C_u_i_integer_bit_setp(C_word n, C_word i)
2854{
2855  if (!(i & C_FIXNUM_BIT)) { /* A bit silly, but might be useful */
2856    return C_u_i_integer_negativep(n);
2857  } else if (i & C_INT_SIGN_BIT) {
2858    barf(C_BAD_ARGUMENT_TYPE_NO_UINTEGER_ERROR, "bit-set?", n, i);
2859  } else {
2860    i = C_unfix(i);
2861    if (n & C_FIXNUM_BIT) {
2862      if (i >= C_WORD_SIZE) return C_mk_bool(n & C_INT_SIGN_BIT);
2863      else return C_mk_bool((C_unfix(n) & ((C_word)1 << i)) != 0);
2864    } else {
2865      C_word nn, d;
2866      d = i / C_BIGNUM_DIGIT_LENGTH;
2867      if (d >= C_bignum_size(n)) return C_mk_bool(C_bignum_negativep(n));
2868
2869      if (C_truep(nn = maybe_negate_bignum_for_bitwise_op(n, d))) n = nn;
2870
2871      i %= C_BIGNUM_DIGIT_LENGTH;
2872      d = C_mk_bool((C_bignum_digits(n)[d] & (C_uword)1 << i) != 0);
2873      if (C_truep(nn)) free_tmp_bignum(nn);
2874      return d;
2875    }
2876  }
2877}
2878
2879void C_ccall CPS_PROC2(C_u_2_integer_bitwise_and, c, self, k, x, y)
2880{
2881  CPS_BODY2(c, self, k, x, y);
2882  if ((x & y) & C_FIXNUM_BIT) {
2883    C_kontinue(k, C_u_fixnum_and(x, y));
2884  } else {
2885    C_word kab[C_SIZEOF_FIX_BIGNUM*2], *ka = kab, negp, size, k2;
2886    if (x & C_FIXNUM_BIT) x = C_a_u_i_fix_to_big(&ka, x);
2887    if (y & C_FIXNUM_BIT) y = C_a_u_i_fix_to_big(&ka, y);
2888
2889    negp = C_mk_bool(C_bignum_negativep(x) && C_bignum_negativep(y));
2890    /* Allow negative 1-bits to propagate */
2891    if (C_bignum_negativep(x) || C_bignum_negativep(y))
2892      size = nmax(C_bignum_size(x), C_bignum_size(y)) + 1;
2893    else
2894      size = nmin(C_bignum_size(x), C_bignum_size(y));
2895
2896    ka = C_alloc(C_SIZEOF_CLOSURE(4)); /* Why faster than static alloc? */
2897    k2 = C_closure(&ka, 4, (C_word)bignum_bitwise_and_2, k, x, y);
2898    C_allocate_bignum(5, (C_word)NULL, k2, C_fix(size), negp, C_SCHEME_FALSE);
2899  }
2900}
2901
2902static void CONT_PROC(bignum_bitwise_and_2, c, self, result)
2903{
2904  CONT_BODY(self, result);
2905  C_word k = C_block_item(self, 1),
2906         x = C_block_item(self, 2),
2907         y = C_block_item(self, 3),
2908         size = C_bignum_size(result), nx, ny;
2909  C_uword *scanr = C_bignum_digits(result),
2910          *endr = scanr + C_bignum_size(result),
2911          *scans1, *ends1, *scans2;
2912
2913  if (C_truep(nx = maybe_negate_bignum_for_bitwise_op(x, size))) x = nx;
2914  if (C_truep(ny = maybe_negate_bignum_for_bitwise_op(y, size))) y = ny;
2915
2916  if (C_bignum_size(x) < C_bignum_size(y)) {
2917    scans1 = C_bignum_digits(x); ends1 = scans1 + C_bignum_size(x);
2918    scans2 = C_bignum_digits(y);
2919  } else {
2920    scans1 = C_bignum_digits(y); ends1 = scans1 + C_bignum_size(y);
2921    scans2 = C_bignum_digits(x);
2922  }
2923
2924  while (scans1 < ends1) *scanr++ = *scans1++ & *scans2++;
2925  C_memset(scanr, 0, C_wordstobytes(endr - scanr));
2926
2927  if (C_truep(nx)) free_tmp_bignum(nx);
2928  if (C_truep(ny)) free_tmp_bignum(ny);
2929  if (C_bignum_negativep(result)) bignum_digits_destructive_negate(result);
2930
2931  C_kontinue(k, C_bignum_simplify(result));
2932}
2933
2934void C_ccall CPS_PROC2(C_u_2_integer_bitwise_ior, c, self, k, x, y)
2935{
2936  CPS_BODY2(c, self, k, x, y);
2937  if ((x & y) & C_FIXNUM_BIT) {
2938    C_kontinue(k, C_u_fixnum_or(x, y));
2939  } else {
2940    C_word kab[C_SIZEOF_FIX_BIGNUM*2], *ka = kab, negp, size, k2;
2941    if (x & C_FIXNUM_BIT) x = C_a_u_i_fix_to_big(&ka, x);
2942    if (y & C_FIXNUM_BIT) y = C_a_u_i_fix_to_big(&ka, y);
2943
2944    ka = C_alloc(C_SIZEOF_CLOSURE(4)); /* Why faster than static alloc? */
2945    k2 = C_closure(&ka, 4, (C_word)bignum_bitwise_ior_2, k, x, y);
2946    size = C_fix(nmax(C_bignum_size(x), C_bignum_size(y)) + 1);
2947    negp = C_mk_bool(C_bignum_negativep(x) || C_bignum_negativep(y));
2948    C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
2949  }
2950}
2951
2952static void CONT_PROC(bignum_bitwise_ior_2, c, self, result)
2953{
2954  CONT_BODY(self, result);
2955  C_word k = C_block_item(self, 1),
2956         x = C_block_item(self, 2),
2957         y = C_block_item(self, 3),
2958         size = C_bignum_size(result), nx, ny;
2959  C_uword *scanr = C_bignum_digits(result),
2960          *endr = scanr + C_bignum_size(result),
2961          *scans1, *ends1, *scans2, *ends2;
2962
2963  if (C_truep(nx = maybe_negate_bignum_for_bitwise_op(x, size))) x = nx;
2964  if (C_truep(ny = maybe_negate_bignum_for_bitwise_op(y, size))) y = ny;
2965
2966  if (C_bignum_size(x) < C_bignum_size(y)) {
2967    scans1 = C_bignum_digits(x); ends1 = scans1 + C_bignum_size(x);
2968    scans2 = C_bignum_digits(y); ends2 = scans2 + C_bignum_size(y);
2969  } else {
2970    scans1 = C_bignum_digits(y); ends1 = scans1 + C_bignum_size(y);
2971    scans2 = C_bignum_digits(x); ends2 = scans2 + C_bignum_size(x);
2972  }
2973
2974  while (scans1 < ends1) *scanr++ = *scans1++ | *scans2++;
2975  while (scans2 < ends2) *scanr++ = *scans2++;
2976  if (scanr < endr) *scanr++ = 0; /* Only done when result is positive */
2977  assert(scanr == endr);
2978
2979  if (C_truep(nx)) free_tmp_bignum(nx);
2980  if (C_truep(ny)) free_tmp_bignum(ny);
2981  if (C_bignum_negativep(result)) bignum_digits_destructive_negate(result);
2982
2983  C_kontinue(k, C_bignum_simplify(result));
2984}
2985
2986void C_ccall CPS_PROC2(C_u_2_integer_bitwise_xor, c, self, k, x, y)
2987{
2988  CPS_BODY2(c, self, k, x, y);
2989  if ((x & y) & C_FIXNUM_BIT) {
2990    C_kontinue(k, C_fixnum_xor(x, y));
2991  } else {
2992    C_word kab[C_SIZEOF_FIX_BIGNUM*2], *ka = kab, size, k2, negp;
2993    if (x & C_FIXNUM_BIT) x = C_a_u_i_fix_to_big(&ka, x);
2994    if (y & C_FIXNUM_BIT) y = C_a_u_i_fix_to_big(&ka, y);
2995
2996    ka = C_alloc(C_SIZEOF_CLOSURE(4)); /* Why faster than static alloc? */
2997    k2 = C_closure(&ka, 4, (C_word)bignum_bitwise_xor_2, k, x, y);
2998    size = C_fix(nmax(C_bignum_size(x), C_bignum_size(y)) + 1);
2999    negp = C_mk_bool(C_bignum_negativep(x) != C_bignum_negativep(y));
3000    C_allocate_bignum(5, (C_word)NULL, k2, size, negp, C_SCHEME_FALSE);
3001  }
3002}
3003
3004static void CONT_PROC(bignum_bitwise_xor_2, c, self, result)
3005{
3006  CONT_BODY(self, result);
3007  C_word k = C_block_item(self, 1),
3008         x = C_block_item(self, 2),
3009         y = C_block_item(self, 3),
3010         size = C_bignum_size(result), nx, ny;
3011  C_uword *scanr = C_bignum_digits(result),
3012          *endr = scanr + C_bignum_size(result),
3013          *scans1, *ends1, *scans2, *ends2;
3014
3015  if (C_truep(nx = maybe_negate_bignum_for_bitwise_op(x, size))) x = nx;
3016  if (C_truep(ny = maybe_negate_bignum_for_bitwise_op(y, size))) y = ny;
3017
3018  if (C_bignum_size(x) < C_bignum_size(y)) {
3019    scans1 = C_bignum_digits(x); ends1 = scans1 + C_bignum_size(x);
3020    scans2 = C_bignum_digits(y); ends2 = scans2 + C_bignum_size(y);
3021  } else {
3022    scans1 = C_bignum_digits(y); ends1 = scans1 + C_bignum_size(y);
3023    scans2 = C_bignum_digits(x); ends2 = scans2 + C_bignum_size(x);
3024  }
3025
3026  while (scans1 < ends1) *scanr++ = *scans1++ ^ *scans2++;
3027  while (scans2 < ends2) *scanr++ = *scans2++;
3028  if (scanr < endr) *scanr++ = 0; /* Only done when result is positive */
3029  assert(scanr == endr);
3030
3031  if (C_truep(nx)) free_tmp_bignum(nx);
3032  if (C_truep(ny)) free_tmp_bignum(ny);
3033  if (C_bignum_negativep(result)) bignum_digits_destructive_negate(result);
3034
3035  C_kontinue(k, C_bignum_simplify(result));
3036}
3037
3038static void bignum_digits_destructive_negate(C_word result)
3039{
3040  C_uword *scan, *end, digit, sum;
3041
3042  scan = C_bignum_digits(result);
3043  end = scan + C_bignum_size(result);
3044
3045  do {
3046    digit = ~*scan;
3047    sum = digit + 1;
3048    *scan++ = sum;
3049  } while (sum == 0 && scan < end);
3050
3051  for (; scan < end; scan++) {
3052    *scan = ~*scan;
3053  }
3054}
3055
3056/* This is ugly but really cleans up the code below */
3057#define RETURN_Q_AND_OR_R(calc_q, calc_r)                 \
3058  if (C_truep(C_and(return_q, return_r))) {               \
3059    C_kontinue2(k, calc_q, calc_r);                       \
3060  } else if (C_truep(return_r)) {                         \
3061    C_kontinue(k, calc_r);                                \
3062  } else {                                                \
3063    C_kontinue(k, calc_q);                                \
3064  }
3065
3066/* Lossy; we could be in "quotient&remainder" or "modulo" */
3067#define DIVREM_LOC ((C_truep(C_and(return_q, return_r))) ? "/" :        \
3068                    (C_truep(return_q) ? "quotient" : "remainder"))
3069
3070static C_regparm void
3071basic_divrem(C_word k, C_word x, C_word y, C_word return_q, C_word return_r)
3072{
3073  if (x & C_FIXNUM_BIT) {
3074    if (y & C_FIXNUM_BIT) {
3075      C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab;
3076      if (y == C_fix(0)) numbers_div_by_zero_error(DIVREM_LOC);
3077
3078      RETURN_Q_AND_OR_R(C_a_u_i_fixnum_quotient_checked(&a, 2, x, y),
3079                        C_u_i_fixnum_remainder_checked(x, y));
3080    } else if (C_immediatep(y)) {
3081      barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, DIVREM_LOC, y);
3082    } else if (C_block_header(y) == C_FLONUM_TAG) {
3083      C_word ab[C_SIZEOF_FLONUM*3], *a = ab;
3084      if (C_flonum_magnitude(y) == 0.0) numbers_div_by_zero_error(DIVREM_LOC);
3085
3086      x = C_a_i_fix_to_flo(&a, 1, x);
3087      RETURN_Q_AND_OR_R(C_a_i_flonum_actual_quotient_checked(&a, 2, x, y),
3088                        C_a_i_flonum_remainder_checked(&a, 2, x, y));
3089    } else if (C_IS_BIGNUM_TYPE(y)) {
3090      integer_divrem(k, x, y, return_q, return_r);
3091    } else {
3092      barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, DIVREM_LOC, y);
3093    }
3094  } else if (C_immediatep(x)) {
3095    barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, DIVREM_LOC, x);
3096  } else if (C_block_header(x) == C_FLONUM_TAG) {
3097    if (!C_truep(C_u_i_fpintegerp_fixed(x))) {
3098      barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, DIVREM_LOC, x);
3099    } else if (y & C_FIXNUM_BIT) {
3100      C_word ab[C_SIZEOF_FLONUM*3], *a = ab;
3101      if (y == C_fix(0)) numbers_div_by_zero_error(DIVREM_LOC);
3102
3103      y = C_a_i_fix_to_flo(&a, 1, y);
3104      RETURN_Q_AND_OR_R(C_a_i_flonum_actual_quotient_checked(&a, 2, x, y),
3105                        C_a_i_flonum_remainder_checked(&a, 2, x, y));
3106    } else if (C_immediatep(y)) {
3107      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, DIVREM_LOC, y);
3108    } else if (C_block_header(y) == C_FLONUM_TAG) {
3109      C_word ab[C_SIZEOF_FLONUM*2], *a = ab;
3110      if (C_flonum_magnitude(y) == 0.0) numbers_div_by_zero_error(DIVREM_LOC);
3111
3112      RETURN_Q_AND_OR_R(C_a_i_flonum_actual_quotient_checked(&a, 2, x, y),
3113                        C_a_i_flonum_remainder_checked(&a, 2, x, y));
3114    } else if (C_IS_BIGNUM_TYPE(y)) {
3115      C_word k2, ab[C_SIZEOF_CLOSURE(3)], *a = ab;
3116      x = flo_to_tmp_bignum(x);
3117      k2 = C_closure(&a, 3, (C_word)divrem_intflo_2, k, x);
3118      integer_divrem(k2, x, y, return_q, return_r);
3119    } else {
3120      barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, DIVREM_LOC, y);
3121    }
3122  } else if (C_IS_BIGNUM_TYPE(x)) {
3123    if (y & C_FIXNUM_BIT) {
3124      integer_divrem(k, x, y, return_q, return_r);
3125    } else if (C_immediatep(y)) {
3126      barf(C_BAD_ARGUMENT_TYPE_NO_NUMBER_ERROR, DIVREM_LOC, y);
3127    } else if (C_block_header(y) == C_FLONUM_TAG) {
3128      if (!C_truep(C_u_i_fpintegerp_fixed(y))) {
3129        barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, DIVREM_LOC, y);
3130      } else if (C_flonum_magnitude(y) == 0.0) {
3131        numbers_div_by_zero_error(DIVREM_LOC);
3132      } else {
3133        C_word k2, ab[C_SIZEOF_CLOSURE(3)], *a = ab;
3134        y = flo_to_tmp_bignum(y);
3135        k2 = C_closure(&a, 3, (C_word)divrem_intflo_2, k, y);
3136        integer_divrem(k2, x, y, return_q, return_r);
3137      }
3138    } else if (C_IS_BIGNUM_TYPE(y)) {
3139      bignum_divrem(k, x, y, return_q, return_r);
3140    } else {
3141      barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, DIVREM_LOC, y);
3142    }
3143  } else {
3144    barf(C_BAD_ARGUMENT_TYPE_NO_INTEGER_ERROR, DIVREM_LOC, x);
3145  }
3146}
3147
3148static void CONT_PROCN(divrem_intflo_2, c, self)
3149{
3150  CONT_BODYN(c, self);
3151  C_word k = C_block_item(self, 1), x, y;
3152#ifndef ARGVECTOR_CHICKEN
3153  va_list v;
3154#endif
3155
3156  free_tmp_bignum(C_block_item(self, 2));
3157 
3158#ifndef ARGVECTOR_CHICKEN
3159  va_start(v, self);
3160#endif
3161  if (c == 2) {
3162    C_word ab[C_SIZEOF_FLONUM], *a = ab;
3163#ifdef ARGVECTOR_CHICKEN
3164    x = __av[1];
3165#else
3166    x = va_arg(v, C_word);
3167    va_end(v);
3168#endif
3169    if (x & C_FIXNUM_BIT) x = C_a_i_fix_to_flo(&a, 1, x);
3170    else x = C_a_u_i_big_to_flo(&a, 1, x);
3171    C_kontinue(k, x);
3172  } else { /* c == 3 */
3173    C_word ab[C_SIZEOF_FLONUM*2], *a = ab;
3174#ifdef ARGVECTOR_CHICKEN
3175    x = __av[1];
3176    y = __av[2];
3177#else
3178    x = va_arg(v, C_word);
3179    y = va_arg(v, C_word);
3180    va_end(v);
3181#endif
3182   
3183    if (x & C_FIXNUM_BIT) x = C_a_i_fix_to_flo(&a, 1, x);
3184    else x = C_a_u_i_big_to_flo(&a, 1, x);
3185    if (y & C_FIXNUM_BIT) y = C_a_i_fix_to_flo(&a, 1, y);
3186    else y = C_a_u_i_big_to_flo(&a, 1, y);
3187    C_kontinue2(k, x, y);
3188  }
3189}
3190
3191static void CONT_PROC(bignum_divrem_fixnum_2, c, self, negated_big)
3192{
3193   CONT_BODY(self, negated_big);
3194   C_word k = C_block_item(self, 1),
3195          return_q = C_block_item(self, 2),
3196          return_r = C_block_item(self, 3);
3197   RETURN_Q_AND_OR_R(negated_big, C_fix(0));
3198}
3199
3200static C_regparm void
3201integer_divrem(C_word k, C_word x, C_word y, C_word return_q, C_word return_r)
3202{
3203  if (!(y & C_FIXNUM_BIT)) { /* y is bignum. */
3204    if (x & C_FIXNUM_BIT) {
3205      /* abs(x) < abs(y), so it will always be [0, x] except for this case: */
3206      if (x == C_fix(C_MOST_NEGATIVE_FIXNUM) &&
3207          C_bignum_negated_fitsinfixnump(y)) {
3208        RETURN_Q_AND_OR_R(C_fix(-1), C_fix(0));
3209      } else {
3210        RETURN_Q_AND_OR_R(C_fix(0), x);
3211      }
3212    } else {
3213      bignum_divrem(k, x, y, return_q, return_r);
3214    }
3215  } else if (x & C_FIXNUM_BIT) { /* both x and y are fixnum. */
3216    C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab;
3217    if (y == C_fix(0)) numbers_div_by_zero_error(DIVREM_LOC);
3218
3219    RETURN_Q_AND_OR_R(C_a_u_i_fixnum_quotient_checked(&a, 2, x, y),
3220                      C_u_i_fixnum_remainder_checked(x, y));
3221  } else { /* x is bignum, y is fixnum. */
3222    C_word absy = (y & C_INT_SIGN_BIT) ? -C_unfix(y) : C_unfix(y);
3223
3224    if (y == C_fix(1)) {
3225      RETURN_Q_AND_OR_R(x, C_fix(0));
3226    } else if (y == C_fix(0)) {
3227      numbers_div_by_zero_error(DIVREM_LOC);
3228    } else if (y == C_fix(-1)) {
3229      C_word *ka, k2;
3230      ka = C_alloc(C_SIZEOF_CLOSURE(4));
3231      k2 = C_closure(&ka, 4, (C_word)bignum_divrem_fixnum_2,
3232                     k, return_q, return_r);
3233      CPS_CALL(C_u_integer_negate, 3, (C_word)NULL, k2, x);
3234    } else if (C_fitsinbignumhalfdigitp(absy) ||
3235               ((((C_uword)1 << (C_ilen(absy)-1)) == absy) &&
3236                C_fitsinfixnump(absy))) {
3237      if (C_truep(return_q)) {
3238        C_word q_negp = C_mk_bool((y & C_INT_SIGN_BIT) ?
3239                                  !(C_bignum_negativep(x)) :
3240                                  C_bignum_negativep(x)),
3241                r_negp = C_mk_bool(C_bignum_negativep(x)),
3242                *ka, k2, size;
3243        ka = C_alloc(C_SIZEOF_CLOSURE(7));
3244        size = C_fix(C_bignum_size(x));
3245        k2 = C_closure(&ka, 7,
3246                       (C_word)bignum_destructive_divide_unsigned_small,
3247                       k, x, C_fix(absy),
3248                       return_q, return_r, r_negp);
3249        C_allocate_bignum(5, (C_word)NULL, k2, size, q_negp, C_SCHEME_FALSE);
3250      } else {
3251        C_word rem;
3252        C_uword next_power = (C_uword)1 << (C_ilen(absy)-1);
3253
3254        if (next_power == absy) { /* Is absy a power of two? */
3255          rem = *(C_bignum_digits(x)) & (next_power - 1);
3256        } else { /* Too bad, we have to do some real work */
3257          rem = bignum_remainder_unsigned_halfdigit(x, absy);
3258        }
3259        C_kontinue(k, C_bignum_negativep(x) ? C_fix(-rem) : C_fix(rem));
3260      }
3261    } else {                    /* Just divide it as two bignums */
3262      C_word ab[C_SIZEOF_FIX_BIGNUM], *a = ab;
3263      bignum_divrem(k, x, C_a_u_i_fix_to_big(&a, y), return_q, return_r);
3264    }
3265  }
3266}
3267
3268static C_regparm void
3269bignum_divrem(C_word k, C_word x, C_word y, C_word return_q, C_word return_r)
3270{
3271  C_word q_negp = C_mk_bool(C_bignum_negativep(y) ?
3272                            !C_bignum_negativep(x) :
3273                            C_bignum_negativep(x)),
3274         r_negp = C_mk_bool(C_bignum_negativep(x)), size;
3275
3276  switch(bignum_cmp_unsigned(x, y)) {
3277  case 0:
3278    RETURN_Q_AND_OR_R(C_truep(q_negp) ? C_fix(-1) : C_fix(1), C_fix(0));
3279  case -1:
3280    RETURN_Q_AND_OR_R(C_fix(0), x);
3281  case 1:
3282  default:
3283    size = C_bignum_size(x) - C_bignum_size(y);
3284    if (C_bignum_size(y) > C_BURNIKEL_ZIEGLER_THRESHOLD &&
3285        size > C_BURNIKEL_ZIEGLER_THRESHOLD) {
3286      try_extended_number("numbers#@bignum-2-divrem-burnikel-ziegler", 5,
3287                          k, x, y, return_q, return_r);
3288    } else if (C_truep(return_q)) {
3289      C_word kab[C_SIZEOF_CLOSURE(6)], *ka = kab, k2;
3290      k2 = C_closure(&ka, 6, (C_word)bignum_divide_2_unsigned, k,
3291                     x, y, return_r, r_negp);
3292      size = C_fix(C_bignum_size(x) + 1 - C_bignum_size(y));
3293      C_allocate_bignum(5, (C_word)NULL, k2, size, q_negp, C_SCHEME_FALSE);
3294    } else { /* We can skip bignum_divide_2_unsigned if we need no quotient */
3295      C_word kab[C_SIZEOF_CLOSURE(7)], *ka = kab, k2;
3296      k2 = C_closure(&ka, 7, (C_word)bignum_divide_2_unsigned_2, k,
3297                     x, y, return_q, return_r, C_SCHEME_UNDEFINED);
3298      size = C_fix(C_bignum_size(x) + 1); /* May need to be normalized */
3299      C_allocate_bignum(5, (C_word)NULL, k2, size, r_negp, C_SCHEME_FALSE);
3300    }
3301  }
3302}
3303
3304static C_word
3305bignum_remainder_unsigned_halfdigit(C_word num, C_word den)
3306{
3307  C_uword *start = C_bignum_digits(num),
3308          *scan = start + C_bignum_size(num),
3309          rem = 0, two_digits;
3310
3311  assert((den > 1) && (C_fitsinbignumhalfdigitp(den)));
3312  while (start < scan) {
3313    two_digits = (*--scan);
3314    rem = C_BIGNUM_DIGIT_COMBINE(rem, C_BIGNUM_DIGIT_HI_HALF(two_digits)) % den;
3315    rem = C_BIGNUM_DIGIT_COMBINE(rem, C_BIGNUM_DIGIT_LO_HALF(two_digits)) % den;
3316  }
3317  return rem;
3318}
3319
3320/* External interface for the above internal divrem functions */
3321void C_ccall CPS_PROC2(C_basic_divrem, c, self, k, x, y)
3322{
3323  CPS_BODY2(c, self, k, x, y);
3324  if (c != 4) C_bad_argc_2(c, 4, self);
3325  basic_divrem(k, x, y, C_SCHEME_TRUE, C_SCHEME_TRUE);
3326}
3327
3328void C_ccall CPS_PROC2(C_u_integer_divrem, c, self, k, x, y)
3329{
3330  CPS_BODY2(c, self, k, x, y);
3331  integer_divrem(k, x, y, C_SCHEME_TRUE, C_SCHEME_TRUE);
3332}
3333
3334void C_ccall CPS_PROC2(C_basic_remainder, c, self, k, x, y)
3335{
3336  CPS_BODY2(c, self, k, x, y);
3337  if (c != 4) C_bad_argc_2(c, 4, self);
3338  basic_divrem(k, x, y, C_SCHEME_FALSE, C_SCHEME_TRUE);
3339}
3340
3341void C_ccall CPS_PROC2(C_u_integer_remainder, c, self, k, x, y)
3342{
3343  CPS_BODY2(c, self, k, x, y);
3344  integer_divrem(k, x, y, C_SCHEME_FALSE, C_SCHEME_TRUE);
3345}
3346
3347void C_ccall CPS_PROC2(C_basic_quotient, c, self, k, x, y)
3348{
3349  CPS_BODY2(c, self, k, x, y);
3350  if (c != 4) C_bad_argc_2(c, 4, self);
3351  basic_divrem(k, x, y, C_SCHEME_TRUE, C_SCHEME_FALSE);
3352}
3353
3354void C_ccall CPS_PROC2(C_u_integer_quotient, c, self, k, x, y)
3355{
3356  CPS_BODY2(c, self, k, x, y);
3357  integer_divrem(k, x, y, C_SCHEME_TRUE, C_SCHEME_FALSE);
3358}
3359
3360/* "small" is either a number that fits a halfdigit, or a power of two */
3361static void CONT_PROC(bignum_destructive_divide_unsigned_small,
3362                      c, self, quotient)
3363{
3364  CONT_BODY(self, quotient);
3365  C_word k = C_block_item(self, 1),
3366         numerator = C_block_item(self, 2),
3367         denominator = C_unfix(C_block_item(self, 3)),
3368         /* return_quotient = C_block_item(self, 4), */
3369         return_remainder = C_block_item(self, 5),
3370         remainder_negp = C_block_item(self, 6);
3371  C_uword *start = C_bignum_digits(quotient),
3372          *end = start + C_bignum_size(quotient),
3373          remainder;
3374  int shift_amount;
3375
3376  bignum_digits_destructive_copy(quotient, numerator);
3377
3378  shift_amount = C_ilen(denominator)-1;
3379  if (((C_uword)1 << shift_amount) == denominator) { /* Power of two?  Shift! */
3380    remainder = bignum_digits_destructive_shift_right(start,end,shift_amount,0);
3381    assert(C_ufitsinfixnump(remainder));
3382  } else {
3383    remainder = bignum_digits_destructive_scale_down(start, end, denominator);
3384    assert(C_fitsinbignumhalfdigitp(remainder));
3385  }
3386
3387  quotient = C_bignum_simplify(quotient);
3388
3389  if (C_truep(return_remainder)) {
3390    remainder = C_truep(remainder_negp) ? -remainder : remainder;
3391    C_kontinue2(k, quotient, C_fix(remainder));
3392  } else {
3393    C_kontinue(k, quotient);
3394  }
3395}
3396
3397
3398/* Full bignum division */
3399
3400static void CONT_PROC(bignum_divide_2_unsigned, c, self, quotient)
3401{
3402  CONT_BODY(self, quotient);
3403  C_word k = C_block_item(self, 1),
3404         x = C_block_item(self, 2),
3405         y = C_block_item(self, 3),
3406         size = C_fix(C_bignum_size(x) + 1),
3407         return_r = C_block_item(self, 4),
3408         r_negp = C_block_item(self, 5),
3409         kab[C_SIZEOF_CLOSURE(7)], *ka = kab, k2;
3410
3411  k2 = C_closure(&ka, 7, (C_word)bignum_divide_2_unsigned_2, k,
3412                 x, y, C_SCHEME_TRUE, return_r, quotient);
3413  C_allocate_bignum(5, (C_word)NULL, k2, size, r_negp, C_SCHEME_FALSE);
3414}
3415
3416/* For help understanding this algorithm, see:
3417   Knuth, Donald E., "The Art of Computer Programming",
3418   volume 2, "Seminumerical Algorithms"
3419   section 4.3.1, "Multiple-Precision Arithmetic".
3420
3421   [Yeah, that's a nice book but that particular section is not
3422   helpful at all, which is also pointed out by P. Brinch Hansen's
3423   "Multiple-Length Division Revisited: A Tour Of The Minefield".
3424   That's a more down-to-earth step-by-step explanation of the
3425   algorithm.  Add to this the C implementation in Hacker's Delight
3426   (section 9-2, p141--142) and you may be able to grok this...
3427   ...barely, if you're as math-challenged as I am -- sjamaan]
3428*/
3429
3430static void CONT_PROC(bignum_divide_2_unsigned_2, c, self, remainder)
3431{
3432  CONT_BODY(self, remainder);
3433  C_word k = C_block_item(self, 1),
3434         numerator = C_block_item(self, 2),
3435         denominator = C_block_item(self, 3),
3436         return_quotient = C_block_item(self, 4),
3437         return_remainder = C_block_item(self, 5),
3438         quotient = C_block_item(self, 6),
3439         length = C_bignum_size(denominator);
3440  C_uword d1 = *(C_bignum_digits(denominator) + length - 1),
3441          *startr = C_bignum_digits(remainder),
3442          *endr = startr + C_bignum_size(remainder);
3443  int shift;
3444
3445  shift = C_BIGNUM_DIGIT_LENGTH - C_ilen(d1); /* nlz */
3446
3447  /* We have to work on halfdigits, so we shift out only the necessary
3448   * amount in order fill out that halfdigit (base is halved).
3449   * This trick is shamelessly stolen from Gauche :)
3450   * See below for part 2 of the trick.
3451   */
3452  if (shift >= C_BIGNUM_HALF_DIGIT_LENGTH)
3453    shift -= C_BIGNUM_HALF_DIGIT_LENGTH;
3454
3455  /* Code below won't always set high halfdigit of quotient, so do it here. */
3456  if (quotient != C_SCHEME_UNDEFINED)
3457    C_bignum_digits(quotient)[C_bignum_size(quotient)-1] = 0;
3458
3459  bignum_digits_destructive_copy(remainder, numerator);
3460  *(endr-1) = 0;    /* Ensure most significant digit is initialised */
3461  if (shift == 0) { /* Already normalized */
3462    bignum_destructive_divide_normalized(remainder, denominator, quotient);
3463  } else { /* Requires normalisation; allocate scratch denominator for this */
3464    C_uword *startnd;
3465    C_word ndenom;
3466
3467    bignum_digits_destructive_shift_left(startr, endr, shift);
3468
3469    ndenom = allocate_tmp_bignum(C_fix(length), C_SCHEME_FALSE, C_SCHEME_FALSE);
3470    startnd = C_bignum_digits(ndenom);
3471    bignum_digits_destructive_copy(ndenom, denominator);
3472    bignum_digits_destructive_shift_left(startnd, startnd+length, shift);
3473
3474    bignum_destructive_divide_normalized(remainder, ndenom, quotient);
3475    if (C_truep(return_remainder)) /* Otherwise, don't bother shifting back */
3476      bignum_digits_destructive_shift_right(startr, endr, shift, 0);
3477
3478    free_tmp_bignum(ndenom);
3479  }
3480
3481  if (C_truep(return_remainder)) {
3482    if (C_truep(return_quotient)) {
3483      C_kontinue2(k, C_bignum_simplify(quotient),
3484                  C_bignum_simplify(remainder));
3485    } else {
3486      C_kontinue(k, C_bignum_simplify(remainder));
3487    }
3488  } else {
3489    assert(C_truep(return_quotient));
3490    C_kontinue(k, C_bignum_simplify(quotient));
3491  }
3492}
3493
3494static void
3495bignum_destructive_divide_normalized(C_word big_u, C_word big_v, C_word big_q)
3496{
3497  C_uword *v = C_bignum_digits(big_v),
3498          *u = C_bignum_digits(big_u),
3499          *q = big_q == C_SCHEME_UNDEFINED ? NULL : C_bignum_digits(big_q),
3500           p,               /* product of estimated quotient & "denominator" */
3501           hat, qhat, rhat, /* estimated quotient and remainder digit */
3502           vn_1, vn_2;      /* "cached" values v[n-1], v[n-2] */
3503  C_word t, k;              /* Two helpers: temp/final remainder and "borrow" */
3504  /* We use plain ints here, which theoretically may not be enough on
3505   * 64-bit for an insanely huge number, but it is a _lot_ faster.
3506   */
3507  int n = C_bignum_size(big_v) * 2,       /* in halfwords */
3508      m = (C_bignum_size(big_u) * 2) - 2; /* Correct for extra digit */
3509  int i, j;                /* loop  vars */
3510
3511  /* Part 2 of Gauche's aforementioned trick: */
3512  if (C_uhword_ref(v, n-1) == 0) n--;
3513
3514  /* These won't change during the loop, but are used in every step. */
3515  vn_1 = C_uhword_ref(v, n-1);
3516  vn_2 = C_uhword_ref(v, n-2);
3517
3518  /* See also Hacker's Delight, Figure 9-1.  This is almost exactly that. */
3519  for (j = m - n; j >= 0; j--) {
3520    hat = C_BIGNUM_DIGIT_COMBINE(C_uhword_ref(u, j+n), C_uhword_ref(u, j+n-1));
3521    if (hat == 0) {
3522      if (q != NULL) C_uhword_set(q, j, 0);
3523      continue;
3524    }
3525    qhat = hat / vn_1;
3526    rhat = hat % vn_1;
3527
3528    /* Two whiles is faster than one big check with an OR.  Thanks, Gauche! */
3529    while(qhat >= (1L << C_BIGNUM_HALF_DIGIT_LENGTH)) { qhat--; rhat += vn_1; }
3530    while(qhat * vn_2 > C_BIGNUM_DIGIT_COMBINE(rhat, C_uhword_ref(u, j+n-2))
3531          && rhat < (1L << C_BIGNUM_HALF_DIGIT_LENGTH)) {
3532      qhat--;
3533      rhat += vn_1;
3534    }
3535
3536    /* Multiply and subtract */
3537    k = 0;
3538    for (i = 0; i < n; i++) {
3539      p = qhat * C_uhword_ref(v, i);
3540      t = C_uhword_ref(u, i+j) - k - C_BIGNUM_DIGIT_LO_HALF(p);
3541      C_uhword_set(u, i+j, t);
3542      k = C_BIGNUM_DIGIT_HI_HALF(p) - (t >> C_BIGNUM_HALF_DIGIT_LENGTH);
3543    }
3544    t = C_uhword_ref(u,j+n) - k;
3545    C_uhword_set(u, j+n, t);
3546
3547    if (t < 0) {                /* Subtracted too much? */
3548      qhat--;
3549      k = 0;
3550      for (i = 0; i < n; i++) {
3551        t = (C_uword)C_uhword_ref(u, i+j) + C_uhword_ref(v, i) + k;
3552        C_uhword_set(u, i+j, t);
3553        k = t >> C_BIGNUM_HALF_DIGIT_LENGTH;
3554      }
3555      C_uhword_set(u, j+n, (C_uhword_ref(u, j+n) + k));
3556    }
3557    if (q != NULL) C_uhword_set(q, j, qhat);
3558  } /* end j */
3559}
Note: See TracBrowser for help on using the repository browser.